Three-Dimensional Digital Image Correlation Based on Speckle Pattern Projection for Non-Invasive Vibrational Analysis

Non-contact vibration measurements are relevant for non-invasively characterizing the mechanical behavior of structures. This paper presents a novel methodology for full-field vibrational analysis at high frequencies using the three-dimensional digital image correlation technique combined with the projection of a speckle pattern. The method includes stereo calibration and image processing routines for accurate three-dimensional data acquisition. Quantitative analysis allows the extraction of several deformation parameters, such as the cross-correlation coefficients, shape and intensity, as well as the out-of-plane displacement fields and mode shapes. The potential of the methodology is demonstrated on an Unmanned Aerial Vehicle wing made of composite material, followed by experimental validation with reference accelerometers. The results obtained with the projected three-dimensional digital image correlation show a percentage of error below 5% compared with the measures of accelerometers, achieving, therefore, high sensitivity to detect the dynamic modes in structures made of composite material.


Introduction
Non-contact optical sensing technologies combined with data processing algorithms play an essential role in non-invasively assessing a structure's dynamic behavior. Structures with non-uniform material properties and/or complex geometry exhibit spatially local structural behaviors and are temporally transient. Therefore, full-field vibration-based measurements at very high spatial resolutions are required for the characterization and analysis of the structural dynamics for accurate identification of the modal parameters (i.e., modal shapes, frequencies and damping ratios) [1,2].
In recent years, optical-based vibration measurements have been proposed to measure structural dynamics, as they overcome the limitations of contact sensors (e.g., accelerometers), provide high-resolution measurements and do not change the structural dynamic behavior during testing. There exist a number of optical techniques for analyzing the dynamic response of a structure, such as laser Doppler vibrometry [3], electronic speckle interferometry [4,5], digital speckle shearography [6,7] and digital image correlation [8][9][10][11], among others. Laser Doppler vibrometry is the most extensively used technique for high spatial resolution vibration measurements owing to its superior resolution, sensitivity, robustness and non-contact nature, but it is an expensive and time-consuming method, which only offers pointwise measurements. Alternatively, scanning and continuous scanning laser Doppler vibrometry integrate a scanning system for full-field vibration measurements; however, different authors described (1) pseudo-vibration due to surface roughness and (2) measurement errors due to defocusing and motion artifacts at some scanning points [12,13].
This study proposes a custom-developed 3D-DIC methodology based on speckle pattern projection and image processing routines for full-field vibrational analysis at high frequency. The potential of the technique has been demonstrated on an Unmanned Aerial Vehicle (UAV) wing made of composite material, followed by experimental validation with reference accelerometers. The vibrational parameters, such as 3D displacement fields, natural frequencies and modal shapes, have been obtained and studied to demonstrate the capability of the developed methodology. Figure 1 illustrates the custom experimental set-up. Basically, the system comprises an optical channel for speckle pattern projection and two high-speed sensors. The optical channel consists of an on-axis illumination module (Broadband halogen fiber optic illuminator, OSL2, Thorlabs, Bergkirchen, Germany), which emits a stable power in the spectral range between 185 and 2000 nm, the optical elements required to collimate the beam and generate a uniform irradiance across the sample include two converging lenses and a printed speckle pattern in acetate. The off-axis imaging channel consists of two highspeed and high-resolution cameras (Mini AX100, Photron, Bucks, UK; 1024 × 1024 pixels at 4000 fps; pixel size: 20 µm; accuracy: 0.01 pixel) provided with a 150 mm focal length objective lens (Irix 150 mm f/2.8, Irixlens, Switzerland) and mounted on the X-Y-axis translation stage with a high precision rotation of 360 degrees (PR01, Thorlabs, Bergkirchen, Germany). iterative algorithm to process the information [44], the sub-pixel interpolation scheme [45,46], the quality and resolution of images [47], the speckle pattern quality [48] and the subset size [24]. This study proposes a custom-developed 3D-DIC methodology based on speckle pattern projection and image processing routines for full-field vibrational analysis at high frequency. The potential of the technique has been demonstrated on an Unmanned Aerial Vehicle (UAV) wing made of composite material, followed by experimental validation with reference accelerometers. The vibrational parameters, such as 3D displacement fields, natural frequencies and modal shapes, have been obtained and studied to demonstrate the capability of the developed methodology. Figure 1 illustrates the custom experimental set-up. Basically, the system comprises an optical channel for speckle pattern projection and two high-speed sensors. The optical channel consists of an on-axis illumination module (Broadband halogen fiber optic illuminator, OSL2, Thorlabs, Bergkirchen, Germany), which emits a stable power in the spectral range between 185 and 2000 nm, the optical elements required to collimate the beam and generate a uniform irradiance across the sample include two converging lenses and a printed speckle pattern in acetate. The off-axis imaging channel consists of two high-speed and high-resolution cameras (Mini AX100, Photron, Bucks, UK; 1024 × 1024 pixels at 4000 fps; pixel size: 20 µ m; accuracy: 0.01 pixel) provided with a 150 mm focal length objective lens (Irix 150 mm f/2.8, Irixlens, Switzerland) and mounted on the X-Y-axis translation stage with a high precision rotation of 360 degrees (PR01, Thorlabs, Bergkirchen, Germany).

Projected 3D-DIC: Fundamentals
For high-frequency vibration analysis in projected 3D-DIC, the grey intensity of the reference image (unloaded state) is compared with the high-speed acquired images under deformation (loaded state). However, DIC is not directly applied between the reference Figure 1. Experimental custom 3D-DIC set-up based on speckle pattern projection for full-field vibration measurements: (a) schematic layout of the experimental set-up and testing conditions (UAV wing made of composite material (test sample) clamped to the shaker); (b) set-up employed in the vibration test; (c) image of the UAV wing fixed to the electromagnetic shaker (lateral excitation); (d) optical channel for speckle pattern projection and high-speed sensors (distance of the tested sample during the measurements: 3.75 m; instantaneous field of view of the sensor: 0.5 mm).

Projected 3D-DIC: Fundamentals
For high-frequency vibration analysis in projected 3D-DIC, the grey intensity of the reference image (unloaded state) is compared with the high-speed acquired images under deformation (loaded state). However, DIC is not directly applied between the reference and the deformed images; the acquired images are processed and divided into subsets ( Figure 2).  Figure 3 describes the developed procedure for structural characterization using projected 3D-DIC. First, before obtaining the calibration images, the high-speed sensors, the optical channel for speckle pattern projection and the sample should be perfectly aligned and uniformly illuminated for 3D-DIC applications. The high-precision rotation stage ensures an accurate orientation of the cameras at 30 degrees in relation to the sample for stereo imaging. This configuration resulted in a scale of 0.5 mm/pixel. Once the speckle pattern is focused, the image acquisition can be performed on unloaded and loaded state images. Then, image processing algorithms and routines were developed to obtain the 3D field displacement and include: (i) the calibration and distortion parameters estimation, (ii) the alignment algorithm to obtain the displacement between subsets (2D-DIC) and (iii) the 3D reconstruction.  Figure 3 describes the developed procedure for structural characterization using projected 3D-DIC. First, before obtaining the calibration images, the high-speed sensors, the optical channel for speckle pattern projection and the sample should be perfectly aligned and uniformly illuminated for 3D-DIC applications. The high-precision rotation stage ensures an accurate orientation of the cameras at 30 degrees in relation to the sample for stereo imaging. This configuration resulted in a scale of 0.5 mm/pixel. Once the speckle pattern is focused, the image acquisition can be performed on unloaded and loaded state images. Then, image processing algorithms and routines were developed to obtain the 3D field displacement and include: (i) the calibration and distortion parameters estimation, (ii) the alignment algorithm to obtain the displacement between subsets (2D-DIC) and (iii) the 3D reconstruction.  Figure 3 describes the developed procedure for structural characterization using projected 3D-DIC. First, before obtaining the calibration images, the high-speed sensors, the optical channel for speckle pattern projection and the sample should be perfectly aligned and uniformly illuminated for 3D-DIC applications. The high-precision rotation stage ensures an accurate orientation of the cameras at 30 degrees in relation to the sample for stereo imaging. This configuration resulted in a scale of 0.5 mm/pixel. Once the speckle pattern is focused, the image acquisition can be performed on unloaded and loaded state images. Then, image processing algorithms and routines were developed to obtain the 3D field displacement and include: (i) the calibration and distortion parameters estimation, (ii) the alignment algorithm to obtain the displacement between subsets (2D-DIC) and (iii) the 3D reconstruction.

Camera Calibration Parameters
The calibration of the high-speed cameras was performed in two different steps: (i) data acquisition for the calibration; and (ii) calculation of the 3D coordinates through an iterative process that requires the estimation of intrinsic (the optical characteristics of the sensor) and extrinsic (position and orientation of the sensor in the set-up) parameters. Figure 4 shows the scheme of the coordinate system. To calibrate a single camera, the method developed by Zhang [49] was implemented, where a set of images of a chessboard with different positions and orientations covering the largest possible area of the analysis region was required. The relationship between 3D points and planar points is defined by Equation (1), using Zhang's notation.
. In this matrix, 0 and 0 are the coordinates of the principal point resulting from the intersection of the image plane with the optical axis, which is represented by in Figure  4. Coefficients and are the focal length and is the skew coefficient. The extrinsic parameter set is constituted by rotation matrix and the translation vector between the world coordinate system and the camera system. Finally, is an arbitrary scale factor.  The following items describe the process of solving the single-camera calibration problem: 1. A set of thirty images of the chessboard with different orientations was recorded. The coordinate of each corner was obtained ( [ 1] ) and it was associated with the three-dimensional point defined by [ 1 ] ; 2. The analytical solution of Equation (1) was calculated using Zhang's method; 3. A nonlinear optimization based on the maximum likelihood criterion was developed.
The procedure followed is explained in detail in Appendix A.  A 3D point is defined by X Y Z 1 T and the coordinates in the image plane are given by u v 1 T . The matrix constituted by the intrinsic parameter set is denoted by A. In this matrix, u 0 and v 0 are the coordinates of the principal point resulting from the intersection of the image plane with the optical axis, which is represented by C in Figure 4. Coefficients α and β are the focal length and γ is the skew coefficient. The extrinsic parameter set is constituted by rotation matrix R and the translation vector T between the world coordinate system and the camera system. Finally, λ is an arbitrary scale factor.
The following items describe the process of solving the single-camera calibration problem:

1.
A set of thirty images of the chessboard with different orientations was recorded. The coordinate of each corner was obtained ( u v 1 T ) and it was associated with the three-dimensional point defined by X Y Z 1 T ;

2.
The analytical solution of Equation (1) was calculated using Zhang's method;

3.
A nonlinear optimization based on the maximum likelihood criterion was developed. The procedure followed is explained in detail in Appendix A.
The two high-speed cameras were synchronized to obtain the images simultaneously for stereoscopic calibration and measurements. Therefore, the rotation matrix R and the translation vector T between the two cameras used for that pairs of calibration images were obtained. Each pair of images correspond with the same chessboard scene. The nonlinear iterative procedure for the stereoscopic calibration is described in Appendix B.

Distortion Correction
Nonlinear distortion is inherent to imaging systems: radial distortion [50] and tangential distortion [51]. The undistorted image coordinates u v are expressed by Equations (2) and (3), where x and y are the coordinates with respect to the principal point of image plane (C in Figure 4), k i is the coefficient of distortion and r is given by r = x 2 + y 2 .
Both Equations (2) and (3) have been integrated into the calibration estimation procedure. In particular, in the nonlinear optimization of single calibration and in the nonlinear optimization employed in the stereoscopic calibration.

Alignment Algorithm
The developed routine to obtain the displacement between two subsets is based on the inverse compositional Gauss-Newton algorithm (IC-GN) [15,52], and it consists of the alignment of a template (deformed) image to a target image. For that, the grey-level intensity of each reference subset is correlated to the grey-level intensity of the corresponding template subset. The mathematical expression of the IC-GN algorithm performed on a zero-mean normalized sum of squared difference (ZNSSD) correlation criteria is shown in Equations (4) to (11); in which functions F and G represent the subset of reference images and the subset of deformed images, respectively. The term H corresponds to the hessian calculated by means of Equation (8), variable p is the pre-computed deformation parameter vector, variable ξ designs the local coordinates of pixel point in each subset and W is the warp function given by Equation (9).
The working principle of the IC-GN algorithm is represented in Figure 5. The initial estimation of p = u 0 0 v 0 0 is usually calculated by means of the direct application of the ZNSSD criteria or another equivalent criteria, such as zero normalized crosscorrelation (ZNCC). This allows us to obtain an initial estimation of the u and v components of the displacement vector between the reference subset and the target subset. After that, an iterative loop is performed until achieving the convergence criteria (represented in Figure 5).

Three-Dimensional Reconstruction
The three-dimensional reconstruction or triangulation is the reconstruction of a 3D point from its 2D projections from two or more cameras. The least square method (LSM) [53] was implemented because of its accuracy and computational efficiency. The relationships between the left image coordinate u l v l , the right coordinate u r v r and the world coordinate system X Y Z are indicated in Equations (12) and (13), respectively [53]. Integrating both equations, the world coordinate system can be expressed by Equation (14), where b is given by (15) and M + is the left-pseudo inverse (M + = M T M −1 M T ) of M matrix (16).

Projected 3D-DIC: Test and Validation
The structure under testing was the horizontal tail plain (HTP) of the DIANA Unmanned Aerial Vehicle (UAV) with a mass of 530 g and the following dimensions: 455 mm (long) × 265 mm (wide). The DIANA UAV was manufactured with HexPly ® M56/35%UD194/ IMA12K prepreg, and it was composed of a main skin made of a symmetric laminate [60,-60,0], a longitudinal stringer and reinforcements in the borders made of a symmetric laminate [60,-60,03]. All mechanical joints were made with adhesive. For validation, three accelerometers were placed in the structure (mass of the accelerometer: 8.6 g). Their positions are described in Table 1 (referenced to the upper left vertex) and illustrated in Figure 6. In this study, the type of accelerometer employed was a triaxial sensor with a measuring range of 50 g, a typical sensitivity of 100 Mv/g and an analysis range of up to 10 kHz. For the series of test, the UAV wing was installed on a slipping table fixed to an electromagnetic shaker (LDS 406) powered by a PA100E power amplifier. The drive signal, as well as the acquisition of the accelerometer signals, was performed by means of an LDS Focus II acquisition system. First, a modal identification survey to identify the natural frequencies of the structure was performed; and then the base of the structure was displaced sinusoidally at each natural frequency in order to validate the projected 3D-DIC with the accelerometers.  fixed to an electromagnetic shaker (LDS 406) powered by a PA100E power amplifier. The drive signal, as well as the acquisition of the accelerometer signals, was performed by means of an LDS Focus II acquisition system. First, a modal identification survey to identify the natural frequencies of the structure was performed; and then the base of the structure was displaced sinusoidally at each natural frequency in order to validate the projected 3D-DIC with the accelerometers.

Results
In order to evaluate the potential of the projected 3D-DIC for long-duration random vibration measurements, a series of trials were conducted focusing on two major results: (i) modal characterization, which means the modal identification survey to obtain the natural frequencies of the structure inherent to the vibration testing; and (ii) validation of the method with excitation of the UAV wing at the detected natural frequencies of the structure. The 3D-DIC results were compared to the accelerometers' measurements.

Modal Characterization
For modal determination, a frequency sweep from 0 to 200 Hz was performed to obtain the most relevant natural frequencies of the structure (accelerometers) by programming a sine test at a rate of 1/8 min with a frequency resolution that increases linearly with the frequency sweep. The accelerometers show a resonance frequency of 10 kHz in order to ensure high accuracy in the tested frequency range. The most relevant natural

Results
In order to evaluate the potential of the projected 3D-DIC for long-duration random vibration measurements, a series of trials were conducted focusing on two major results: (i) modal characterization, which means the modal identification survey to obtain the natural frequencies of the structure inherent to the vibration testing; and (ii) validation of the method with excitation of the UAV wing at the detected natural frequencies of the structure. The 3D-DIC results were compared to the accelerometers' measurements.

Modal Characterization
For modal determination, a frequency sweep from 0 to 200 Hz was performed to obtain the most relevant natural frequencies of the structure (accelerometers) by programming a sine test at a rate of 1/8 min with a frequency resolution that increases linearly with the frequency sweep. The accelerometers show a resonance frequency of 10 kHz in order to ensure high accuracy in the tested frequency range. The most relevant natural frequencies obtained from the UAV wing were: 42.8 Hz (bending), 57.5 Hz (twisting), 83.1 Hz (twisting) and 97.5 Hz (twisting). Hence, the high-speed cameras were programmed for high-frequency sampling (1000 Hz) at maximum resolution (1024 × 1024 pixels). Three different subset sizes were considered: 20 × 20 pixels, 30 × 30 pixels and 40 × 40 pixels. Figure 7 shows the cross-correlation coefficients for each subset size. It should be noted that as the size of the subset decreases, the probability that the cross-correlation provides erroneous displacements increases mainly due to two reasons: (i) decreasing the subset size increases the magnitude of the residual correlation coefficients, which increases the probability of obtaining wrong displacements; (ii) the smaller the size of the subset, the greater the probability that any of the subsets is empty, with a low number of speckle points in its interior. In the case of applying the cross-correlation to an empty subset, the displacement obtained would be a random number. For the subset size of 20 × 20 pixels, the cross-correlation residues reached values comparable to the maximum value of crosscorrelation coefficients. The subset sizes of 30 × 30 and 40 × 40 pixels showed reliable results, so for further analysis, the configuration with the smaller subset size (30 × 30 pixels) was considered. points in its interior. In the case of applying the cross-correlation to an empty subset, the displacement obtained would be a random number. For the subset size of 20 × 20 pixels, the cross-correlation residues reached values comparable to the maximum value of crosscorrelation coefficients. The subset sizes of 30 × 30 and 40 × 40 pixels showed reliable results, so for further analysis, the configuration with the smaller subset size (30 × 30 pixels) was considered.  Figure 8 shows the amplitude of the displacements for each natural frequency. It should be mentioned that the UAV wing was excited to each of these natural frequencies by means of a sine load of the natural frequency. The amplitude of the displacement was calculated using expression (23): = √ 2 + 2 + 2 (23) where , and represent the displacement of each component and represents the total displacement. Moreover, the sign of displacement was taken into account. The parts of the wing with positive displacement represent an area of the UAV wing that has deformed in the opposite direction to areas with negative displacement. The parts of the wing with null displacement represent a wing area that has not been deformed. Figure 8, Videos S1 and S2 show the raw data (unprocessed images of camera 1 and camera 2) and the temporal evolution of the displacements for the first natural frequency (Figure 9 shows a screenshot of each visualization, the videos can be found as supplementary material).  Figure 8 shows the amplitude of the displacements for each natural frequency. It should be mentioned that the UAV wing was excited to each of these natural frequencies by means of a sine load of the natural frequency. The amplitude of the displacement was calculated using expression (23): where d X , d Y and d Z represent the displacement of each component and D represents the total displacement. Moreover, the sign of displacement was taken into account. The parts of the wing with positive displacement represent an area of the UAV wing that has deformed in the opposite direction to areas with negative displacement. The parts of the wing with null displacement represent a wing area that has not been deformed. Figure 8, Videos S1 and S2 show the raw data (unprocessed images of camera 1 and camera 2) and the temporal evolution of the displacements for the first natural frequency (Figure 9 shows a screenshot of each visualization, the videos can be found as supplementary material).     (a) Raw data acquisition (Video S1) (b) Displacement field (Video S2) Figure 9. Screenshots of supplementary videos: (a) raw data acquisition (Video S1) and (b) the displacement field corresponding to the post-processing of raw data from two different perspective points (Video S2). The color bars represent the displacement field (mm).

Experimental Validation
The displacements obtained with the projected 3D-DIC for the first vibration mode (42.8 Hz) were compared with the results obtained with the accelerometers. The values of the accelerometers were normalized with respect to accelerometer 1, the ratio between accelerometers 2 and 1 was 2.75 and the ratio between accelerometers 3 and 1 was 3.25. The displacements obtained with the projected 3D-DIC were also normalized with respect to displacement corresponding to coordinates of accelerometer 1. For that, an offset is applied to all displacements in order to fix the displacement of accelerometer 1 to 1 pixel. After that, the displacements were contrasted with respect to the ratios of accelerometers 2 and 3. Assuming that the accelerometers were not perfectly placed in the coordinates indicated in Table 1, a sweep was performed from −1 subset to +1 subset around the position of both accelerometers (accelerometers 2 and 3). The result of this operation is shown in Figure 10. For Figure 10a, the errors ranged from 1.8% (center) to 3.38 % (bottom left), while for Figure 10b, the error ranged from 3.7% (center) to 6.3% (bottom left).
plied to all displacements in order to fix the displacement of accelerometer 1 to 1 pixel. After that, the displacements were contrasted with respect to the ratios of accelerometers 2 and 3. Assuming that the accelerometers were not perfectly placed in the coordinates indicated in Table 1, a sweep was performed from −1 subset to +1 subset around the position of both accelerometers (accelerometers 2 and 3). The result of this operation is shown in Figure 10. For Figure 10a, the errors ranged from 1.8% (center) to 3.38 % (bottom left), while for Figure 10b, the error ranged from 3.7% (center) to 6.3% (bottom left).
(a) Percentage of error corresponding to accelerometer 2 area (b) Percentage of error corresponding to accelerometer 3 area

Discussion
This study proposed an accurate method for full-field vibration measurement. The method is based on projected 3D-DIC and is capable of non-invasively identifying the full-field mode shapes and natural frequencies of a UAV wing made of composite material, which could potentially benefit structural dynamics and health monitoring applications [54]. Previous studies described fringe patterns [55,56], laser speckle [41,57] or a combination of both [40] to create a non-invasive pattern for full-field vibration measurement in DIC applications. In a recent study, Felipe-Sesé et al. [40] proposed a combination of a fringe pattern, a laser speckle and a single-camera DIC to measure the three-dimensional shapes and the associated field of displacements of a plate excited by a shaker. The authors validated their results with theoretical predictions and obtained a maximum error near the edges of the structure of 4.95%. Moreover, Pang et al. [58], evaluated the flexural behavior of concrete sleepers with a laser speckle imaging sensor. The authors vali-dated the methodology with foil strain gauges and described a maximum error of 7.15% in their measurements.
In this study, a speckle pattern projection has been chosen not only as an alternative to the conventional invasive painted speckle patterns but also to other non-invasive and complex options (i.e., laser speckle or fringe patterns). Overall, the study holds similarities with those recently reported in the state-of-the-art and includes the novelty of the analysis on a composite material at different vibration modes, showing high sensitivity to detect the temporal evolution of the vibration modes. Accurate 3D mode shapes were obtained for different resonances, the results show a percentage of error below 3% and 5% compared with reference accelerometers 2 and 3, respectively. Accelerometer 3 was placed near the edge of the structure, and the discontinuity at this region might be the reason for a reduced accuracy [59,60].
Although the measurements were performed for a sampling frequency of 1000 Hz, one of the limitations of the study is that the more relevant vibration modes of the UAV structure occurred in the range from 0 to 100 Hz. Another limitation is associated with the out-of-plane motion tracking since the projected speckle pattern does not move tightly together with the surface of the structure. However, the latter might not be considered a weakness for vibrational analysis because the excitation of the structure is applied longitudinally on each axis, and the whole out-of-plane displacement measurements of the measured area are fundamental for the evaluation of the vibration modes of the structure. The advantages of the proposed methodology are the high-frequency capacity (up to 4000 Hz under the same configuration) and great versatility to (i) customize the speckle projected pattern for canceling the residues of the cross-correlation coefficients and (ii) configurate the set-up for measurements in a selected spectrum range (i.e., near-infrared). This methodology could be directly applied for vibrational analysis of optical surfaces, such as solar cells, glass structures and lenses, or even for biomechanical analysis in human tissues (e.g., transparent tissues such as the cornea).

Conclusions
3D-DIC based on speckle pattern projection is an accurate method to characterize the structural dynamics of a composite material. The proposed method successfully obtained full-field vibration measurement with high spatial resolution and is capable of reliably extracting the shape and intensity, as well as the out-of-plane displacement fields and mode shapes of the structure.

Patents
where, • A is the matrix of the intrinsic parameters constituted by α, β, γ, u 0 and v 0 according to Zhang's nomenclature [49]. • k represents a vector with radial and tangential distortion parameters. • R j is the rotation matrix, which corresponds to image j.
• T j is the translation vector, which corresponds to image j. • P(w) l is the vector constituted by w parameters, which corresponds to iteration l.  where, • is the matrix of the intrinsic parameters constituted by , , , 0 and 0 according to Zhang's nomenclature [49].
• represents a vector with radial and tangential distortion parameters. • is the rotation matrix, which corresponds to image . • is the translation vector, which corresponds to image . • ( ) is the vector constituted by parameters, which corresponds to iteration . • ( ) is a vector with the error estimation of each corner. • ( ) −1 is the inverse matrix of the hessian matrix defined in Equation (A10). • Ƞ is the smoothing factor. • , is the two-dimensional corner coordinates of the chessboard, they must be obtained with Harris corner detector.

Appendix B
The stereoscopic calibration parameters are constituted by the rotation matrix from the left camera system to the right camera system (R r←l ), which defines the orientation between both cameras and the translation vector between the left camera system to the right camera system (T r←l ). Just as the estimation of the single camera parameters, the solution was calculated by minimizing the difference between the two-dimensional known coordinates (Harris corner detector) and the two-dimensional estimated coordinates of the corners. However, it is important to highlight that before estimating the stereoscopic parameters, the parameters of each camera should be calculated. As shown in Equation (A11), the number of parameters has grown significantly by incorporating the left single-camera parameters (denoted with superscript l), the right single-camera parameters (defined with superscript r) and the estereoscopic parameters (without superscript). The iterative refinement starts with an initial estimation of ρ parameters. The majority of the parameters were previously calculated by means of single-camera calibration and the initial estimation of R r←l and T r←l were calculated according to the mean values obtained via Equations (A13) and (A14). ρ = A l , A r , k l , k r , R r j , R l j , R r←l , T r j , T l j , T r←l (A11) The iterative loop (see Figure A2) was defined by the following steps: • Estimation of the projected pixel coordinates corresponding to the corners of the images (left camera) by using Equation (A4). • Estimation of the rotation matrix from the world system to the right camera system (R r←w ) by using Equation (A13) and the translation vector from the world system to the right camera system (T r←l ) by means of Equation (A14).

•
Estimation of the projected pixel coordinates corresponding to the corners of the images (right camera) by using Equation (A4).

•
Calculation of the vector error as the difference between the two-dimensional known coordinates and the estimated coordinates of the corners of all the images (both cameras, left and right).