A global digital volume correlation algorithm based on higher-order ﬁnite elements: Implementation and evaluation

We propose a DVC technique that is based on higher-order ﬁnite-element discretization of the displacement ﬁeld and a global optimization procedure. We use curvature penalization to suppress non-physical ﬂuctuations of the displacement ﬁeld and resulting erroneous strain concentrations. The performance of the proposed method is compared to the commercial code Avizo using trabecular bone images and found to perform slightly better in most cases. In addition, we stress that the performance of a DVC method needs to be evaluated using double scans (zero strain), virtual deformation (imposed deformation) and real deformation. Double scans give insight into the presence of noise and artifacts whereas virtual deformation benchmarks allows evaluation of the performance without noise and artifacts. Investigation of the performance for actual deformed heterogeneous materials is needed for evaluation with noise, artifacts and non-zero strains. We show that both decreasing the resolution of the displacement ﬁeld (increasing subvolume size) as well as (increasing) curvature penalization (regularization) have a similar effect on the performance of evaluated DVC methods: Decreasing the detrimental effect of noise, artifacts and interpolation errors, but also decreasing the sensitivity of a DVC method to displacement peaks, discontinuities and strain concentrations. The needed amount of regularization is a trade-off between accuracy and precision of the estimated strain ﬁelds and their resolution. The obtainable accuracy and precision of the estimated displacement ﬁelds are inﬂuenced by interpolation errors in the DVC procedure and the relative amount of detail, noise and artifacts in the images. Errors in the displacement ﬁeld are typically magniﬁed during the strain calculation. Based on the tests and subvolume sizes (16–50 voxels) in this study, the expected order of magnitude of the accuracy and precision is 0.1 micro-voxels and 1 milli-voxels for the displacements and 0.1 and 1 milli-strains of the strain ﬁelds. © 2019 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
High-resolution Computed Tomography (CT) imaging has made it possible to obtain three-dimensional images of materials at the microscale.Sequential images taken at different times during deformation allow for estimation of local deformation fields using Digital Volume Correlation (DVC) ( Bay et al., 1999 ), a logical extension of Digital Image Correlation techniques (see for instance the review paper by Pan et al. (2009) ) to three-dimensional data.Fast acquisition using synchrotron radiation allows obtaining CT images at time intervals of only a few seconds and tracking the deformation mechanisms closely over time.Although significant steps for-  DVC and DIC techniques compare sets of images, a reference image and a deformed image, to find matching features and estimate the corresponding displacement field.The deformation field, i.e. strain, is then typically obtained in a post-processing step by calculating the gradient of the displacement field.To estimate the displacement field, an image is typically divided into many (potentially overlapping) subvolumes.For each subvolume, a matching subset in the next image (and its corresponding displacement) is sought that maximizes their correlation or minimizes their differences in gray scale values.There are several different ways to perform this task.The accuracy (mean error) and precision (standard deviation) of the resulting strain fields depend greatly on the choice of algorithm and algorithmic parameters such as subvolume size ( Liu and Morgan, 2007;Dall'Ara et al., 2014 ).
Images can be compared directly using Fast-Fourier Transform (FFT), which results, in its basic form, in an integer translation for each subvolume and does not account for local rotation or deformation (also called rigid registration) ( Forsberg et al., 2010 ).It is, however, possible to obtain sub-voxel accuracy and treat large deformation problems in an iterative FFT-based algorithm based on incremental deformation of the images and filtering techniques ( Bar-Kochba et al., 2015 ).
For sub-voxel accuracy and displacement fields including deformation (also called elastic registration), it is more common to formulate the comparison of images as a nonlinear minimization problem that needs to be solved iteratively as first done by Bay et al. (1999) using only translational degrees of freedom.For improved accuracy and precision, additional degrees of freedom need to be added to the subvolume displacement field.A six degree of freedom model was developed by Smith et al. (2002) , adding three rotations to the description.A further six degrees of freedom, thus twelve in total, are necessary for a full trilinear interpolation of the displacement field, see e.g.Gates et al. (2011) .An initial guess for the displacement field can be obtained using an integer-displacement FFT-based technique.
The previously described methods for comparing images are typically used on the subvolume level and can be characterized as local ( Roberts et al., 2014 ).A local approach uses independent parameters for the displacement interpolation of each subvolume.These parameters are then obtained in an iterative procedure for each subvolume separately, see e.g. ( Forsberg et al., 2010 ).As a consequence, it is possible that the identified displacements at the boundaries of subvolumes (and also in their interior) are completely different, leading to inaccurate and ambiguous results, see Fig. 1 .Alternatively, a global approach uses a globally shared set of parameters for the displacement interpolation and a solution procedure that uses information from all subvolumes simultaneously ( Madi et al., 2013;Roberts et al., 2014 ).Due to the continuity of the global displacement interpolation, the accuracy of the results of global DVC methods is generally higher than that of local ones ( Dall'Ara et al., 2014 ).
We would like to emphasize the connection between the displacement and the strain fields.The obtained displacement field from a DVC technique needs to be differentiated spatially in order to obtain the strain field.For local DVC approaches, a global displacement field needs to be constructed based on the information on the subvolume level.Often, the displacements are interpolated and smoothed before translation into strain fields, see e.g.Forsberg et al. (2010) .Note that, if displacement gradient information is available on the subvolume level, it is dubious to directly interpolate the displacement gradient due to potential discontinuities in the displacement of neighboring subvolume, see Fig. 2 .For global approaches, the displacement field can be differentiated directly to obtain the strain.Note that a nonlinear strain definition such as the Green-Lagrange strain may be necessary for finite deformation and rotation.
Similar to DIC techniques, the quality of the results of DVC methods generally depend on the amount of detail in the source images.In regions with little detail, there is not enough information to reliably determine a displacement field.If unchecked, this can lead to large errors in the calculated results.To improve the results, different kinds of regularization can be applied ( Leclerc et al., 2012;Barber and Hose, 2005 ).A common way to improve the performance is to remove the results of subvolumes with low correlation ( Gillard et al., 2014 ) or outliers ( Bar-Kochba et al., 2015 ).
Anomalies in the source images also have a major effect on the performance of DVC techniques.These anomalies include intrinsic noise and potential artifacts such as beam hardening ( Barrett and Keat, 2004 ).They represent the other extreme; too much (nonphysical) detail.If no measures are taken, then a DVC procedure may overfit the data and generate a non-physical displacement field.Regularization can also help alleviate problems in this case.
Another error source is the interpolation of gray-scale values that is inevitable for sub-voxel translations, rotations and deformation ( Madi et al., 2013 ).For integer-translations, one can compare the gray-scale values in a CT image set directly.In all other cases, interpolated values are needed.Again, regularization can be used to try and alleviate the effect on the performance of DVC techniques.This effect will be investigated in the present paper.
Regularization is usually aimed at smoothing the displacement field to remove non-physical strain concentrations.Barber and Hose ( 2005) apply smoothing based on the Laplacian operator that, in fact, reduces the curvature of the displacement field.One should be aware that, in particular for large amounts of regularization, this also affects the displacement field in regions where there is physical deformation.This will reduce the performance of DVC techniques in regions with strain (and stress) concentrations.These effects are illustrated in Fig. 3 .
For validation of DVC techniques the accuracy and precision of two types of benchmark tests have been used.Both have their advantages and disadvantages.
• Double scans: The same object is scanned twice, with or without repositioning, see for instance the work by Dall'Ara et al. (2014) .This test includes the effect of noise and artifacts but does not include actual deformation.In the presence of regularization that directly minimizes variations in the strain field ( i.e. curvature penalization of the displacement field), this type of benchmark needs to be supplemented with non-zero deformation benchmarks for a more realistic evaluation of the performance.• Virtual deformation: A known translation, rotation, deformation or a combination of these is imposed on a CT image, see for instance the work by Bar-Kochba et al. (2015) .It is particularly useful to be able to examine the performance of a DVC technique for varying strains.Unless manually added to the data, this test does not include noise or artifacts allowing for very good accuracy and precision.It is also possible to use a double scan for the construction of the deformed image, introducing more realistic noise and artifacts.Note that the interpolation necessary for creating the deformed image constitutes one of the error sources.
A thorough evaluation should include both types of benchmarks to investigate the performance in the presence of noise and artifacts and non-zero strains.
It is not yet clear what choices of implementation lead to the best performance for DVC.In particular, the accuracy and precision of the calculated displacement and deformation fields depend on the chosen techniques and some algorithmic parameters.Furthermore, the evaluation of the performance of a DVC method does not always include both double scans and virtual displacement fields in the literature, which renders any comparison inconclusive in the opinion of the authors.In addition to presenting a novel DVC approach, this study therefore aims to advance the understanding of the techniques used by investigating the effect of, in particular, the gray-scale interpolation and the amount of regularization.We investigate the order of magnitude of the performance measures based on the different error sources and realistic resolutions of the displacement-field estimation.
Additionally, we propose a DVC technique that features an increased number of degrees of freedom of the displacement field per subvolume and a global optimization with curvature penalization.The global displacement field is interpolated using higher- order (triquadratic) finite elements, also called 27-node bricks.This choice is also convenient for the implementation of overlapping subvolumes for more robust solutions in a global DVC method.After obtaining an initial guess for the displacement field using a local FFT-based approach, the elastic registration problem is solved in a global optimization procedure where the nodal displacements of the finite-element discretization are the design parameters.Additionally, curvature penalization can be used for regularization purposes.Finally, the global displacement field is processed into a global strain field.First, a globally smooth spline interpolation for the displacement field is constructed using all degrees of freedom after which the corresponding finite strain field can be calculated (valid for finite deformation and rotations).
Reliably estimating displacement fields with a high resolution remains challenging for trabecular bone since the struts of the trabecular bone are large compared to the desired resolution.Therefore, the performance of the presented DVC method is tested on a number of different benchmark problems using two sets of CT images of trabecular bone including an actual compression test.For illustration, a slice of the reference volume of these data sets is shown in Fig. 4 .
• Tozzi-2017: This is a relatively low resolution set of images (Specimen-1-VOI4) of the ones used in the work of Tozzi et al. (2017) .The complete size of this set of images is 152 × 152 × 432 of which only a part, see Fig. 4 a, was used in all analyses except the double scan.Due to its relatively low resolution, the features in this set of CT images are small compared to the size of the voxels, which is good for relative accuracy and precision of the estimated displacement field (when measured in voxels) and strains.Our results have not been compared to those of Tozzi et al. (2017) since we do not have access to their implementation or specimen-specific results.This data set is used for the following benchmark problems.
-Virtual Gaussian displacement field: non-zero strain -Virtual sub-voxel translation: zero strain -Double scan (no repositioning): zero strain • UU-2018: This is a new set of high-resolution CT images of a trabecular bone sample obtained at the TOMCAT beamline at Paul Scherrer Institut (PSI), Villigen, Switzerland.The complete sample size is 7 mm in diameter and 11 mm in height, but only part of this has been used, see Fig. 4 b.Due to its high resolution, the features of this set of CT images are relatively large compared to the size of the voxels.For this reason, relatively large subvolumes will be necessary for good accuracy and precision.Additionally, the data also contains significant artifacts which will affect the performance.This data is used in combination with the following benchmarks.
-Double scan (no repositioning): zero strain -Real deformation: non-zero strain First we introduce the methodology of the proposed DVC approach in Section 2 followed by the results of the proposed method in Section 3 .This includes a comparison with a state-of-the-art DVC technique and a discussion of the effect of the subvolume size and the amount of regularization.This work ends with the conclusions and recommendations in Section 4 .

Theory
The proposed DVC approach, which we will refer to as GDVC-UU, uses a global interpolation of the displacement field using a finite-element mesh of overlapping 27-node brick elements.In a preprocessing step, an integer-valued initial displacement field is obtained using a local DVC algorithm based on phase correlation and a Fast-Fourier Transform (FFT), similar to for instance the work of Forsberg et al. (2010) .The nodal degrees of freedom of the displacement interpolation are then sequentially updated in a nonlinear optimization by maximizing the normalized cross-correlation between the gray scale values of the reference (before the deformation increment) and current (after the deformation increment) volumes.Instead of comparing the entire reference and current volumes at once, the procedure is divided into multiple comparisons of overlapping subvolumes to increase the sensitivity to local differences.The resolution of the global displacement field is increased stepwise by refining the finite element mesh and using the interpolated results of the previous coarser mesh as the initial displacement field.To suppress non-physical deformation due to noise, artifacts and/or too little detail in the images, a curvature penalty can be used as regularization.After the nodal displacements have been obtained, a unique displacement field is constructed using a global spline interpolation and the corresponding strain fields can be calculated.This procedure is illustrated in Fig. 5 .

Gray-scale and displacement interpolation
The estimation of local deformation using CT images typically consists of a sequential comparison of three-dimensional images taken at certain time intervals.In every step two images are compared: A reference (before deformation increment) and a current (deformed) volume.The corresponding input consists of threedimensional arrays of gray-scale values, one for the reference volume, S IJK , and one for the current volume, s ijk , where ( I, J, K ) and ( i, j, k ) are the corresponding indices or voxel coordinates.In the proposed approach, corresponding features and the related incremental displacement field are simultaneously identified by maximizing the normalized cross-correlation differences between the gray scale values of current volume and the transported gray scale values of the reference volume.
For such a comparison, we interpret the point-wise voxel information as regular samplings of continuous gray scale functions S ( X ) and s ( x ), where X and x are the reference and current coordinates, respectively.We can interpolate these functions based on the gray-scale values of the input data as,

S( X
where X IJK and x ijk are the coordinates of the centers of the voxels.In the current paper, third-order spline interpolation (locally quadratic) has been used to interpolate the gray scale values.
Following common notation in continuum mechanics, the reference and current coordinates are related by the displacement field u as,  where the Lagrangian description has been used (in terms of the reference coordinates X ), see Fig. 6 .
This displacement field needs to be described by a set of parameters which will serve as the design variables in an optimization.In this paper, we employ a global finite-element interpolation using standard isoparametric finite elements.Consequently, the displacement at an arbitrary position in the reference volume X can be expressed as, where n is the number of nodes of the finite-element mesh and N ki and q i are the corresponding finite-element basis functions and nodal displacement vectors ( i.e. the design variables), respectively.Note that finite-element basis functions typically are sparse, i.e.only nonzero in the direct neighborhood of the corresponding nodes.This allows us to write the same relation for the displacement of a point in the reference volume in terms of only the nodal displacement vectors in the containing finite element/subvolume as, where q e is a vector with 3 n e components, N is a 3 × 3 n e displacement interpolation matrix and n e is the number of nodes of a finite element/subvolume.In the actual implementation, this relation is evaluated for many position vectors at the same time for efficiency.
For a given displacement field u ( X ), that means a given set of nodal displacements q , the interpolated gray scale functions can be used to compare the reference and deformed volumes voxel by voxel.The most straight-forward way to do this is by calculating the current coordinates of all voxels of the reference domain from the relation, where u IJK ( q ) = u ( X IJK ) which depends on the nodal displacements q .Then the corresponding gray-scale values of the current volume are obtained as, These values can then be directly compared to the gray scale values of the reference volume.This procedure is illustrated in Fig. 7 .

Global optimization
The most direct approach to estimate the displacement field is a global minimization of the differences between the gray-scale values of the reference volume and the interpolated gray-scale values of the deformed volume as a function of the displacement field, see for instance the work of Pan et al. (2012) .It is convenient to replace the three indices, I, J and K of the three-dimensional data by a single index, α, so that, S α = S IJK and A direct minimization can then be formulated as a least-squares problem in terms of the vectors S and s ( q ) as, where • 2 indicates the Euclidean norm.This minimization problem can be solved directly, see e.g. the works of Pan et al. (2012) and Madi et al. (2013) .However, instead of directly minimizing the differences in Eq. ( 8) , it is more common in DVC and DIC approaches to use normalized cross-correlation for comparison ( Roberts et al., 2014 ).Introducing the following short-hand notation of the normalized gray scale values of the reference and current volumes, ˆ S = S S and ˆ the normalized cross-correlation, C , is defined as, The normalized cross-correlation is a direct measure of how much the vectors S and s have the same direction.If C = 1 then S and s have exactly the same direction.If C = 0 (or C < 0) then they are perpendicular to each other (or even in opposite directions).The normalized cross-correlation is, however, insensitive to the lengths of the vectors S and s .
Using normalized cross-correlation, the displacement identification can be formulated as a global minimization problem in the following fashion, min q { −C( q ) } (11) In our experience, using the normalized cross-correlation as an objective in an optimization leads to faster convergence than a direct minimization of the differences between gray-scale values.
Instead of using Eq. ( 11) , the global minimization is reformulated as minus the average normalized cross-correlation of all subvolumes to take advantage of the finite-element discretization and increase the sensitivity of the optimization to local differences, min where C e ( q e ) is the normalized cross-correlation for a subvolume e calculated analogously as in Eq. ( 10) , q e is the corresponding element displacement vector and N is the total number of subvolumes.Note that the average correlation of all subvolumes does not exactly equal the correlation of the entire volume, see also the right graph in Fig. 10 which would otherwise be a straight line.The optimization problem is solved using standard unconstrained optimization software.

Regularization
To suppress the effect of intrinsic noise and artifacts, as well as alleviate potential conditioning problems due to too little detail in certain areas of the source images, regularization can be added to the optimization process.Here, we have included the following curvature penalty for the displacement field in the objective function per subvolume, i.e. , P e ( q e ) = where e is the domain of a subvolume.The main integral in this equation is evaluated using standard 3 × 3 × 3 Gauss integration.The magnitude of the penalty depends on the element displacement vector q e through the displacement field, For non-overlapping subvolumes, this formulation results in a constant amount of regularization for refinement of the finiteelement-based displacement interpolation.
The objective function including regularization then becomes, min where p is a penalty factor for adjusting the amount of regularization.One can attempt to relate the magnitude of this penalty factor to the wave length of spurious oscillations in the resulting displacement field, see the work of Leclerc et al. (2011) , or minimize the condition number of the Hessian matrix of the optimization, similar to the work of Barber and Hose (2005) .

Implementational aspects
In the proposed DVC approach we have selected a finiteelement mesh with standard 27-node brick elements.For this type of finite element, the displacement field is interpolated using triquadratic polynomials ( i.e. the terms with highest order are x 2 y 2 z 2 ) allowing for more freedom in the displacement field.A conventional finite-element discretization does not include overlap and has a unique definition of the displacement everywhere in the volume.However, the usage of 27 nodes allows an easy implementation of overlapping subvolumes which, in the authors' experience, increases the robustness of the approach.
Maximum overlap, i.e. , overlaying meshes that are displaced half a subvolume size in all combinations of the three spatial directions, would result in approximately eight times as many subvolumes and corresponding computational effort.Instead, we use a limited overlap with approximately twice as many subvolumes by only adding a second mesh that is displaced half a subvolume size in all three spatial directions, as illustrated in Fig. 8 a.Using overlap implies that the displacement is uniquely defined at the nodes, but not in between nodes of overlapping subvolumes, see Fig. 8 b.This can help prevent overfitting to noise and artifacts in the source images.In a postprocessing step, a single unique global displacement field is constructed from the calculated nodal displacements using a third-order spline interpolation.
To increase the sensitivity of the normalized correlation, it is very important to center the data sets around zero.If one would omit a centering procedure, the correlation between images with all-positive gray-scale values will always be positive, even for visually 'uncorrelated' images, see Fig. 9 .For the current application, it makes most sense to use the same average gray scale value for centering for the two data sets, since the gray scale values of consecutive scans were obtained using the same parameters in the reconstruction process.Therefore, we use the following procedure, where S indicates the average of S .For other applications, it may be more suitable to use the individual averages of the two data sets separately.It can also be advantageous to center the subvolume data sets to increase the sensitivity to differences between the images even more.Note however, that this may lead to errors for small sub-  volumes in the presence of significant noise.This is illustrated in Fig. 10 .For decreasing subvolume size, noise in the data sets causes the average correlation to quickly go down when using centering of subvolumes.This in turn could lead to undesired behavior in the optimization process.
Termination criteria are needed for the optimization process, which can be based on the decrease in objective value, first-order optimality (magnitude of the gradient vector) or change in parameters.In this work we formulated the termination criterion in terms of a change in parameters as, where N is the total number of degrees of freedom.
For good results, a good initial displacement field is essential.It is well known that for larger subvolume sizes, the solutions found by DVC techniques become more reliable.Therefore, it makes sense to increase the resolution of the displacement field (decrease the subvolume size) successively, see Fig. 5 .In this case, a refined finite-element mesh is constructed after completion of the global optimization of a lower resolution displacement field.The results can then be transferred from the old mesh by straight forward interpolation.
It can also be advantageous to exclude a narrow band around the edges of the reference image to avoid issues with missing data in the current volume.This option has also been implemented in GDVC-UU.

Global displacement field and strain post-processing
The procedure to estimate local deformations from CT images typically consists of a sequential comparison of three-dimensional images taken at n time intervals, see Fig. 11 .For a time interval i , two consecutive images are compared, the reference volume at time t i −1 , and the current volume at time t i .
For each time interval, the global optimization produces an estimated displacement field that is characterized by the nodal displacements q i , see Eq. (2) .To obtain a unique global displacement field in the case of overlapping subvolumes, a global thirdorder spline interpolation is constructed for the current coordinates x i ( X i ) and the corresponding displacement field u i ( X i ) based on the computed nodal displacements q i , see also Fig. 8 b.
Assuming that the images are taken at the same coordinates, i.e. coinciding references coordinates X i , the current coordinates at the last time interval x t n of arbitrary reference coordinates at the beginning of the first time interval X t 0 can be computed  recursively as, and the total displacement field becomes, This allows the computation of the total displacements of every voxel in the reference configuration u tot, α .
Finally, the strain fields can be calculated in a post-processing step by taking the spatial derivatives of the displacement field.The most common strain measure used in finite strain theory is the standard Green-Lagrange strain E defined as, where In practice, the necessary spatial derivatives can be calculated by finite differences of the displacements of all voxels in the reference volume.

Results and discussion
In this section, we present the results of the benchmark problems for the proposed method GDVC-UU.We compare its results to that of the commercial high-performance 3D visualization and analysis software Amira-Avizo, in particular the included global DVC method based on the work of Madi et al. (2013) .We will refer to this method as Avizo.
First we demonstrate that GDVC-UU has excellent performance for virtual displacements.In this case, the data can be considered perfect and the only error source is the interpolation of gray-scale values giving an idea of a DVC technique's best case performance.Then we apply both DVC procedures on double scans illustrating the effect of the other error sources, such as the background (void), intrinsic noise and artifacts.In that context, the need for regularization is emphasized.Finally, the results are presented for a CT image including real deformation.For all benchmark tests we use one of the two data sets introduced in Section 1 .Here we start with providing some details regarding the evaluation of the results of the different DVC techniques and their comparison.
A common measure for comparing results is the accuracy acc , defined as the mean error, and the precision prec , defined as the sample standard deviation of the error, where a α is the voxel-wise quantity of interest (for instance displacement in x -direction) and n is the number of voxels.As a measure of the total accuracy and precision of the displacement (vector) or strain (second-order tensor) the following definitions are used here.The measure for the total accuracy is chosen such that a single large error has a relatively large impact.
Global DVC methods provide a continuous description of the global displacement field.Therefore, the consistent strains can be obtained directly by differentiating the displacement field and applying Eq. ( 20) .This is the approach that we use in GDVC-UU.When comparing the strains provided by Avizo with those obtained by Eq. ( 20) , we found significant differences, probably due to additional inter-and/or extrapolation of nodal and/or elemental strains.Although this may be a good idea for practical applications, it complicates the comparison of the results.Instead, we compute the strain field for Avizo using Eq. ( 20) after obtaining the displacement gradients directly from the displacement field with standard finite differences.
A margin near the edges of the reference volume, in which no displacements are calculated, is used for both Avizo and GDVC-UU unless otherwise specified.All presented performance measures are calculated in exactly the same part of the volume such that displacement and strain results are available for both methods.

Virtual displacement fields
A set of CT images without noise or artifacts is ideal to test the performance of a DVC technique.Such a set can be created by imposing a predefined displacement field on any reference volume of choice.We use the Tozzi-2017 data set for a benchmark with a Gaussian displacement field and the UU-2018 set for a subvoxel translation displacement field.It is important to realize that imposing an arbitrary displacement field on a reference volume involves interpolation of the gray-scale values due to sub-voxel translation, rotation and stretching.This is different from the procedure presented in the work of Bar-Kochba et al. (2015) where particles are only displaced with a Gaussian displacement field.
The approach used here consists of two steps, illustrated in Fig. 13 .In the first step, the current coordinates of the centers of all 'reference' voxels are calculated based on the predefined displacement field u ( X ).In the second step, the gray-scale values are calculated of the regularly spaced centers of all 'current' voxels by linear interpolation.This linear interpolation is based on an un-Fig.13.Illustration of the procedure to impose a displacement field on a CT image.structured set of data points and, therefore, we use standard Delaunay triangulation.Note that the construction of the current volume introduces an error source, since the gray-scale values need to be interpolated.Exceptions are integer-voxel translations.

Gaussian displacement field
The Gaussian displacement field is a sum of four multidimensional Gaussian displacement field as, where all operations are component-wise and the constants used are listed in Table 1 .This particular choice was obtained by tuning the height and width of the individual Gaussians to obtain a contribution to the strain of around 15%.Note that these are very large deformations, much larger than one would expect in reality in the absence of strong material nonlinearities such as fracture.
The displacement field for this constructed data set was estimated for subvolume sizes of 48, 32 and 16 voxels using Avizo and GDVC-UU with different amounts of regularization p = 10 .0 , 1.00 and 0.00.In Table 2 and Fig. 14 the accuracy and precision of the resulting displacement fields are listed and presented in a bar plot.The heights of the columns indicate the accuracy (mean error) and the error bars the precision (standard deviation).
For a fair comparison, the differences between Avizo and GDVC-UU in the spacing of the nodes (or degrees of freedom) needs to be taken into account.For GDVC-UU, the node spacing is in fact half the subvolume size, see Fig. 8 a.For this reason, we compare Avizo with subvolume size 16 to GDVC-UU with subvolume size 32 and observe that GDVC-UU outperforms Avizo for this benchmark problem.This is also visible in Fig. 15 where a contour plot of a cross-section of the middle of the data is shown and in Figs. 17 a and b where the estimated displacement fields are compared with the exact solution along two lines A and B (see Fig. 15 ).The triquadratic basis functions of GDVC-UU appear to capture the Gaussian shape of the displacement field better than the linear basis functions used in Avizo ( Madi et al., 2013 ).We can observe some interesting trends in the results.For noise and artifact-free CT images, the accuracy and precision of the results increase with decreasing subvolume size for all methods.This is caused by the inability of a coarse interpolation of the displacement field to describe the locally high curvature of the Gaussian displacement field.This is clearly illustrated in Figs. 17 c and d.In this case, the accuracy and precision of the results can be improved by increasing the resolution of the displacement interpolation.This is the opposite of what is generally observed for double scans.
The regularization penalizes the curvature of the displacement field, i.e. reduces variations of strain.This effect is best illustrated by the results for GDVC-UU with subvolume size 16, see Fig 17 e  and f.Without regularization, the estimated displacement field captures the small Gaussian quite well.For increasing amounts of regularization, the curvature of the displacement field is pushed down leading to underestimation of the displacement peak and the corresponding strain.The effect of the regularization (up to p = 10 .0 ) on the accuracy and precision in Fig. 14 appears quite limited.However, the performance of the proposed method will definitely be affected for much larger amounts of curvature penalization as illustrated in the later examples.
An important observation concerns the accuracy and precision of the strain fields.Even thought the displacement field is captured quite accurately for a small subvolume size, the precision of the associated strain fields is still relatively low with a standard deviation of up to 1% strain.Since the strain is directly related to the spatial derivative of the displacement field, errors in the estimated displacement field are typically magnified during the strain calculations.Due to the continuity of the displacement field, the accuracy can still be high (small mean error), but this is definitely not the case for the precision (standard deviation of the error), see Fig. 16 .In other words, any fluctuations in the displacement field will lead to large errors in strain.Therefore, regularization is essential to suppress fluctuations as much as possible.Preferably, one would like to obtain at least a standard deviation of the errors of the strain field as low as 1%.Based on the results presented here, one should consider using a larger amount of regularization.

Sub-voxel translation
In the absence of noise and artifacts, the only error source present is the interpolation of the gray-scale values.The theoretical performance of the proposed method seems very good judging by the Gaussian displacement field benchmark.However, interpolation errors need to be investigated further using a benchmark that is closer to a worst-case scenario.Such a scenario is a subvoxel translation which introduces large interpolation errors near extremes of the gray-scale values, see Fig. 18 .Therefore, the next benchmark considers a constant translation of 2.5 voxels in all directions which is imposed on the reference volume of the Tozzi-2017 data set.Note that interpolation errors are introduced both when constructing the current volume by imposing a virtual displacement field on the reference volume, as well as in the actual displacement estimation of the DVC.
The proposed method is tested without regularization, i.e. p = 0 , giving an accuracy in the order of 100 micro-voxels and a precision in the order of 10 milli-voxels for the displacement field, see Table 3 .The accuracy and precision of the estimated strains are in the order of 100 micro-strains and 10 milli-strains, respectively.For this particular data set, these values can be interpreted

Table 3
The accuracy and precision of the calculated displacement and strain fields obtained by GDVC-UU for a virtual sub-voxel translation applied to the Tozzi-2017 data set.The results of normal analyses and an analyses where the exact solution was supplied as initial displacement fields are virtually identical.as the worst-case accuracy and precision of the estimated displacement and strain fields only due to interpolation errors (for perfect, artifact-free input data).This worst-case performance depends on the data set and subvolume size used.In a similar test case for the UU-2018 data set, slightly larger interpolation errors for the displacement field and similar values for the strain field were found.These orders of magnitude are listed in Table 4 .

GDVC-UU
One might think that these errors are the result of a bad solution (clearly suboptimal local minimum) in the optimization process.To check this, we included the results of the same analysis where the exact solution is given as the starting point of the optimization process leading to the same solution with the same accuracy and precision, see Table 3 .In other words, the average crosscorrelation of the final result of the optimization is higher than that of the exact solution.This clearly illustrates the influence of interpolation errors.
Regularization discourages oscillations in the displacement field and, therefore, suppresses interpolation errors.In Table 3 the results are included with the maximum regularization used before, i.e. p = 10 .0 .In particular, the precision (standard deviation) is improved for increasing amounts of regularization.The results suggest that even larger values of the penalty factor would be necessary for an imposed translation of 2.5 voxels in all directions.Indeed, this benchmark case really is a worst-case scenario and much better performance is achieved for the Gaussian displacement field.
The accuracy and precision also improve for increasing subvolume size, since the effect of interpolation errors on the crosscorrelation becomes lower.

Double scans
Besides interpolation errors in the DVC procedure itself, double scans also contain intrinsic noise and potential artifacts such as beam hardening.The Tozzi-2017 data set contains few clearly visible artifacts (as compared with the size of the features), whereas the UU-2018 data set contains clear ring artifacts, see Fig. 4 .Both data sets contain background noise, which is most apparent in void areas.Therefore, the background noise will be more problematic for the UU-2018 data set which has larger voids.In this subsection, we will compare the performance of Avizo and GDVC-UU for both data sets, several subvolume sizes and amounts of regularization.
The results of all methods agree on the average translation of both double scans.For instance, for the Tozzi-2017 the displacements found by all methods are within two hundredth of a voxel, see Fig. 19 .Also the standard deviation of the displacements  estimated displacement field.However, voids only contain background intensities and artifacts.During deformation, only the material domain is transported and deformed whereas the background intensities and artifacts are not.This can cause a discrepancy between the estimated displacement of neighboring material and void regions and lead to erroneous strains near the corresponding interface.This problem is most severe for large voids relative to the subvolume size.Due to these additional difficulties, larger steps in the amount of regularization have been chosen for the UU-2018 double scan.The order of magnitude of the accuracy and precision of the strains is similar to that of the virtual translation benchmark.
From Fig. 20 , it is clear that the accuracy and precision deteriorate for decreasing subvolume size (as opposed to the benchmark with the Gaussian displacement field).The same effect can be observed in Fig. 21 c and d: Decreasing the subvolume size increases the freedom of the displacement field which causes the DVC procedure to overfit to noise, artifacts and interpolation errors.Although not apparent from the benchmark problems, it is known that higher-order polynomial basis functions are more susceptible to overshoot than linear ones.From Table 5 , we observe that a large amount of regularization, more than p = 10 .0 for the Tozzi-2017 double scan and p = 10 4 for the UU-2018 double scan, is necessary to push the erroneous strains down to values in the order of 1%.This is also illustrated in Figs.21 e and f for the UU-2018 double scan.The problem is that increasing the regularization will also push down (variations in) the actual physical strains.For large amounts of regularization, the estimated strain fields are less likely to show large spatial differences in strain.For low amounts of regularization, the estimated strain fields will be dominated by large erroneous strains.Therefore, the amount of regularization that is needed is a difficult trade-off between accuracy and precision of the estimated strain fields and their resolution.
In most cases, the subvolume size and the amount of regularization have a very similar effect on the estimated displacement and strain fields.Increasing the amount of regularization and the subvolume size tend to decrease the amplitude of the displacement fluctuations and corresponding strains.
To illustrate the regularizing effect of introducing overlap in the displacement interpolation, a comparison of different amounts of overlap is included in Fig. 22 .For increasing amounts of overlap there is a decreasing trend in the accuracy and precision of GDVC-UU for the UU-2018 double scan.This effect is not very large though and introducing overlap significantly increases the amount of computational effort required to perform the displacement identification.For this reason, we only used limited overlap, i.e. , two overlapping finite element meshes, for the other results included in this paper, see also Section 2.4 .

Real deformation
The performance of DVC approaches for actual deformation is difficult to analyze in the absence of an exact solution.However, one can compare the results for different numbers of degrees of freedom and different amounts of regularization.The results are shown in Von Mises strain (a type of effective strain similar in definition to the Von Mises stress) and plotted for only material voxels (relative gray scale value higher than 0.30).Based on the experimental setup, an average vertical strain of around 1.5% is expected.
As mentioned before, relatively large regions of void lead to discrepancies in the estimated displacement field and erroneous strains near the void-material interfaces.For this reason, it can be advantageous to exclude the cross-correlation maximization for subvolumes that hardly contain any material.The results of GDVC-UU presented in this subsection have been obtained excluding subvolumes where more than 99% of the gray scale values are below 35 0 0 0 on a scale between 0 and 65 535.
In Fig. 23 , the resulting Von Mises strain fields of Avizo (calculated from the displacement field output) are shown for subvolume sizes (node spacings) 100 and 50.The tetrahedron displacement interpolation is clearly visible from the piecewise constant strain field.The results using subvolume size 100 and 50 are very  different.It is not immediately clear if any variations in real strain are suppressed in Fig. 23 a or if the strain concentrations in Fig. 23 b are non-physical and caused by error sources.There also appear to be edge effects causing non-physical strain concentrations of over 4% at the image boundaries.
In Figs.24 a and b, the corresponding results of GDVC-UU are shown for subvolume sizes 100 and 50 (node spacings of 50 and 25) and p = 100 .We observe similar behavior.In Fig. 24 a physical information may have been suppressed, but, more likely, strain concentrations of over 4% (present for both subvolume sizes) are non-physical.When increasing the amount of regularization to p = 10 4 , variations in strain are discouraged, see Figs. 24 c and d.The Von Mises strain fields for both subvolume sizes are now similar and much more realistic.This is supported by a rough estimate of the applied experimental displacement, which would result in an average strain in the z -direction around 1-2%.However, some physical strain concentrations may have been suppressed for this increased level of regularization.

Conclusions
The performance of a DVC method needs to be evaluated using double scans, virtual deformation (using a single scan or double scans) and real deformation.Double scans give insight into the presence of noise and artifacts whereas virtual deformation benchmarks using a single scan evaluate the performance for a known deformation and without noise and artifacts.Finally, using real deformation is important for investigating the performance for actual deformed heterogeneous materials.
In this paper, a DVC technique has been proposed that is based on higher-order finite-element discretization of the displacement field and a global optimization procedure.Curvature penalization has been used to suppress non-physical fluctuations of the displacement field and resulting erroneous strain concentrations.The proposed method is compared to the commercial code Avizo and performs slightly better in most cases, most likely due to the regularization.
For the Gaussian virtual deformation benchmark using a single image (noise and artifact-free images), regularization in fact reduces the performance of the DVC methods by suppressing variations in the strain fields.Decreasing the amount of regularization improves the performance until the interpolation errors start to dominate.Similarly, the performance is improved for decreasing subvolume size.The increased number of degrees of freedom allows a better approximation of the displacement field.
For the double scans (with noise and artifacts), the performance of the DVC methods improves for increasing but moderate amounts of regularization.The regularization suppresses erroneous fluctuations in the displacements and corresponding strain field caused by overfitting to noise and artifacts in the images.The accuracy and precision are influenced by interpolation errors in the DVC procedure and the relative amount of detail, noise and artifacts in the images.
The estimated strain fields are directly related to the spatial derivatives of the displacement field.Therefore, errors in the displacement field are typically magnified during the strain calculation.Based on the tests and subvolume sizes in this study, the expected order of magnitude of the accuracy and precision is 0.1 micro-voxels and 1 milli-voxels for the displacements and 0.1 and 1 milli-strains of the strain fields.
High-resolution images allow a detailed investigation of microand nano-structures, but analyses with both included DVC methods become more challenging as local features become more sparse (relatively large void regions) and CT images are more likely to contain noise and artifacts.For GDVC-UU much regularization is necessary to eliminate large fluctuations in the estimated displacement fields and corresponding large strains.The amount of regularization necessary is a trade-off between accuracy and precision of the estimated strain fields and their resolution.
Additional benchmark tests involving virtual deformation using a double scan are recommended to investigate the performance under non-zero deformation and in the presence of noise and artifacts.
where Eq. ( 10) has been used to simplify the expression.The partial derivative of the gray scale values of the current volume with respect to the displacement is obtained by differentiating Eq. ( 6) , The second-order derivative of the normalized cross-correlation with respect to the nodal displacements, the Hessian of the optimization with a minus sign, is more complex.(A.5)All the first-order partial derivatives in the Hessian have already been treated.After some manipulations, the second-order derivative of the normalized cross-correlation with respect to the gray scale values of the current volume can be expressed as, .6)where δ αβ is the Kronecker delta, equivalent to the an identity matrix.
In implementation it is not a good idea to build the large intermediate full matrix in Eq. (A.6) because of the excessive memory requirements.Instead one can obtain the first term of the significantly smaller Hessian in Eq. (A.5) directly by combining appropriate terms of the following form, where Q αi is defined as, 3 C t i t j − T i t j + t i T j − C Q γ i Q γ j (A.9) The second-order derivative of the gray scale interpolation with respect to the displacement is obtained by directly differentiating Eq. (A.3) one more time, .10)The second term in Eq. (A.5) can then be calculated in a straight forward manner.
The contribution of the regularization to the gradient and Hessian are obtained by differentiating Eq. ( 13) .The contribution to the gradient is,  .11)and the contribution to the Hessian becomes, d 2 P e d q i d q j = 2 p e d e 3 j,k,l=1 .12)Again, the main integrals in these equations are evaluated by Gauss integration.

Fig. 1 .
Fig. 1.Illustration of the difference between local and global DVC methods.

Fig. 3 .
Fig. 3. Illustration of the positive and negative effects of regularization (smoothing of the strain field).

Fig. 4 .
Fig. 4. Cross sections of the CT images used in this work.

Fig. 6 .
Fig.6.Illustration of the current and reference position vectors, x and X , and their difference, the displacement, u .

Fig. 7 .
Fig. 7. Illustration of the comparison of gray scale values.

Fig. 9 .
Fig. 9. One-dimensional illustration of the importance of centering image data around zero to avoid positive correlation of uncorrelated sets.

Fig. 10 .
Fig. 10.One-dimensional illustration of the effect of subdividing the volumes into several subvolumes (the same total number of subvolumes can be obtained for several subvolume sizes).

Fig. 11 .
Fig. 11.Illustration of the recursive relation of the current coordinates throughout all time steps.

Fig. 12 .
Fig. 12. Illustration of the definition of the performance measures accuracy and precision.

Fig. 14 .
Fig. 14.The accuracy and precision of the calculated displacement fields in [voxels] obtained by Avizo and GDVC-UU for a virtual Gaussian displacement field applied to the Tozzi-2017 data set.

Fig. 15 .
Fig. 15.Contour plot of the x -component of the displacement fields u ( x, y ) in [voxels] at z = 70 voxels obtained by Avizo and GDVC-UU for the virtual Gaussian displacement field and the Tozzi-2017 data set.The GDVC-UU with a subvolume size of 32 and penalty factor p = 10 .0 and Avizo with subvolume size 16 are compared since they have the same node spacing (lines A and B are for later use, see Fig.17).
Fig. 15.Contour plot of the x -component of the displacement fields u ( x, y ) in [voxels] at z = 70 voxels obtained by Avizo and GDVC-UU for the virtual Gaussian displacement field and the Tozzi-2017 data set.The GDVC-UU with a subvolume size of 32 and penalty factor p = 10 .0 and Avizo with subvolume size 16 are compared since they have the same node spacing (lines A and B are for later use, see Fig.17).
Fig. 15.Contour plot of the x -component of the displacement fields u ( x, y ) in [voxels] at z = 70 voxels obtained by Avizo and GDVC-UU for the virtual Gaussian displacement field and the Tozzi-2017 data set.The GDVC-UU with a subvolume size of 32 and penalty factor p = 10 .0 and Avizo with subvolume size 16 are compared since they have the same node spacing (lines A and B are for later use, see Fig.17).

Fig. 16 .
Fig.16.Illustration of a good approximation of a displacement field for which the strain is very accurate (correct average) but not precise (large standard deviation).

Fig. 17 .
Fig. 17.A cross section of the resulting displacement fields u ( x, y ) in [voxels] along two lines A and B (see Fig. 15 ) obtained by Avizo and GDVC-UU for a virtual Gaussian displacement field applied to the Tozzi-2017 data set.Figures (a) and (b) compare the different methods, (c) and (d) different subvolume sizes and (e) and (f) different amounts of regularization.

Fig. 21 .
Fig. 21.A cross section of the resulting strain fields εxx ( x ) and εyy ( y ) in [%] along two lines A and B (see Fig. 15 ) obtained by Avizo and GDVC-UU for the UU-2018 double scan.Figures (a) and (b) compare the different methods, (c) and (d) different subvolume sizes and (e) and (f) different amounts of regularization.

Fig. 22 .Fig. 23 .
Fig. 22.The accuracy and precision of the calculated strain fields in [%] obtained by GDVC-UU for the UU-2018 double scan and varying amounts of overlap.

Fig. 24 .
Fig. 24.Results in terms of Von Mises strain as obtained by GDVC-UU for real deformation estimated from the UU-2018 CT images (cropped to 600 × 600 × 300).
necessary partial derivatives can be derived from the corresponding relations presented in Section 2 .The partial derivative of the correlation with respect to the gray scale values of the current volume is obtained as, where ∂ s α /∂ x k is directly related to the spline interpolation of the gray scale function.Finally, the derivative of the displacement with respect to the nodal displacements is obtained from Eq. (4) ,∂u k ∂q i = N ki ( X α ) (A.4)

Table 1
Settings for the applied virtual Gaussian displacement field on the Tozzi-2017 data set.The coordinates used are the voxel coordinates, i.e. { x, y, z} ∈ [ 1

Table 2
The accuracy and precision of the calculated displacement and strain fields obtained by Avizo and GDVC-UU for a virtual Gaussian displacement field applied to the Tozzi-2017 data set.The values in italics are best compared, since they have the same node spacing.

Table 4
Orders of magnitude for the worst-case accuracy and precision of the calculated displacement and strain fields due to interpolation errors obtained by GDVC-UU for subvolume sizes in the range 16-50 voxels.