Toward fully anisotropic viscoelastic material models using an automated high-speed optical rig

This paper describes some of the algorithms that are employed in the postprocessing of the data obtained from an automated high-speed optical rig. The setup, entirely designed by some of the authors, is used for measuring the full field displacement in time of anisotropic viscoelastic solids undergoing controlled dynamic excitation. The setup main goal is to provide complex transfer functions relating the three-dimensional sinusoidal displacement on the surface of the sample to the sinusoidal motion of the foundation, suitable for the inverse estimation of the dynamic properties of the material under analysis. The complete automation of the measurement process aims at maximizing repeatability and signal-to-noise ratio. In particular, a stereo high-speed camera system is responsible for the actual measurement, while a laser Doppler vibrometer is used for a preliminary autoranging procedure and for data validation. The camera system does not directly output the high-quality data that is needed for a successful, physically meaningful inverse estimation. Hence, here the algorithms described in this paper come into play. Since digital image correlation (DIC) may fail in some frames or in some regions of the pictures, a mask is applied and calculated for every measurement run to preserve information only where the success rate of the DIC algorithm is higher than a specified threshold. The sample is mounted on a foundation whose sides are flat: random sample consensus (RANSAC) and principal component analysis (PCA) are used to respectively eliminate outliers and find the best fitting plane to the foundation side – the measured displacement field can then be expressed in a local coordinate system, truly relative to the foundation. The transfer functions are now obtained in two steps: first an optimizer is used to fit a cosine to all the measured displacement time-series, then the complex amplitudes representing the motion of the sample faces are divided by the average of the complex amplitudes representing the motion of the foundation.


Introduction
3D digital image correlation (3D-DIC) proved itself as a viable alternative to more conventional accelerometer-and laser-based measurements to validate finite element models and to provide qualitative information about the resonant behavior of simple structures [1,2,3,4,5]. Typical limitations of the method (e.g. the evaluation of displacements that the human eye cannot perceive when observing the image sequence) are being overcome with the integration of phasebased video magnification techniques [6]. Still, if the data obtained from 3D-DIC is meant to be the target of an inverse estimation process aimed at determining material parameters, it is necessary to develop a robust post-processing chain able to automatically deal with the mishaps 2 that inevitably occur when repeating a measurement a large number of times: failures of the DIC algorithm for some regions of the image or a misalignment of the test object with respect to the coordinate system used by the measurement equipment.
In this particular case, a cubic sample of a viscoelastic anisotropic material is to be dynamically characterized. The measurement of the displacement undergone by the outer faces of the sample, as well as the displacement of the driving foundation on which the sample is mounted, is computed using 3D-DIC and a pair of monochromatic high-speed cameras. The setup is briefly described in section 2, and yields three-dimensional full field displacement information in time for both the four sample side faces and the foundation. The data output by the 3D-DIC algorithm and in particular the flaws that render it not directly usable for a computation of the transfer functions are discussed in section 3. The solutions to such issues are discussed in sections 4 to 6 and represent the main contribution of this work. Results are shown for each algorithm, using both synthetic data for ease of understanding and real data to show the efficacy of the methods.

Brief setup description
The setup has been extensively described in [7], a brief overview is here given for the sake of self-containedness. Many details considered irrelevant for the current work are here omitted, the reader is again pointed at [7] for a much more detailed exposition. A cubic sample of the material to be tested is glued to a foundation with speckled sides, driven in the vertical direction by a uniaxial shaker imposing a sinusoidal motion to the foundation. A tuned seismic mass is glued to the top of the sample, and a couple of high-speed cameras record the motion of the sample and of the foundation for the whole duration of the experiment. The sample, along with the foundation and the seismic mass, is represented in figures 1 and 2.  The whole assembly (shaker, foundation, sample and seismic mass) is placed in a vacuum chamber on top of a rotation stage that allows for the rotation of the sample without opening the vacuum chamber. This avoids the need to re-calibrate the camera system each time the sample is rotated, since the pressure difference deforms the window through which the cameras look at the sample. Once the sample has been set up and the vacuum chamber depressurized, the measurement proceeds as follows: (i) the sample is rotated in such a way that both cameras see one of its faces; (ii) the shaker is driven with a pure sine of a randomly chosen frequency among the frequencies of interest; (iii) the amplitude is chosen by an iterative auto-ranging procedure (explained in [7]); (iv) once the system reaches steady state, the cameras acquire a stereo video recording of the sample and the foundation; (v) points (ii) to (iv) are repeated until all frequencies of interest are spanned; (vi) points (i) to (v) are now repeated for the three remaining faces of the sample.

Data obtained from the 3D-DIC camera system
The 3D-DIC is performed by a proprietary software bundled with the cameras. For each region of the image on which 3D-DIC is performed (two in our case, the foundation and the face of the sample), two outputs are obtained: • N files, one per stereo frame, containing the surface height map of the measured area; • N files, one per stereo frame, containing the coordinates of each measurement facet as computed in the reference image along with its three-dimensional displacement with respect to the aforementioned reference.
In order to compute the sought after transfer functions, it is necessary to know the amplitude and phase of each facet per each frequency and per each coordinate -assumed a linear elastic behavior of the sample and a three-dimensional cartesian coordinate system such that: • the y axis is vertical, parallel to the direction in which the shaker excites the foundation; • the face of the foundation currently observed entirely lies in the xy plane; • the origin of such coordinate system is not important, since zero-mean signals are of interest for the current application.
Sadly, the DIC algorithm does not necessarily return valid results for all frames nor for all measurement facets. This can happen between the frames in a stereo pair (e.g. specular reflections, appearing as extremely bright points in only one of the two frames) as well as between stereo frames captured at different times (a facet laying on the boundary of the sample would be particularly tricky to match). To automatically obtain the transfer functions it is then necessary to exclude the facets that are not consistently tracked. A procedure for doing so is described in section 4.
Due to the need of measuring in vacuum, a thick sheet of acrylic glass separates the cameras from the sample. The calibration routine adopted by the software forces the coordinate system to be aligned with the calibration plate [8,9]. This implies that if the face of the sample under observation is not perfectly parallel to the calibration plate (which is the real scenario, given the unlikeliness of achieving a perfect alignment) then the displacement fields obtained need to be expressed in a new coordinate system before they can be used to compute transfer functions. This scenario is dealt with in section 5.
Once the valid displacement time series have been transformed in the new coordinate system, the amplitude and phase of each facet per each frequency has to be evaluated. This is done by fitting a cosine to each displacement time series, as described in section 6.

Discarding invalid measurement facets
The algorithm here described is operated on the data output of each individual video recording, that is for each individual frequency and for each face of the sample and of the foundation. A description of the algorithm in natural language follows, an example from real measurement data being presented in figure 3. The choice of a rectangular domain is done for mere convenience: on such a regular domain it is easier to perform a variety of operations (e.g. a one-and two-dimensional FFT), should it be needed in the future. Moreover, as the size of the facets is reduced, the shape of the rectangular domain better approximates the shape of the sample.

Changing the coordinate system
As it can be observed in figure 3, the origin of the coordinate system does not overlap with the physical center of the sample face. While this would not be a problem in itself, it is a symptom of the misalignment of the sample and the calibration plate. The foundation on the other hand (painted with a speckle pattern and clearly visible in figure 2) is aligned with the sample and provides for a useful reference to express the displacement data in a more convenient measurement system. In particular the foundation height map is a point cloud that should ideally fit into a plane: such a plane would provide the xy-plane of the new coordinate system, relative to the foundation and with the y-axis aligned with the direction of excitation. In reality, the foundation height map is not a perfect plane because of two different phenomena: (i) the resolving power of the cameras is limited (as it is for all digital equipment), hence the true location of the points oscillates around the true value; (ii) outliers, due to imperfect manufacturing of the foundation and errors in the matching algorithm, pollute the point cloud. Keeping these phenomena in mind, it is necessary to remove the outliers before using the point cloud to determine the best fitting plane hence deciding on a coordinate system.
The outlier detection is operated using the random sample consesus (RANSAC) paradigm, an iterative method well suited to estimate model parameters from data that is contaminated with outliers [10]. RANSAC requires the setting of some confidence parameters, of which the most critical is the threshold value for determining whether a data point fits the model or not. The determination of this threshold value is the only operation that requires some human intervention, since it is sensitive to the surface roughness of the material observed and to the quality of the camera calibration. After RANSAC has determined the set of inliers, principal component analysis (PCA) can be used to determine the eigenvectors that will transform the displacement fields to the new coordinate system [11]. The procedure is illustrated using synthetic data for the sake of simplicity in figure 4.

From time series to amplitudes and phases
There are many different ways to extract an amplitude and a phase from an ideally sinusiodal signal, with the Discrete Fourier Transform being among the most common tools in the arsenal. Still, the main drawback of such a method is that (all the other parameters being fixed) high frequency resolution requires a high number of samples [12], a computationally expensive luxury when dealing with 3D-DIC. The approach used relies instead on solving the nonlinear leastsquares problem of fitting a cosine to the time series, a well-posed problem given reasonable bounds and initial guesses [13]. A one-dimensional linear interpolation is used to recover information for missing data points, where the DIC might have failed. Then amplitude, phase, offset and frequency are the four parameters that are obtained for each uni-directional displacement of each facet. An example with real measurement data is provided in figure 5. While such an approach might mask any non-linear behavior of the system, the material model used in subsequent stages for the estimation of the dynamic properties is a linear one: as such, it would be unable to incorporate nonlinear phenomena in its current state. The estimated covariance matrices of the optimal values resulting from the solution of the nonlinear least-squares are nevertheless an indicator of the quality of the fit. Foundation vertical displacement at 50.03 Hz, tile (0,0) discrete data from DIC initial guess best fitting cosine Figure 5: Vertical displacement of one of the foundation facets: the DIC data from 50 stereo frames, the initial guess and the best fitting cosine are shown. The method works even with non-periodic data, and in the case of frequencies of longer time period than the duration of the observation.

Conclusions
Part of a much broader project, the work presented in this paper allows for the fully automatic post-processing of the data produced by a 3D-DIC system meant for the characterization of viscoelastic anisotropic materials: the outcome of the process is a series of three-dimensional displacement transfer functions. Future work may include (but not be limited to): • uncertainty analysis of the measured quantities; • an automatic, smart choice of the RANSAC threshold value (e.g. MLESAC [14]); • a high-performance, highly parallel implementation of the above routines; • development of an own 3D-DIC solver.