3D LEAST SQUARES MATCHING APPLIED TO MICRO-TOMOGRAPHY DATA

The paper introduces 3D least squares matching as a technique to analyze multi-temporal micro-tomography data in civil engineering material testing. Time series of tomography voxel data sets are recorded during an in-situ tension test of a strain-hardening cement-based composite probe at consecutive load steps. 3D least squares matching is a technique to track cuboids in consecutive voxel data sets minimizing the sum of the squares of voxel value differences after a 12-parameter 3D affine transformation. For a regular grid of locations in each voxel data set of the deformed states, a subvoxel-precise 3D displacement vector field is computed. Discontinuities in these displacement vector fields indicate the occurrence of cracks in the probes during the load tests. These cracks are detected and quantitatively described by the computation of principal strains of tetrahedrons in a tetrahedral mesh, that is generated between the matching points. The subvoxel-accuracy potential of the technique allows the detection of very small cracks with a width much smaller than the actual voxel size.


INTRODUCTION
In materials research, several measurement techniques are used. Classical instruments, that allow the observation of surfaces of specimens, are for example strain gauges, inductive displacement transducers or inclinometers. In the last years, photogrammetric camera systems also became more and more part of the facilities of research labs with the advantage of being contactless and having a high spatial resolution and accuracy. Most photogrammetric contributions in this field use image sequences of camera systems and apply digital image correlation (DIC) techniques in order to compute displacement vector fields and use them for further analysis (Hampel and Maas, 2003, Maas and Hampel, 2006, Hampel and Maas, 2009, Sutton et al., 2009, Barazzetti and Scaioni, 2010, Koschitzki et al., 2011, Liebold and Maas, 2016, Liebold and Maas, 2018, Liebold et al., 2019, Liebold and Maas, 2020, Liebold et al., 2020a, Liebold et al., 2020b. Early applications of DIC were based on the cross-correlation method (Barnea and Silverman, 1972). Later, gradient based techniques were developed (Lucas and Kanade, 1981). These techniques used an iterative least squares algorithm to compute the displacements. (Ackermann, 1984) and (Grün, 1985) extended the mathematical model to an affine transform where rotation, scaling and shearing were considered. However, DIC techniques applied to image sequences only allow the observation of a specimen's surface. A further interesting point is the view inside the specimen. Therefor, X-ray tomography can be used. During load tests, a specimen can be observed in an in-situ experiment where different states of deformation are recorded in a tomograph. The above- * Corresponding author mentioned 2D techniques can be extended to 3D and are then called digital volume correlation (DVC). (Bay et al., 1999) applied 3D cross-correlation to X-ray tomography data in order to compute 3D displacement vector fields. While 3D crosscorrelation determines three shift parameters between consecutive cuboids of voxel data, (Maas et al., 1994) described a subvoxel-precise 12-parameter 3D least-squares cuboid tracking applied to a sequence of 3D images taken from two mixed fluids. The mathematical model was an affine transform that included translation, rotation, scale as well as shear parameters. The paper at hand is based on the work of (Maas et al., 1994) and will extend the mathematical model with radiometric parameters. The algorithm will be applied to data sets of different deformation states from an in-situ tension test with a fiber-reinforced composite probe (strain-hardening cementbased composite, SHCC). The data was recorded by a micro-CT instrument and the experiment is described in (Lorenzoni et al., 2020). Fig. 1 shows one slice of the voxel data of the reference state as well as one slice of a deformed state from the experiment. In the deformed state, cracks are clearly visible.
The paper is structured as follows: Section 2 gives an overview of the 3D least squares matching algorithm. Then, section 3 shows the analysis of the displacement fields using the strain analysis. At the end, a conclusion is given.

Mathematical Model
The 3D least squares matching (3D LSM) algorithm is applied to two volume data sets. In the reference volume, the integer The aim is to find the corresponding subvolume with the coordinates x Def , y Def , z Def in the second volume of the deformed state. Eq. 2 shows the relationship between the gray values of the first (reference) and second (deformed) volume, also taking into account radiometric parameters (r0, r1): Eq. 2 is valid for all voxels in the cuboid. Between the two states, a coordinate transformation is applied. Similar to (Ackermann, 1984) and (Maas et al., 1994), a 3D affine transformation is used containing displacements, rotations, shear and scaling, see Eq. 3. The shifts are represented by a0, b0 and c0. x where ai, bi, ci = affine parameters Thus, the vector of unknowns is: p has 14 unknowns and the cuboid contains at least 3 × 3 × 3 = 27 voxels (usually more), so that there is an over-determination. Thus, the parameters of p are computed using the least squares method. Therefore, a residuum is added to Eq. 2 in order to get the observation equation: The linearization of the observation equation is done using the first terms of the Taylor series: where x Def,0 = x computed with initial values of ai y Def,0 = y computed with initial values of bi z Def,0 = z computed with initial values of ci r0,0, r1,0 = initial values of r0, r1 The gray value g(x Def,0 , y Def,0 , z Def,0 ) is computed by trilinear interpolation at the position (x Def,0 , y Def,0 , z Def,0 ). Eq. 6 also contains derivations: the volume gradients of the gray values.
They can be obtained by a combination of numerical differentiation gx,i, gy,i, gz,i at the integer positions xi, yi, zi (central differences, Eq. 8) and tri-linear interpolation.
gx,i ≈ 0.5 · (g(xi + 1, yi, zi) − g(xi − 1, yi, zi)) gy,i ≈ 0.5 · (g(xi, yi + 1, zi) − g(xi, yi − 1, zi)) gz,i ≈ 0.5 · (g(xi, yi, zi + 1) − g(xi, yi, zi − 1)) The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B2-2021 XXIV ISPRS Congress (2021 edition) The differentials dx Def , dy Def and dz Def are computed by: Thus, the vector of corrections to the unknowns is: The vector of unknowns is composed of the vector of the initial values p 0 and the vector of the corrections to the unknowns: The linearized observation equations can be written in matrix notation: where A = Jacobian matrix l = reduced observation vector v = residual vector The computation of the Jacobian matrix is shown in Eq. 13 and Eq. 14 shows the reduced observation vector.
In the Gauss-Markov model, the sum of the squares of the residuals is minimized: To achieve this, the normal equations are solved for dp: The corrections to the unknowns are added to the initial values using Eq. 11. The steps from Eq. 13 to Eq. 16 (and p = p 0 +dp) are repeated until the process converges.

 
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B2-2021 XXIV ISPRS Congress (2021 edition) The standard deviation of the unit weight is: where v = A · dp − l = residual vector n = number of observations u = number of unknowns In addition, the standard deviation of the i th parameter is: where Q xx = (A T · A) −1 = cofactor matrix of unknowns

Initial Values
Due to the non-linearity of the gray value distribution in the volume data, initial values have to be obtained for the 3D LSM algorithm if the movements between the epochs exceed the dimensions of the cuboid. In case of predominant translations and small rotations, 3D cross-correlation can be used to compute initial shifts.

Application in the Experiment
As explained in the introduction part, the 3D least squares matching method is applied to a sequence of voxel data sets of a reference and six consecutive load levels. First, a regular 3D grid of points is defined, and points with insufficient contrast in their neighborhood are excluded. Before applying the least squares method, initial values are computed using 3D cross-correlation (Bay et al., 1999). Then, the displacement vector fields are computed using the 3D least squares matching algorithm. Fig. 2 shows a sequence of six displacement vector fields between the reference state and the deformed states. With higher loads (higher steps), cracks appear and grow. This also leads to higher differences between parts of the displacement vector fields.
To evaluate the influence of the affine and the radiometric parameters, the algorithm is applied introducing 3 (only translation), 12 (affine transformation) and 14 parameters (affine transformation and radiometric parameters) in the least squares process. The mean standard deviations of unit weight (Eq. 17) of the matching results are compared for the three versions and for each epoch, see Table 1. The differences in the s0 values are less than 2 %. Between the epochs in the experiment, there are only very small rotations and also the radiometric conditions are almost the same such that there is no significant advantage to use 12 or 14 parameters. However, for other experiments, it may improve the results.
In addition, the precision is analyzed:

DEFORMATION ANALYSIS
In our material testing approach, it is important to find the discontinuities in the 3D displacement vector field. The computation of strains allows to identify such areas. The method that will be presented in Sec. 3.1 requires at least four points and the corresponding shifts to calculate a strain. Because of that, the points of the displacement vector field are triangulated into a tetrahedral mesh and for each tetrahedron, strains are computed. Fig. 3 shows a tetrahedral mesh of the probe of Fig. 1.

Strain Analysis
This section gives a short overview of the calculation of principal strains as well as volume strains in tetrahedral meshes as an extension of the 2D triangle analysis approach presented in (Liebold and Maas, 2016). The principal strains of tetrahedrons are computed using the coordinates of its vertices in the reference (x Ref ,y Ref ,z Ref ) and in the deformed state (x Def ,y Def ,z Def ). First, the parameters of an affine transformation between the coordinates are calculated using the four vertices of the tetrahedron.
where ti,fij = affine parameters The deformation gradient F is composed of the nine parameters fij. F can be decomposed into a product of a symmetric and rotation matrix (Becker and Bürger, 1975).
where F = deformation gradient tensor R = rotation matrix V = left stretch tensor (symmetric) The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B2-2021 XXIV ISPRS Congress (2021 edition) Step 0→1 Step 0→2 Step 0→3 Step 0→4 Step 0→5 Step 0→6  The left Cauchy-Green deformation tensor V 2 is calculated in order to compute the polar decomposition: In the next step, an eigenvalue decomposition of left Cauchy-Green deformation tensor is performed: where C = eigenvector matrix (orthogonal matrix) Λ = eigenvalue matrix (diagonal matrix) λi = i th eigenvalue, diagonal element of Λ and λ1 ≤ λ2 ≤ λ3 The greatest principal strain 3 (technical strain) is derived from the eigenvalue λ3 (Eq. 23). It is a dimensionless quantity.
From F, it is also possible to derive a volume strain (Ogden, 1997, Altenbach, 2018: A crack crossing a tetrahedron will cause an extension of the tetrahedron and will thus lead to a larger value of 3 as well as v . The corresponding eigenvector (column of C) gives the direction of the strain 3. Strain values greater than 0 indicate a tension whereas a value of 0 indicates a stable element.

Application to the Experimental Data
In Fig. 4, the s3 strains (Eq. 23) are visualized for the six different load steps (cross sections of the voxel data blended with a cross section of the color-coded tetrahedrons). The evolution of cracks is visible. The corresponding stress-displacement curve is shown in the work of (Lorenzoni et al., 2020).
Especially, the visualization of step 5 shows some fluctuation in the strain values at the lower crack due to varying sizes of the tetrahedrons. It is a consequence of some missing matching results for points with matching cuboids crossed by the crack where LSM fails. In the lower center, some artefacts in the strain map are visible due to some matching errors at the air hole.

CONCLUSION
The paper presents the algorithm and application of 3D least squares matching to crack detection in multi-temporal microtomography voxel data sets. The affine parameters in the model allow to regard rotations and scaling between states. Changes in illumination are considered by radiometric parameters. Displacement vector fields are computed with standard deviations of the translation parameters of less than 0.02 vx. Strain values are calculated from the displacements to show deformed areas and to detect cracks.
Future work could concentrate on the determination of the accuracy of the 3D least squares matching.