Stereo Digital Image Correlation in MATLAB

: Digital Image Correlation (DIC) has found widespread use in measuring full-ﬁeld displacements and deformations experienced by a body from images captured of it. Stereo-DIC has received signiﬁcantly more attention than two-dimensional (2D) DIC since it can account for out-of-plane displacements. Although many aspects of Stereo-DIC that are shared in common with 2D DIC are well documented, there is a lack of resources that cover the theory of Stereo-DIC. Furthermore, publications which do detail aspects of the theory do not detail its implementation in practice. This literature gap makes it difﬁcult for newcomers to the ﬁeld of DIC to gain a deep understanding of the Stereo-DIC process, although this knowledge is necessary to contribute to the development of the ﬁeld by either furthering its capabilities or adapting it for novel applications. This gap in literature acts as a barrier thereby limiting the development rate of Stereo-DIC. This paper attempts to address this by presenting the theory of a subset-based Stereo-DIC framework that is predominantly consistent with the current state-of-the-art. The framework is implemented in practice as a 202 line MATLAB code. Validation of the framework shows that it performs on par with well-established Stereo-DIC algorithms, indicating it is sufﬁciently reliable for practical use. Although the framework is designed to serve as an educational resource, its modularity and validation make it attractive as a means to further the capabilities of DIC.


Introduction
Digital Image Correlation (DIC) is an optical metrology technique capable of determining full-field displacements and deformations experienced by a body from images captured of the body's surface.The accuracy of DIC results is greatly influenced by the quality of the speckle pattern [1], the internals of the DIC algorithm [2,3], the appropriateness of the chosen processing parameters [4] and the quality of the calibration process [5][6][7], among others.However, despite this, DIC, a non-interferometric technique, is more robust in withstanding vibrations and ambient light variations than interferometric techniques, such as Electronic Speckle Pattern Interferometry and Moiré interferometry, allowing it to be used outside the confines of a laboratory.Additionally, its non-contact, full-field nature; its ability to function at a range of scales (from scanning electron microscopy images [8,9] to satellite images [10]); and the rapid improvement of cameras has resulted in it becoming one of the most favored techniques in experimental solid mechanics.
Following the introduction of two-dimensional (2D) DIC in 1982 by Peters and Ranson [38] it was extended to three dimensions (Stereo-DIC) by Luo et al. [39] in 1993 and to digital volume correlation (DVC) by Bay et al. [40] in 1999.The correlation aspect of DIC has undergone significant refinement over the years.Current state-of-the-art DIC algorithms as identified by Pan [41] make use of: (i) bi-quintic b-spline interpolation which was shown by Schreier et al. [3] to produce the most accurate sub-pixel displacements; (ii) the first and second-order shape functions (SFs) proposed by Peters et al. [42] and Lu and Carry [43] respectively; (iii) the zero-mean normalized sum of squared difference (ZNSSD) correlation criterion suggested by Tong [44]; (iv) Gaussian pre-filtering of images to reduce bias in the displacements caused by high frequency noise in the images as proposed by Pan [45]; and (v) the inverse compositional Gauss-Newton (IC-GN) optimization method introduced by Baker and Matthews [46].Although the IC-GN optimization method is theoretically equivalent to the Newton-Raphson method [46], in practice it has greater efficiency, accuracy and robustness to withstand noise [47] making it preferable.
Although DIC is a complex process, the papers proposed by Pan et al. [48], Gao et al. [49], Solav et al. [50], and Blaber et al. [51] provide an in-depth explanation of the theory of the DIC process, which coupled with the Good Practices Guide for DIC [52], enables newcomers to the field of DIC to gain the necessary knowledge to effectively apply it in established use cases.However, in order for a newcomer to contribute to the field of DIC by developing new algorithms which improve its capabilities or adapting it for new applications, he or she must gain a deep understanding of DIC.Gaining such a deep understanding is cumbersome due to the lack of resources which directly bridge the gap between the theory of DIC and its coded implementation.More specifically, although some publications release code that is consistent with the theory presented [50,51], these codes focus on robustness and ease of use which, despite making them more suitable for real world applications, results in a complex code which is ineffective as a learning resource.This lack of resources acts as a barrier to newcomers intending to further the capabilities of DIC, thereby limiting the development rate of the field.
It is for this reason that the authors published a paper bridging the gap between the theory of a 2D DIC framework (ADIC2D) and its practical implementation as a modular 117 line MATLAB code [53].However, 2D DIC, being limited to determining in-plane displacements, is susceptible to errors in the presence of out-of-plane motion, which is generally unavoidable especially outside the confines of a laboratory [54].As such, Stereo-DIC has become favored since it accounts for out-of-plane motion, allowing for more reliable use across a broader range of applications and fields.This is supported by the significantly larger number of publications of Stereo-DIC relative to 2D DIC.
This paper builds upon the work of the preceding paper by extending ADIC2D to Stereo-DIC with the aim of bridging the gap between the theory and implementation of Stereo-DIC.This is done by first presenting the theory of a subset-based, Stereo-DIC framework that is predominantly consistent with current state-of-the-art techniques.How this theory is implemented as a modular 202 line MATLAB code is then discussed.These theory and implementation sections focus on aspects of Stereo-DIC which differ from that of ADIC2D.Thereafter the framework is validated using image sets provided by the Stereo-DIC Challenge [55].More specifically, ADIC3D's results are compared to those of LaVision's StrainMaster software and DICe [56] in order to make observations about its capabilities.Finally, the proposed framework's modularity, limitations, and overall performance are discussed.
The framework, referred to as ADIC3D, was developed in MATLAB since its simple syntax lets the reader focus on the core mathematics of the code.Furthermore, MATLAB's efficient built-in functions are leveraged to simplify the code and reduce its computation time.The code is designed to be modular, so that the link between the theory and code is clear and enables readers to systematically develop an understanding of the code.Additionally, this modularity allows the code to be easily altered, urging readers to further the capabilities of DIC.

Framework Theory
Stereo-DIC is fundamentally an extension of 2D DIC, which relies upon two or more cameras simultaneously capturing images of the specimen (where specimen refers to the body being tracked) from different perspectives in order to determine its in-plane and out-of-plane displacements.As such, there are aspects of Stereo-DIC which are identical to those of 2D DIC; namely calibration and correlation.Stereo-DIC relies upon epipolar geometry to determine in-plane and out-of-plane displacements in the real world from inplane displacements determined within images.Section 2.5 discusses how calibration and correlation are employed by Stereo-DIC based on epipolar geometry and how triangulation methods are used to perform displacement transformations.
Although Stereo-DIC can employ more than two cameras [50], ADIC3D is limited to the conventional two-camera implementation.Variables that need to be defined on a per camera basis have the subscript j to indicate that they are associated with the j th camera (i.e., camera 1 and 2).Within this work, an image series is defined as the collection of images captured by a single camera and an image set refers to the two corresponding image series.Thus, an image pair refers to the two images, one from each image series, which were captured simultaneously.

Homogeneous Coordinates
Calibration and triangulation use homogeneous coordinates [57], to represent coordinates in projective space, such that affine and projective transformations can be implemented through matrix multiplication.Appending a scaling variable of unity to a n dimensional vector in Euclidean space converts it to an n + 1 dimensional vector in projective space.Converting back to Euclidean space, represented as function Φ, entails dividing the vector's elements by the scaling variable before eliminating it to reduce the vector to n dimensions.Underlined variables denote homogenous coordinate vectors.

Calibration
Calibration, explained in terms of a single camera since each camera is calibrated separately, determines the parameters of the camera model, which involves the four coordinate systems (CSs) depicted in Figure 1: (i) the three-dimensional (3D) world CS having arbitrary orientation and position; (ii) the 3D camera CS with its origin at the aperture of the camera and its z-axis aligned with the optical axis; (iii) the 2D ideal sensor CS which is coincident with the plane of the charge coupled device (CCD) representing the image in the absence of distortion; and (iv) the 2D distorted sensor CS which is coincident with the ideal sensor CS but represents the actual image by accounting for radial distortion.No symbols have been assigned to the camera CS since it is not directly involved in DIC.
The camera model uses the pinhole camera model to project a 3D coordinate of a point on the specimen in the world CS to its corresponding 2D location in the ideal sensor CS before the radial distortion model determines its location in the distorted sensor CS.Although this camera model, and subsequent calibration process, is rather simple, it performs sufficiently well as shown in Section 4.However, the accuracy the displacements determined could be improved by using a more sophisticated camera calibration method which provides a more complete description of the camera optics such as the generic camera calibration method [58,59].

Pinhole Camera Model
The projection matrix Q j , consisting of intrinsic (K j ) and extrinsic (V j ) parameter matrices, relates a coordinate in the world CS, xw = xw ŷw ẑw 1 T , to a coordinate in the ideal sensor CS, xs j = xs j ŷs j 1 T , as Here α j is the arbitrary scaling variable of the homogenous coordinates.First, the rotation matrix, R j , and translation vector, T j , transform xw to the camera CS.Thereafter this 3D coordinate is projected to the ideal sensor CS using K j which performs four operations: (i) projects the 3D coordinate to a 2D coordinate on the plane of the CCD; (ii) scales the coordinate from metric units to units of pixels using ratios ξ x j and ξ y j in the x-and ydirections respectively; (iii) applies a translation, c x j and c y j in the x-and y-directions respectively, placing the origin of the ideal sensor CS at the top left of the image; and (iv) converts from the orthogonal camera CS to a skew ideal sensor CS using c s j .In this work c s j = 0 since an orthogonal ideal sensor CS is assumed.K j is only dependent on the camera system and does not change if the camera position changes.In contrast, V j depends on the relative position and orientation between the camera and the world CS.

Radial Distortion Model
The radial distortion model is sufficient to account for distortion in the images caused by imperfections of the lens system [60].First the normalized, ideal image coordinates xn j = xn j ŷn j T are computed as Thereafter, the unit-less radial distortion parameters, κ 1 j and κ 2 j , are used to compute the normalized, distorted image coordinates, x n j = x n j y n j T , as [61] x Finally, the distorted coordinates in the distorted sensor CS, x s j = x s j y s j T , are obtained as

Determining Calibration Parameters
Prior to performing calibration, a calibration image series must be captured of the calibration plate.The calibration plate is a planar object containing distinct, point-like features, called calibration targets (CTs), at known 3D locations which define the position and orientation of the world CS.A corresponding set of 2D coordinates in the distorted sensor CS are obtained by locating the CTs in the calibration images.These sets of 3D and 2D coordinates are used to solve for the parameters of the camera model, which describe the relationship between the two.Using at least 50 calibration images is recommended since this inverse problem is sensitive to noise in the 3D and 2D coordinates [41,62].
Calibration is performed using two steps.First, initial estimates for K j and V j are obtained using the method proposed by Zhang [61] while κ 1 j and κ 2 j are set to zero.Secondly, non-linear, least-squares optimization is used to optimize all the parameters by minimizing the total projection error, E proj , given as ( CT in the distorted sensor CS determined from the calibration image series and predicted by the camera model respectively.There is a total of L many CTs visible within the calibration image series.

Correlation
ADIC3D uses correlation to track a subset, a cluster of pixels, between images.The specimen's surface requires an isotropic, random pattern with high information content such that subsets are tracked reliably.Correlation operates on two images: a reference image, F, representing the specimen in its relaxed state and a deformed image, G, representing the specimen after experiencing displacement and/or deformation.Correlation determines how a subset of the reference image, called a reference subset (f), must displace and deform such that it matches a corresponding subset in the deformed image, called an investigated subset (g).This is referred to as the correspondence problem, which is mathematically defined as the objective function.For a more detailed discussion on the correspondence problem, the reader is referred to the authors' previous paper [53].
Note that correlation operates solely within the distorted sensor CS and is unconcerned with which camera the distorted sensor CS belongs to.Thus, in Sections 2.3 and 3.2, the subscript s j is dropped from x s j (i.e., x) to simplify the mathematics.
To understand how the objective function is comprised of the SF, the interpolation method and correlation criterion first consider that the i th pixel position of f, x i = x i y i T , is represented by the reference subset center position, x o = x o y o T , and the distance from x o to x i , ∆x i = ∆x i ∆y i T , as The SF, W, determines how the pixel positions of f displace as a result of the displacement and deformation of f, which is quantified by the shape function parameters (SFPs), P. The SF does this by modifying ∆x i as The pixel positions of f after displacement and deformation, x i = x i y i T , (referred to as query points) are determined as Thus ∆x i represents the distance from the center of f to x i .The light intensity values of G are interpolated at x i to obtain g.Thus, g represents a portion of G at a location described by the displacement SFPs with the deformation defined by the SFPs having been removed.The correlation criterion quantifies the degree of similarity between f and g.

Shape Function
The main SFs are the zero (W SF0 ), first (W SF1 ), and second-order SFs (W SF2 ) given as [43] Their sets of SFPs are given as Here u and v represent the displacement of f in the x-and y-directions respectively.Furthermore, u x , u xx , v y , and v yy define stretching of f, whereas u y , v x , u yy , v xx , u xy , and v xy define shearing.Although higher-order SFs more reliably determine displacements in complex displacement fields by accounting for more complex deformation, the DIC process is only concerned with the displacements determined by correlation.

Interpolation
ADIC3D employs the bi-cubic b-spline interpolation method detailed in the work of Hou et al. [63].Letting F and G represent interpolation functions, the light intensity values of pixel i of the reference, f i , and investigated subsets, g i , are interpolated as Images are filtered using a Gaussian low-pass filter, detailed by Pan [45], since otherwise the sensitivity of bi-cubic b-spline interpolation to high frequency noise would lead to bias in the displacements.The Gaussian filtering parameters, being the window size (β) and the standard deviation of the Gaussian function (σ g ), should be chosen to reduce bias without increasing variance.

Correlation Criterion
The objective function employs the more computationally efficient ZNSSD correlation criterion (C ZNSSD ), having a range of {C ZNSSD ∈ R|0 ≤ C ZNSSD ≤ 4} where smaller values indicate higher similarity, given as where Here there are I many pixels per subset.The correlation coefficient is reported as the zero-mean normalized cross correlation criterion (ZNCC), C ZNCC , since its range, {C ZNCC ∈ R|−1 ≤ C ZNCC ≤ 1} where larger values indicate higher similarity, is more intuitive.These correlation criteria are related as [64]

Objective Function
The objective function, based on Equation (12), is given as Therefore, correlation seeks the SFPs that define g to contain a light intensity pattern similar to f thereby optimizing C ZNSSD .However, the IC-GN optimization method requires modifying the objective function in order to derive an iterative optimization equation from it.More specifically, the current estimate of the SFPs, P, is applied to the investigated subset while the iterative improvement of the SFPs, ∆P, is applied to the reference subset as

Optimization Equation
The optimization equation derived from Equation (17), derivation detailed in [53], is given as Here H is the Hessian, defined in Equation (19), and the remaining terms (within the summation) form the Jacobian, J. H, being independent of P, is precomputed prior to beginning iterations.
is the light intensity gradient of f at pixel i. ∂W i ∂P is the Jacobian of the SF for pixel i and based on the SF order is given as and Although, each iteration of Equation (18) attempts to determine how f should displace and deform (according to ∆P) such that it matches g, f is not displaced or deformed in practice.Instead ∆P is used to determine the updated SFPs, P update , which serves as the current estimate in the next iteration.The stopping criterion deems when further iterations are redundant.

Updating the SFPs
P update is obtained by composing the inverse of ∆P with P as Here ω, being dependent on the SF order, uses the SFPs to populate a square matrix as [43] where

Stopping Criterion
The stopping criterion, ∆P , determines the magnitude of the change in the SFPs between iterations, based on ∆P, as [49] ∆P SF0 = ∆u 2 + ∆v 2 0.5 , and Here ζ is the largest value of ∆x i or ∆y i for the subset.Iterations cease once ∆P falls below a predefined value called the stopping criterion value.

Epipolar Geometry
Epipolar geometry is the projective geometry between the ideal sensor CSs, i.e., between camera 1 and 2. Consider the following definitions aided by Figure 2: (i) the baseline is formed between the two origins of the camera CSs; (ii); the epipolar plane is formed between the 3D coordinate and the baseline; (iii) the epipolar lines, λ j , are the intersection of the epipolar plane and the x-y plane of the ideal sensor CSs; and (iv) the epipoles, given as e s j = e 1 s j e 2 s j e 3 s j T in homogeneous coordinates, are located at the intersection of the baseline and the x-y plane of the ideal sensor CSs.For 2D coordinates in the ideal sensor CSs of the first, xs 1 , and second cameras, xs 2 , to correspond to the same xw requires that these coordinates satisfy the epipolar constraint defined as xT More specifically, since the fundamental matrix, B, maps a point in one ideal sensor CS to its corresponding epipolar line in the other ideal sensor CS, the 2D coordinates must lie on the epipolar lines of their respective ideal sensor CSs.B encapsulates the intrinsic and extrinsic relation between the ideal sensor CSs.
Computing B requires applying a homography, Ω, to the projection matrices of the cameras to transform them to canonical form (Q c j ), where the world CS coincides with the camera CS of the first camera, as [65] where Here I is a 3 × 3 identity matrix.Q c 2 , which consists of a 3 × 3 matrix M and a 3 × 1 vector m, transforms a coordinate from the camera CS of the first camera to its location in the ideal sensor CS of the second camera.B is given as where [m] × is the skew-symmetric matrix of m = m 1 m 2 m 3 T defined as

Stereo-DIC Overview
The core operations of Stereo-DIC are subset matching, calibration, and displacement transformation.For a point on the specimen, subset matching determines the two sets of in-plane, pixel displacements it experiences in each distorted sensor CS which, along with the calibration parameters determined by the calibration process, are used by displacement transformation to determine the in-plane and out-of-plane (thus 3D) metric displacements it experiences in the world CS.Note that with regards to subset matching and displacement transformation, the term in-plane refers to the x-y plane of the distorted sensor CSs and world CS respectively.
Subset matching is composed of stereo and temporal matching which utilize correlation differently.First, subset matching breaks up the first image of the first image series (FIS1) into subsets.Second, stereo matching determines the location of the corresponding subsets in the first image of the second image series (FIS2).Such subsets, one from each image of an image pair which correspond to the same point on the specimen, are termed a subset pair.Lastly, temporal matching, processing each image series separately, determines the displacements experienced by these subsets through their respective image series.This is achieved by sequentially stepping through an image series, starting at the second image, and at each step performing correlation by treating the current image as the deformed image while the reference image is dictated by the reference image strategy as discussed in Section 2.6.1.Thus, stereo matching is unique to Stereo-DIC while temporal matching is identical to that of 2D DIC.
The calibration method is identical to that of 2D DIC.However, in order to determine the epipolar geometry from the calibration parameters requires that the two calibration image series, used to calibrate the cameras separately, form part of the same calibration image set.
Displacement transformation, illustrated in Figure 3, steps through image pairs of the image set similarly to temporal matching.More specifically, considering a single reference subset pair with known locations, x o s j , it uses the in-plane displacements determined by temporal matching (u j and v j in the x-and y-directions) between the reference images (of the reference image pair) and the deformed images (of the deformed image pair) to determine the locations of the investigated subset pair, x d s j .In contrast to 2D DIC, distortion correction alone is insufficient to determine the ideal coordinate of a subset in the ideal sensor CS.More specifically, discretization of light to form an image introduces noise in the images [66] and displacements measured by correlation, which are used to determine the locations of subsets, inherently contain some noise.Thus, correcting x o s j and x d s j for radial distortion results in the measured coordinates of the reference, xo s j , and investigated subset pair, xd s j , in the ideal sensor CS respectively.Generally, these measured coordinates of a subset pair, represented as xs j , do not satisfy the epipolar constraint of Equation ( 24) since their back-projection rays do not intersect.Thus, triangulation, which determines the 3D coordinate in the world CS, xw , from xs j , is performed in two steps.First the polynomial triangulation method, detailed in Section 2.7, determines the ideal coordinates of a subset pair, xs j , which lie as close to xs j as possible while satisfying the epipolar constraint.Thereafter the linear triangulation method, detailed in Section 2.8, determines xw to which xs j correspond as the intersection of their back-projected rays.Performing triangulation on xo s j and xd s j results in xo w and xd w respectively.Finally, xo w is subtracted from xd w to determine the displacement experienced by the point on the specimen, in the x-, y-, and z-directions represented as ûw , vw , and ŵw respectively, between the time stamps at which the reference and deformed image pairs were captured.Note that subscript w indicates coordinates or displacements in the world CS, whereas variable w indicates displacement in the z-direction.This is done for all subset pairs across all image pair combinations, analyzed by temporal matching, to determine the full-field metric displacements experienced by the specimen throughout the image set.

Subset Matching
The optical flow determined by temporal and stereo matching is a result of the displacement and/or deformation experienced by the specimen (over time) and the perspective change between the cameras respectively.They differ predominantly in the method used to obtain SFP initial estimates.

Temporal Matching
There are two reference image strategies for temporal matching: absolute and incremental.The absolute and incremental strategies define the reference image as either the first image of the image series or the previous image, relative to the current deformed image, respectively.The incremental strategy more reliably tracks more severe deformations between images but incurs accumulative errors when determining the displacement relative to the first image.
Temporal matching uses the Phase Correlation Method (PCM) to determine initial estimates for the displacements between two subsets.PCM efficiently computes the correlation coefficients in the frequency domain for a range of integer displacements of up to half the subset size in each direction between the two subsets.The integer displacements with the best correlation coefficient are identified and used as initial estimates for the displacement SFPs while the deformation SFPs are set to zero.Consult the work of Foroosh et al. [67] for a discussion on PCM.For the absolute strategy, PCM is only used for the first correlation run.Thereafter the SFPs of the previous correlation run, of the same image series, are used as an initial estimate for the SFPs of the current correlation run.The incremental strategy uses PCM for each correlation run.

Stereo Matching
Determining reliable out-of-plane displacements in the world CS requires an angle between the optical axes of the cameras of between 15 • and 35 • [5,68].As such, stereo matching uses the second-order SF since subsets experience complex deformation between the FIS1 and FIS2 due to this perspective change [41].
Traditionally, template matching, which utilizes image rectification, was used to determine SFP initial estimates for stereo matching [69].However, image rectification is susceptible to errors in the calibration parameters, can cause aliasing in the rectified images and interpolation during image rectification can lead to errors in the determined SFPs [70].
As such, feature matching methods are becoming favored [41].ADIC3D's feature matching method is based on that proposed by Zhou et al. [71].It uses the scale-invariant feature transform (SIFT) feature matching algorithm [72] to identify corresponding keypoints between the images (where a keypoint is a unique feature in an image) and the affine transformation mapping model to mathematically relate the coordinates of these keypoint pairs between the images.
The SIFT algorithm is used because it is robust in translation, rotation, image scaling, variation in illumination, affine transformation and moderate perspective change.SIFT identifies keypoints that are scale invariant and computes a descriptor, in the form of a 128-dimensional vector, for each keypoint based on the light intensity gradient information in the vicinity of the keypoint.
For each keypoint of the FIS1, the Euclidean distance between its descriptor and that of every keypoint of the FIS2 is determined.The two keypoints of the FIS2 with the smallest descriptor distances are identified.If the ratio of the smallest to second smallest descriptor distance is less than 0.8, the keypoint of the FIS2 with the smallest descriptor distance is designated as the matching keypoint of the keypoint of the FIS1 [72].Otherwise, this keypoint of the FIS1 is discarded.
The affine transformation relates a keypoint pair as Here T are the locations of the k th keypoint in the FIS1 and FIS2 respectively while a 1 through a 6 are the affine transformation parameters.Given three or more keypoint pairs, the affine transformation parameters can be solved either as a system of linear equations or through linear least-squares respectively.The K many x k s 1 that are nearest to the subset's center position are used to determine the affine transformation parameters of the subset.However, since SIFT often returns false keypoint pairs, the m-estimator sample consensus (MSAC) method [73] is used to remove these false keypoint pairs when determining the affine transformation parameters.False keypoint pairs are identified as those which have a squared error distance greater than a predefined threshold, τ.The squared error distance of the k th keypoint pair, E dist k , is calculated as ADIC3D uses K = 20, τ = 1 pixel 2 and a confidence, specifying the probability that MSAC finds the maximum number of inliers, of 99.5%.Second-order SFP initial estimates are determined from the affine transformation parameters and the subset position in the FIS1, Once correlation has optimized these SFPs, the corresponding subset position in the FIS2, x s 2 = x s 2 y s 2 T , is determined as

Polynomial Triangulation Method
The polynomial triangulation method, proposed by Hartley and Sturm [66], determines xs j by minimizing the geometric error cost function, D, given as Here function δ determines the Euclidean distance between two points.Determining xs j by minimizing the distance between the coordinates xs j and xs j within the ideal sensor CSs ensures that the polynomial triangulation method is projective-invariant. Projectiveinvariance implies that the projective frame within which xw is defined does not affect the values of xs j [65].The polynomial method assumes that the calibration parameters are known with greater accuracy than the measured coordinates [66].Furthermore, it is assumed that neither xs j nor xs j coincides with an epipole since this would lead to an unfavorable solution for the 3D coordinate.
It is known that a pair of ideal coordinates which satisfy the epipolar constraint lie on a pair of corresponding epipolar lines.Additionally, the epipolar plane can only rotate about the baseline since the camera CSs, and thus epipoles, are fixed.Thus, there is a family of epipolar lines upon which the idea coordinates can lie.Furthermore, the projective-invariance of the polynomial triangulation method allows applying a separate rigid transformation to each ideal sensor CS in order to parameterize the family of epipolar lines in terms of a single variable t.As such, homogeneous coordinates are utilized and Equation (32) becomes where λ j (t) is a vector representing an epipolar line in the ideal sensor CS as a function of t while δ xs j , λ j (t) represents the distance from xs j to its orthogonal projection on λ j (t).
The rigid transformations first apply a translation matrix, T pm j , to each ideal sensor CS translating xs j to the origin of its CS.
The original fundamental matrix, B, determined using Equation ( 27) needs to be updated as The epipoles e s 1 and e s 2 of the first and second cameras are computed, using singular value decomposition, as the right and left null-space of B 1 respectively.Thereafter they are normalized such that e 1 The fundamental matrix is updated as Since B 2 1 0 e 3 s 1 T = 1 0 e 3 s 2 B 2 = 0 0 0 T the elements of the fundamental matrix are represented as An epipolar line of the first camera, λ 1 (t), passing through point 0 t 1 T on the y-axis of the ideal sensor CS is represented as The corresponding epipolar line of the second camera, λ 2 (t), is computed as 2 (ϕ 3 t+ϕ 4 ) 2 and so Equation (34) has the form Finding the optimal value of t involves taking the derivative Equation (42) in terms of t, collecting terms over a common denominator and equating the resulting numerator to 0 as The 6 roots of Equation ( 43) are solved for and Equation ( 42) is evaluated at the real parts of these roots to determine the value of t which corresponds to a minimum.The asymptotic value of Equation (42) should also be evaluated as t → ∞ which corresponds to an epipolar line in the first image of e 3 s 1 xs 1 = 1 [65].The optimal t is used to determine xs 1 and xs 2 as and Note that the inverse of the rigid transformations are applied during Equations ( 44) and (45) such that xs 1 and xs 2 are in the original ideal sensor CSs.

Linear Triangulation Method
Denoting the n th row of Q c j as q n c j the projection of xw to xs 1 is described by the following three equations Similarly, for xs 2 this is given as Using the third relation of Equations ( 46) and ( 47) to remove the scaling variables, α 1 and α 2 , the following relation is obtained.
xw is computed, using singular value decomposition, as the unit eigenvector of S T S associated with the smallest eigenvalue [65].

Displacement Transformation
Consider a single subset pair and its two sets of in-plane displacements determined by temporal matching between the reference and deformed image pairs under consideration.
Thereafter x o s j and x d s j are undistorted to obtain xo s j = xo
ADIC3D returns a structured array ResultData containing four main fields: structured variable Stereo holding stereo-matching data as detailed in Table 2; structured arrays ProcData1 and ProcData2, detailed in Table 3, holding the processing data for temporal matching of subsets of the first and second image series respectively; and structured array DispTrans containing displacement transformation data detailed in Table 4.

FileNames1
Cell array of character vectors containing the image file names of the first image series.All images need to be the same size.

FileNames2
Cell array of character vectors containing the image file names of the second image series.All images need to be the same size.

Mask
Logical matrix, which is the same size as the images, indicating which pixels should not be analyzed during correlation.

Variable Variable Description
P(b,q) SFPs (b = 1 for u and b = 7 for v).

Iter(q)
Number of iterations until stopping criterion is satisfied (maximum of 100 iterations).

StopVal(q)
Final stopping criterion value for subset q.

Variable Variable Description
ImgName Deformed image name.

ImgFilt(b)
Standard deviation (b = 1) and window size (b = 2) for the Gaussian filter respectively in pixels.

Xos(b,q)
Reference subset position in the distorted sensor CS of the relevant camera (b = 1 for x o and b = 2 for y o ).

P(b,q)
SFPs (b = 1 for u and b = 7 for v).

Iter(q)
Number of iterations until stopping criterion is satisfied (maximum of 100 iterations).

Table 4.
Accessing the displacement transformation data for image d (contained in ResultData.DispTrans(d)) for subset q.

Variable Variable Description
Xow(b,q) Reference subset position in the world CS (b = 1 for xo w , b = 2 for ŷo w , and b = 3 for ẑo w ).
Refer to Appendix B for brief discussion of utilizing parallel processing to reduce computation time of ADIC3D.

ADIC3D Function
The main function ADIC3D, outlined in Table 5, sets up the DIC problem and calls the appropriate subroutines.The output variables are preassigned in lines 9-12 to save input data and efficiently store computed data.ADIC3D allows subsets to be analyzed in a flexible manner by assigning Xos, SubSize, SubShape, and SFOrder on a per subset and per image basis.However, despite being capable of this, ADIC3D assigns the same SubSize, SubShape, and SFOrder to each subset as is conventional for DIC.Subsets containing invalid pixels, identified by Mask, are eliminated.The subroutine StereoMatch is called to perform stereo matching.StereoMatch has input variables n (the total number of images in each image series), ResultData, FileNames1, FileNames2, and StopCritVal.It returns the subset positions in the FIS2 stored as ProcData2(d).Xos and stereo matching correlation results P, C, Iter, and StopVal stored in Stereo.
Temporal matching is performed on each image series separately in lines 15 and 17 using the ImgCorr subroutine.ImgCorr's inputs are n, ProcData1 or ProcData2, FileNames1 or FileNames2, RefStrat and StopCritVal.It returns P, C, Iter, and StopVal for each subset throughout the image series stored in ProcData1 or ProcData2 for the first and second image series respectively.
Finally, subroutine CSTrans determines coordinates and displacements in the world CS based on inputs n, ResultData, WorldCTs, ImgCTs, and RefStrat.CSTrans returns Xow, Uw, and CamParams stored in DispTrans.Note that within the subroutines ResultData is shortened to RD and ProcData is shortened to PD where appropriate.

Correlation Implementation
Both stereo and temporal matching make use of correlation, detailed in Section 2.3, which is implemented using the three subroutines: SubShapeExtract, SFExpressions, and SubCorr.These subroutines are discussed briefly since they are identical to those used in ADIC2D [53].
SubShapeExtract, detailed in Table 6, determines f i , ∇f i and ∆x i of a subset based on inputs SubSize, SubShape, Xos, F, ∇F, and SubExtract.Note that variables with subscript i refer to the full set of this variable for the subset (i.e., f i refers to f i ∀i ∈ I) throughout Section 3.2.∇F is the light intensity gradient of the reference image whereas SubExtract is an anonymous function which extracts a square subset from a matrix based on the size and position of the subset.For SubShape defined as 'Circle', SubSize specifies the diameter.Compute ∆x i using SubSize; Line 13 Determine mask of elements that fall within the circular subset; Line 14-16 Use mask to extract appropriate data for circular subset; Line 17 end switch SFExpressions returns anonymous functions for expressions that are dependent on the SF order, as outlined in Table 7.These anonymous functions are as follows: W defines Equation (9); dFdWdP defines ∇f i ∂W i ∂P ; SFPVec2Mat defines Equation ( 22); Mat2SFPVec extracts the SFPs from SFPVec2Mat; and StopCrit defines Equation (23).ADIC3D allows for the zero, first, or second-order SFs to be used by setting SFOrder to 0, 1, or 2 respectively.SubCorr performs the core operations of correlation for a subset as outlined in Table 8.Its inputs are the interpolation coefficients of the deformed image, f i , ∇f i , SubSize, SFOrder, Xos, ∆x i , initial estimates for P and StopCritVal.SubCorr returns P, C, Iter, and StopVal.The size of vector P is consistent with that of the second-order SF.Thus, the elements of P SF2 , not relevant to the specified SF order, remain zero.SubShapeExtract is called prior to calling SubCorr in order to determine data used by SubCorr while SFExpressions is called within SubCorr to define the appropriate anonymous functions on a subset basis.Thus, the effect of the subset shape and SF order as well as the core mathematical operations of correlation can be considered separately.

Stereo Matching Implementation
Excluding the subroutines used to perform correlation, stereo matching is implemented using two additional subroutines: FeatureMatch and StereoMatch.FeatureMatch returns SFP initial estimates, determined according to the feature matching method outlined in Section 2.6.2, based on inputs ProcData1, d (the image to be analyzed), F, G, and SubExtract.StereoMatch refines these SFPs through correlation and uses them to determine the subset positions in the FIS2.

StereoMatch Function
StereoMatch, detailed in Table 9, loads the FIS1 and FIS2 as F and G, respectively, and performs Gaussian filtering on these using MATLAB's imgaussfilt function.Lines 4 and 6 can be changed to use alternative filtering methods.Bi-cubic b-spline interpolation coefficients are determined for the deformed image using MATLAB's griddedInterpolant function.In line 7 'spline' can be replaced with 'linear' or 'cubic' for bi-linear or bi-cubic polynomial interpolation respectively.Moreover, line 7 can be replaced with an alternative interpolation method such as MATLAB's spapi function.Line 9 of SubCorr may need to be updated if alternative interpolation methods are employed.Lines 11 to 16 cycle through each subset of the FIS1, for which FeatureMatch in line 9 determined SFP initial estimates, performing correlation using subroutines SubShapeExtract and SubCorr.These displacement SFPs are used to determine the reference subset positions in the second image series.
Subset pairs which achieve a ZNCC coefficient below 0.6, or which fall outside the bounds of the FIS2, are assumed to have failed stereo matching.These subsets are determined in lines 21-22 and their subset positions in both image series are set to NaN (not a number) in lines 23-26 such that they are not analyzed during temporal matching or displacement transformation.
StereoMatch returns ResultData with the subset positions in the FIS2 assigned to ProcData2(d).Xos; the correlation results for stereo matching (P, C, Iter, StopVal) are assigned to Stereo.

FeatureMatch Function
FeatureMatch, outlined in Table 10, implements SIFT feature matching using codes available from the VLFeat open-source library of computer vision algorithms.More specifically, VLFeat's vl_sift function is used to determine the keypoint locations and descriptors of F and G in lines 3-4.The keypoints of F which fall within the perimeter of at least one subset (square subsets equivalent in size to that specified for the subset under consideration are used) are identified in line 5.Only these keypoints of F are retained, in line 6, in order to reduce the computation time of VLFeat's vl_ubcmatch function in determining matching keypoint pairs between F and G in lines 7-8.

Line Numbers Task Performed
Line 2 if VLFeat library is not setup, then return an error; Line 3-4 Compute keypoint locations and descriptors for F and G using vl_sift; Line 5 Determine vector KptsInVacinity identifying the keypoints of F which fall within the perimeter of the square subsets, equivalent in size to that specified for the subset under consideration, of F; Line 6 Eliminate keypoints (and their associated descriptors) of F which do not fall within the perimeter of any of the subsets; Line 7-8 Determine matching keypoints using vl_ubcmatch; Line 9 Determine the 20 nearest keypoints for each subset stored in matrix relevantKpts; Line 10-11 Define anonymous functions RansacModel and RansacError for determining affine transformation parameters (of Equation ( 29)) and evaluating Equation (30) respectively; Line 12 Initialise P as NaNs; Line 13 for subset number q = 1 to number of subsets, do Line 14 try Line 15 Use MSAC to determine affine transformation parameters for subset q from its relevant keypoints; Line 16 Convert affine transformation parameters to second-order SFPs using Equation (31); Line 17 end try Line 18 end for Line 19 Display information for SIFT feature matching; For each subset, the 20 closest keypoints to its center position are identified in line 9.The for loop of line 13 cycles through the subsets using MSAC to determine reliable affine transformation parameters, which are used by Equation (31) to compute the corresponding second-order SFPs.
MSAC is implemented using MATLAB's ransac function.Two anonymous functions are used in the call to ransac.RansacModel, defined in line 10, determines the affine transformation parameters from the keypoint locations by solving Equation (29).RansacError, defined in line 11, evaluates Equation (30) to determine each keypoint's squared error distance (for keypoints identified to be relevant to the subset).The fifth argument of ransac sets τ = 1 pixel 2 and the final argument specifies the confidence that the maximum number of inliers are determined.The fifth and sixth arguments of ransac, as well as the number of nearest keypoints used for each subset (line 9), can be altered to suit the reader's needs.The try statement of line 14 catches instances where ransac fails to identify three keypoint inliers.As such, the SFP initial estimates of these subsets, stored in P, remain NaN as initialized on line 12.

Temporal Matching Implementation
Excluding the subroutines used to perform correlation, temporal matching is implemented using PCM and ImgCorr.ADIC3D's implementation of temporal matching is discussed briefly since it is almost identical to that of ADIC2D [53] with the exception that PCM and ImgCorr are modified to avoid analyzing subsets which fail stereo matching.

PCM Function
PCM determines an initial estimate of the displacement SFPs of a subset based on inputs F, G, SubSize, the subsets position, and SubExtract as outlined in Table 11.Subset positions are passed as separate x and y positions as required by MATLAB's arrayfun function which is used to call PCM.The if statement of line 2 avoids analyzing subsets which contain NaN positions and returns NaN displacements in this case.12, makes use of a nested for loop to perform temporal matching.The outer loop cycles through the images of an image series, starting at the second image, preparing data for correlation.This includes loading the appropriate reference and deformed images; filtering the images; determining the interpolation coefficients of G; determining the reference image gradient; shifting the reference subset positions by the displacement SFPs determined in the previous correlation run (these displacement SFPs are rounded down such that the reference subset need not be interpolated [75]); and using PCM to determine initial estimates of the SFPs for the current correlation run.In the case of the absolute reference image strategy, the reference subset positions are not shifted and PCM is only used to obtain an initial estimate of the displacement SFPs for the first correlation run.For subsequent iterations, the SFPs of the previous correlation run are used as initial estimates.
The inner for loop cycles through the subsets calling SubShapeExtract and SubCorr to perform correlation of the subsets.The if statement of line 19 ensures that only subsets which pass stereo matching are analyzed within the for loop.

Displacement Transformation Implementation
Displacement transformation is implemented using subroutines CSTrans and Triangulation.

CSTrans Function
CSTrans, outlined in Table 13, uses MATLAB's estimateCameraParameters function to determine the calibration parameters according to Section 2.2 and the fundamental matrix according to Section 2.4.The canonical form of the projection matrices are used such that the world CS is aligned with the camera CS of the first camera.
Lines 6-18 cycle through image pairs of the image set performing displacement transformation according to Section 2.9.Line 9 identifies subsets which do not contain NaN values in either their position or SFPs (subsets which pass stereo matching) which are to be analyzed during displacement transformation.MATLAB's undistortPoints function is used to remove distortion from the subset positions while the subroutine Triangulation uses polynomial and linear triangulation to determine the 3D coordinate in the world CS corresponding to a subset pair.The if statement of line 10 avoids recomputing the reference subset positions for the absolute reference image strategy.
CSTrans is susceptible to particular issues which can cause a fatal error.Refer to Appendix C for a discussion of a function which improves the robustness of CSTrans.

Triangulation Function
Triangulation, outlined in Table 14, cycles through the subset pairs determining their ideal coordinates in the ideal sensor CSs using the polynomial triangulation method according to Section 2.7, in lines 3-20 (based on the code of Lourakis [76]), before performing linear triangulation according to Section 2.8, in lines 21-22, to determine their 3D coordinate in the world CS.DICe was run using the same parameters as ADIC3D.StrainMaster, a streamlined version of LaVision's DaVis software, was used to analyze Samples 1 and 2. Its streamlined nature prevents modifying some correlation parameters such that it employed circular subsets, a stopping criterion of 0.01, and a maximum of 200 iterations.Although this prohibits its results from being directly comparable to ADIC3D, it is included to reflect the capabilities of a commercial code.All codes used identical masks for each sample.For ADIC3D, subsets achieving a ZNCC criterion below 0.8 for stereo matching were eliminated while DICe and StrainMaster have their own methods of eliminating subsets which fail stereo matching.The final sets of valid subset pairs for ADIC3D and DICe were determined in two steps such that the results of the two codes are directly comparable.
First, a preliminary set of valid subset pairs was determined for each code as the subsets which achieved a ZNCC criterion above 0.8 for temporal matching.Additionally, for Samples 1 and 2 of DICe, Grubbs's test [79] was used to remove outliers, in terms of displacement error, to reduce outliers occurring in its preliminary set.Secondly, the final set of valid subset pairs for both codes was determined as the intersection of these preliminary sets based on the subset positions in the FIS1.
Successive images of all samples exhibited a jump in displacement such that PCM failed to determine SFP initial estimates.As such, subroutine FeatureMatch was used to obtain SFP initial estimates during temporal matching by inserting "[PD(d).P]=FeatureMatch (PD,d,F,G, SubExtract);" after line 16 in subroutine ImgCorr.DICe also used feature matching to determine SFP initial estimates.
The provided calibration image sets, of dot-grid-style calibration plates, were used by each code to perform calibration.ADIC3D used the STEP1_CalcDLTparameters function of MultiDIC, an open-source Stereo-DIC software proposed by Solav et al. [50], to locate the CTs in the distorted sensor CSs.All experimental image sets were captured using the FLIR Grasshopper 2 model Gras-50S5M cameras with (Edmund Optics) lenses having a focal length of 35 mm.

Samples 1 and 2
Sample 1 captures images of a complex specimen geometry, consisting of a flat base plate with protruding 3D features, as it is displaced in increments of 10 mm in the xand z-directions according to Table A2 in Appendix D. The specimen is displaced by an Aerotech nano-positioning stage (model ANT130-160-XY) with a high-precision linear encoder and position feedback control such that the true displacements are known with high certainty.Sample 2 simulates the experiment of Sample 1 using the synthetic image generator documented by Balcaen et al. [62] to minimize error sources in the image set.
The world CS of each DIC code differs from that within which the displacements were imposed.Thus, for each code Pythagoras's theorem is used to compute the true, μtrue w q , and calculated, μcalc w q , displacement magnitudes of each subset q.The resulting errors in displacement magnitude are quantified as bias, computed as the mean of the absolute error (MAE), and variance, computed as the standard deviation of the absolute error (σ).
Here Q is the number of valid subset pairs in the image pair analyzed.For Samples 1 and 2, ADIC3D and DICe had 71,715 and 66,850 valid subsets, respectively, while StrainMaster had 112,021 and 108,085 valid subsets respectively.Tables 15 and 16 report errors metrics for Samples 1 and 2, respectively, for steps 3, 7, and 15 which correspond to displacement extremes.The improved error metrics of Sample 2 relative to Sample 1 for each code, with the exception of StrainMaster's variance for step 15, reflects the impact of the reduced error sources of the synthetic image set.For Sample 1, ADIC3D performs on par with DICe (up to 27% higher bias) and StrainMaster (at least 12% lower bias).In contrast, for Sample 2 ADIC3D shows improved error metrics relative to DICe and StrainMaster (at least 69% and 54% lower bias respectively).StrainMaster's higher bias and variance are a result of it analyzing more subsets in the vicinity of protruding features, which are challenging to track, which are removed for ADIC3D and DICe.

Sample 5
Sample 5 captures images of an ASTM tensile specimen made of Aluminium 2024-T6 as it is loaded in tension.The displacements computed by ADIC3D and DICe are directly comparable since both align their world CS with the camera CS of the first camera.
Figure 4 shows the ûw displacement field of image pair 1000 superimposed on the FIS1.The mean and maximum of the percentage difference of the ûw , vw , and ŵw displacements, between ADIC3D and DICe, are reported in Table 17 for selected image pairs.The cameras being positioned apart horizontally in the x-direction is reasoned to cause the larger percentage difference of vw and ŵw relative to ûw which is consistent with the findings of Balcaen et al. [5].The low percentage difference (remaining below 2% mean and 4% maximum) indicates that ADIC3D performs consistently with DICe.

Discussion
ADIC3D's modularity is threefold.Firstly, the main tasks are performed by separate subroutines, closely linked to the mathematical theory presented, allowing readers to progressively build up their understanding of ADIC3D.For instance, correlation is separated into three subroutines enabling the reader to consider: (i) how data is prepared based on the subset shape (SubShapeExtract); (ii) the SF order's influence on expressions for variables used during correlation (SFExpressions); and (iii) the core operations of correlation (SubCorr) separately.Furthermore, subroutines ImgCorr and StereoMatch show how correlation is utilized differently by temporal and stereo matching.
Secondly, the subset shape, SF order, interpolation method, and Gaussian filtering parameters can be readily changed.Although the influence of these on displacement results is well documented in experimental solid mechanics applications [3,4,45], this allows investigating their influence in novel applications to fine-tune ADIC3D's setup.
Thirdly, the SF order, subset size, and subset shape can be changed on a per subset and per image basis enabling straightforward integration with adaptive strategies.Adaptive strategies, which are gaining interest in the field of DIC, attempt to assign optimal parameters to each subset such that the resulting displacements are, theoretically, as reliable as possible and independent of the user.
In order for the code to retain its simplicity, four aspects, identified by Pan [41] as stateof-the-art requirements, were omitted.First, calibration does not use bundle adjustment, proposed by Furukawa and Ponce [80], to refine the calibration parameters.Additionally, Vo's method [81] of combing the frontal image plane and DIC to more accurately locate the CTs is not implemented.It should be noted that although the camera model presented in Section 2.2 only accounts for radial distortion, accounting for tangential distortion in ADIC3D is straightforward, as outlined in Appendix E.
Secondly, although the incremental reference image strategy is included, ADIC3D does not automatically employ it if the ZNCC coefficient falls below a certain threshold as recommended [82].Thirdly, ADIC3D uses bi-cubic b-spline interpolation whereas bi-quintic b-spline interpolation is recommended [3] since increasing the degree of the interpolation method reduces the errors in the "ultimate error regime" [83] especially for smaller subsets.
Lastly, ADIC3D does not employ Pan's [84] reliability guided displacement tracking strategy (RGDT).ADIC3D's temporal matching uses the SFPs of each subset of the previous correlation run as initial estimates in the current correlation run.In contrast, RGDT only does this for the subset achieving the best ZNCC in the previous correlation run.This subset's optimized SFPs are used as an initial estimate to correlate its neighboring subsets.Of the correlated subsets, having uncorrelated neighbors, the one with the best ZNCC has its SFPs used as an initial estimate to correlate its neighbors.This is repeated to correlate all subsets in the current correlation run.As noted with ADIC2D [53], this exclusion makes ADIC3D susceptible to propagating spurious SFPs of a subset throughout the image series.
Furthermore, RGDT is excluded such that ADIC3D is path-independent (where each subset is correlated independently of its neighbors).Thus, this avoids the need to specify seed points and the issue path-dependent methods have, with complex specimen geometries, of transferring SFPs across subsets which experience different deformation modes.This also allows ADIC3D to leverage MATLAB's parfor loop to utilize parallel processing to reduce computation time (as discussed in Appendix B).
Additionally, although not a state-of-the-art requirement, ADIC3D's feature matching method can be improved.Sample 1 and 2 highlight its limitation as it determines invalid SFP initial estimates for subsets in the vicinity of protruding features that experience a large perspective change.Iniyan et al. [85] recommend using the DeepFlow [86] feature matching algorithm combined with the homography transform and MSAC.
Despite these omissions, ADIC3D's performance is on par with established codes.More specifically, Samples 1 and 2 show that ADIC3D determines displacement magnitudes in the world CS with bias and variance consistent with that of DICe and StrainMaster, while Sample 5 shows that ADIC3D determines displacement vectors within a 4% maximum difference of those determined by DICe.Although Samples 1 and 2 reveal limitations of ADIC3D's stereo matching, it performs sufficiently well in the presence of a complex specimen geometry.This, coupled with ADIC2D's validation of temporal matching and thus the correlation method of ADIC3D indicates that ADIC3D performs sufficiently well for practical use in the field of experimental solid mechanics.Moreover, due to its modularity ADIC3D can be readily adapted to a wide range of applications across many fields.Furthermore, validation of this modular framework makes it both attractive as an education resource and as a basis to further the capabilities of DIC.

Conclusions
The theory of a subset-based, Stereo-DIC framework, that is predominantly consistent with current state-of-the-art techniques, and its implementation as a modular 202 line MATLAB code, was presented.This framework is an extension of the previously presented ADIC2D framework and as such is also comprised of square or circular subset shape selection, assigning of Gaussian filtering parameters, altering of the interpolation method, the zero, first and second-order SFs, reference image strategy selection, the Phase Correlation Method to determine SFP initial estimates for temporal matching and calibration through using MATLAB's image calibration toolbox.ADIC3D, being capable of performing Stereo-DIC to determine both in-plane and out-of-plane displacements in the world CS, is additionally comprised of stereo matching, SIFT feature matching to determine SFP initial estimates for stereo matching and displacement transformation using the polynomial and linear triangulation methods.Furthermore, ADIC3D allows for assignment of SF order, subset shape, and subset size on a per image and per subset basis enabling straightforward integration with adaptive strategies.The framework is modular, enabling the reader to gain a deeper understanding of DIC as well as encouraging readers to adapt ADIC3D for new applications.Validation of the full framework coupled with its modularity makes it appealing not only as an educational resource, but also as a starting point to develop the capabilities of DIC.Projection matrix of the j th camera Q c j Projection matrix of the j th camera in canonical form q n c j

Glossary of Symbols
The n th row of the canonical projection matrix of the j

Figure 1 .
Figure 1.Schematic diagram of coordinate systems employed by the camera model.

T
are the locations of the l th

Figure 2 .
Figure 2. Schematic illustration of the epipolar geometry.

Figure 3 .
Figure 3. Schematic illustration of the displacement transformation process.Note that a projection line represents the combined process of distortion correction and triangulation.

s 1 2 + e 2 s 1 2 = 1 and e 1 s 2 2 + e 2 s 2 2 = 1 . 1 T
The second step of the rigid transformation applies rotations R pm 1 and R pm 2 to the ideal sensor CSs of the first and second cameras respectively in order to place their epipoles on the x-axes of their respective CSs as R pm 1 e s 1 = 1 0 e 3 s and R pm 2 e s 2 = 1 0 e 3 s 2 T .
is performed using non-linear, least-squares optimization since an exact solution for the inverse ofEquation (3)  does not exist because it requires determining the roots of a polynomial of degree greater than four[74].Triangulation, represented by function Ψ, determines the 3D coordinates in the world CS to which the reference subset pair, displacements ûw , vw , and ŵw are determined as

Figure 4 .
Figure 4. Displacement in the x-direction in the world CS superimposed on the FIS1 of the ASTM tensile specimen[55].

TTTTk
th camera Q Number of subset pairs per image pair R j Rotation matrix for pinhole camera model of the j th camera R pm j rRotation matrix of the j th camera for the polynomial triangulation method S Matrix for the linear triangulation method σ g Standard deviation of Gaussian function σ Standard deviation of the absolute error T j Translation vector for pinhole camera model of the j th camera T pm jTranslation matrix of the j th camera for polynomial triangulation method t Variable for parameterization of epipolar lines τ Squared error distance threshold u Displacement in x-direction in the distorted sensor coordinate system u x , u xx , u y , u yy and u xy Derivatives of the x-direction displacement ûw Displacement in x-direction in the world coordinate system V j Extrinsic parameter matrix of the j th camera v Displacement in y-direction in the distorted sensor coordinate system v x , v xx , v y , v yy and v xy Derivatives of the y-direction displacement vw Displacement in y-direction in the world coordinate system W Shape function ∂W i ∂P Jacobian of the shape function, in terms of the shape function parameters, for pixel i ŵw Displacement in z-direction in the world coordinate system xw = xw ŷw ẑw T Ideal coordinate in the world coordinate system xs j = xs j ŷs j T Ideal coordinate in the ideal sensor coordinate system of the j th camera xn j = xn j ŷn j TNormalized ideal image coordinates of the j th camerax n j = x n j y n j T Normalized distorted image coordinates of the j th camera x s j = x s j y s j T Coordinate in the distorted sensor coordinate system of the j th camera Location of l th calibration target in distorted sensor coordinate system predicted by the camera modelx i = x i y i TPixel position of the i th pixel of the reference subset in distorted sensor coordinate systemx o = x o y o TCenter of reference subset in distorted sensor coordinate system∆x i = ∆x i ∆y i TDistance from the reference subset center to i th pixel position of reference subset∆x i = ∆x i ∆y i TDistance from the reference subset center to i th pixel position of investigated subsetx i = x i y i T i th pixelposition of the investigated subset in the distorted sensor coordinate system Investigated subset position in the distorted sensor coordinate system of the j th camera xo s j = xo s j yo s j Measured position of the reference subset in the ideal sensor coordinate system of the j th camera xd s j = xd s j yd s jMeasured position of the investigated subset in the ideal sensor coordinate system of the j th camera xs j = xs j ys j T Measured coordinate in the ideal sensor coordinate system of the j th keypoint location in the first image of the image series of the j th camera ξ x j and ξ y j Scaling of metric units to pixels for the j th camera ζ Maximum distance, along a single axis, from the center of the reference subset to a pixel in the reference subset

Table 1 .
Required input variables of the ADIC3D framework.

Table 2 .
Accessing the stereo matching data contained in ResultData.Stereo for subset q.

Table 3 .
Accessing the processing data for temporal matching (ResultData.ProcData1(d) and ResultData.ProcData2(d)) for image d and subset number q.
Call subroutine ImgCorr to perform temporal matching of first image series; Line 17 Call subroutine ImgCorr to perform temporal matching of second image series; Line 18 Call subroutine CSTrans to perform displacement transformation from the distorted sensor CSs to the world CS;

Table 17 .
Percentage difference of computed displacements between ADIC3D and DICe.
Gaussian filtering window size c x j and c y j Translation applied to sensor coordinate system of the j th camera c s j Skew of ideal sensor coordinate system of the j th camera C ZNSSD Zero-mean normalized sum of squared difference correlation criterion C ZNCC Zero-mean normalized cross correlation criterion C ObjFun Objective function of correspondence problem D Geometric error cost function δ Function to determine the Euclidean distance between two points e s j = e 1