Discrete volumetric digital image correlation for the investigation of granular type media at microscale : accuracy assessment

. The recent development of e ﬃ cient 3D imaging tools such as X-Rays computed microtomography combined with the extension to volumetric images of Digital Image Correlation (DIC) techniques provide new insights on the analysis of materials and structures. Among many other possible ﬁelds of application, geomaterials are good candidates for such investigations, owing to their relative transparency to X-rays and the presence in many samples of a natural contrast suitable for deformation mapping. However, these materials often deform discontinuously at microscale, for instance in the form of the development of a networks of microcracks. Discontinuity is even the dominant rule in granular-type materials such as sand in which the contribution to overall deformation of the microcontinuous phenomena -elastic strains inside grains- are negligible. To investi-gate deformation at the scale of these discontinuous mechanisms, speciﬁc DIC algorithms are required, which override the assumption of continuity of the transformation at the scale of the correlation windows. The recent so-called Discrete-DIC procedure (Hall et al, 2010) is a possible answer. We recall here its general principles and focus on its potential accuracy, from both theoretical and practical points of view. We show that the position and the rotation of individual grains with an average diameter of 500 µ m can be determined from images recorded with a laboratory microCT scanner, with a 15 µ m voxel size, with an accuracy of the order of 1 µ m and 0,1 degree, respectively.


Introduction
The availability of 3D imaging tools such as computed X-Rays tomography and microtomography, confocal optical microscopy or magnetic resonance imaging, combined with the increased power of nowadays computers allowing to process large quantities of data at reasonable costs, made it recently possible to extend to three dimensions the now almost standard 2D image mapping techniques such as Digital Image Correlation (DIC).This opens new perspectives in the analysis of materials and structures and more and more applications of Volumetric Digital Image Correlation are reported in, e.g., biomechanics or rock mechanics.We refer to [1] for a review.Two or three-dimensional image correlation techniques consist in solving a set of optimization problems, providing an evaluation Φ(α min , •) of the unknown mechanical transformation Φ 0 (•) that links a reference image characterized by its grey levels f (•) to a deformed image characterized by g(•) over a subset D where α is a set of scalar parameters that characterize uniquely the transformation and C is a measure of the difference between the grey level distribution over D of image f (•) and back-convected image g(Φ(α, •)).
In standard implementations, the subsets D are regularly shaped (squares or rectangles in 2D, cubes or parallelepipeds in 3D) and distributed over a region of interest in the reference image, while Φ(α, •) is a continuous local map, i.e. such that the displacement u(X) = Φ(α, X) − X is a continuous function of position X in D. Classical shape functions are truncated Taylor expansions of the actual transformation at various order: rigid translation (zeroth order), homogeneous transformations (affine, or first order), quadratic (second order), etc.Such assumptions are satisfactory for most applications where actual local transformations are indeed continuous.They however may fail when strong local discontinuities are observed.Examples of such situations are cracked media, interfacial glide, plastic slip along crystallographic slip planes and granular media.The presented work aims at modifying standard DIC algorithms to cope with the latter class of materials, of which sand is the main example, others being suspensions or composites made of a very soft matrix reinforced with rigid particles.Such materials are made of many grains that evolve in an essentially rigid manner, since the contribution to overall strain of their elastic deformation is in general negligible.Rigid body motions can be very different between neighbor grains, generating strong displacement discontinuities at their contacts, with either rolling, sliding, separation or contacting.The kinematical analysis of such motions at the scale of individual grains can not be performed with standard DIC techniques, even though they may apply at a larger scale [4], where displacement can again be considered as continuous, as neighbor grains keep close during motion.

Discrete Volumetric Digital Image Correlation
In order to measure local motions in granular media, DIC procedures have been modified in two ways.First, the sets D are no longer regularly shaped and distributed but are linked to each individual rigid grain.When the latter exhibits a sufficient internal grey level contrast, the set D coincides exactly with the grain; if not an additional layer has to be added, in order to use the border, i.e., the shape of the grain, as its local signature in the images.Second, the transformation has to be restricted to rigid body motions, described by a translation vector and a rotation (i.e. 6 scalars in 3D, 3 in 2D).The practical implementation for volumetric images of this "Discrete DIC" technique comprises the following four consecutive steps: (i) The image of the undeformed specimen is segmented in order to identify and label individual grains.This can be performed using a watershed algorithm and commercial or academic image processing packages are available for such a task.Note that "oversegmentation", with individuals grains cut into several subgrains is not really a problem, but grains need to be separated from each others.(ii) A mask is defined for each grain, covering the grain plus, if necessary, a typically three-voxel wide layer around the grain.(iii) Standard DIC procedures making use of windows slightly larger then the grains are used to construct a first, non accurate, evaluation of the translation of each grain at a larger scale.(iv) Starting from these initial estimates, the translation and rotation of each grain are determined by solving optimization problem (1) separately for each grain, making use of the subsets defined in step (ii).
In the present study, step (i) has been performed with the commercial software Visilog.Steps (ii) to (iv) have been implemented in the CMV-3D code [2], using the generic registration algorithms available in the image processing library ITK [3].All volumetric images of a sequential mechanical test can be processed by applying incrementally steps (iii) and (iv), the analysis at loading increment (n) taking advantage of the results of increment (n-1).Note that the same reference image can be used for all increments, avoiding the segmentation of all images of the sequence.The procedure has for instance been successfully applied to a triaxial test on a centimetric sand sample made of about 60000 grains with an average 280µm diameter, imaged with 14µm resolution synchrotron computed X-Rays tomography at ESRF ID15 at Grenoble.The detailed description of the results, giving a new insight on the pre-and post-stress peak localization of deformation in such materials are given elsewhere [4,5].We focus here on the evaluation of the experimental errors, from both a theoretical (section 3) point of view and a practical one (section 4).

Theoretical evaluation of errors of discrete DIC due to image noise
The present analysis is restricted to the influence of noise in the images on the accuracy of the rigid grain registration.This does of course not cover all sources of errors.Additional ones are related to reconstruction artifacts (rings, beam hardening effects, phase contrast superimposed on absorption contrast, etc), possible geometric errors in image generation, geometric instability of the CT setup (e.g.magnification variation due to uncontrolled motion of the sources or detectors in laboratory CT scanner), and errors of the DIC algorithms themselves, such as inappropriate grey level interpolation.The proposed theoretical error evaluation should thus be considered as a lower bound on actual errors.Since noise levels are often high in CT images, this analysis is however useful and provides a first realistic quantification of the error, with the additional advantage that it can take into account the shape of an individual grain to quantify the potential error on the rigid body motion of this particular grain.Clearly, the rotation of an almost spherical grain is much more difficult to determine than that of a cubic-like one.A theoretical error model, which might be useful to qualify the accuracy of local measurements, has thus to take into account the shape of the grains.
The approach is essentially an adaptation to discrete DIC of the model developed by [6] for DIC algorithms base on FEM-like shape functions.It is based on the assumption that the grey levels of the voxels carry a white noise, independent from one voxel to the other, but with, for simplicity, the same standard deviation σ n , identical in reference and deformed image, and that noiseless grey level distributions do not have too strong spatial gradients, so that first order expansions are sufficient to evaluate the grey level variations.The last assumption might be questionable in case of sharp images of granular materials.However, the interface between grains and voids is often spread over several voxels, because of limited spatial resolution of the imaging device, so that it might apply anyway.The analysis is developed for a classical "sum of squared differences" correlation criterion, 2 dX but the results are expected to be very similar with other criteria.One obtains the following relation which links the covariance matrix α i α j of the statistical noise on the components of the parameter vector α to be determined, to the image noise level and the image characteristics: where p is the pixel size, d is the dimension of the image space (d = 3 for 3D imaging), and the matrix M is given by and depends on the gradients ∇g of the grey levels distribution in the noiseless deformed image computed for the exact unknown transformation with parameters α e .In this relation Φ ,α p is the partial derivative of Φ with respect to parameter α p .The symmetric positive matrix M, with dimension n × n can be diagonalized in the n−dimensional vectorial space of variations of the parameters α.Let µ i , i ∈ [1, n] be its eigenvalues and U i the corresponding set of eigenvectors with components U i p , or- thonormal with respect to some appropriate inner product.Rewriting the noises α in this new parametric representation, in the form α = i β i U i , and the matrix M as M pq = i µ i U i p U i q one gets: This equality on n × n tensors leads in terms of components relative to the orthonormal basis U n ⊗ U p to the following result, where δ i j is Kronecker's symbol: The noise on the new variables β i are thus statistically independent and the standard deviation of variable β i , given by √ µ i , is proportional to image noise and inversely proportional to the square root of the corresponding positive eigenvalue of the matrix M. In particular, it blows up when this eigenvalue tends to zero, which means that the parameter β i can not be identified with the available image contrast in D.
The error is thus essentially governed by the properties of matrix M defined by equation (3), which needs to be evaluated for the particular case of a 3D rigid transformation decomposed into a translation and a rotation, for which various parametrizations are possible, among which Euler angles, unitary quaternions or orthogonal matrices.Making use of the latter one, we write where R T is the transpose of R and where X 0 is some arbitrary center and T its translation.In the current implementation of discrete DIC, X 0 is the center of the bounding box of the grain.A variation dR of a rotation matrix is such that dR.R T is a skew-symmetric matrix where the dω i are the components of an infinitesimal rotation vector dω such that dΩ.X = dω× X.The errors α in this specific representation are thus parametrized by two infinitesimal vectors, dT giving the error on the translation and dω.The matrix M takes then the form with: Use has been made of the condition of perfect convection of grey levels by the exact rigid transformation (no shape function mismatch error) g(Φ(R e , T e , X)) = f (X), which induces ∇g(Φ(R e , T e , X)) = R e .∇f (X).Note that the above expressions can be rewritten with respect to translation and rotation vector back-convected into the reference configuration (i.e.considering variables [R e ] T .dT,[R e ] T .dωinstead of (dT, dω).With such a choice, the exact rotation matrix R e simply needs to be replaced by the identity operator in relations (8).These expressions can be computed analytically for some representative ideal grain shapes.Consider first a spherical grain with radius a, inner uniform grey level f 1 and surrounding uniform grey level f 2 , with a linear interfacial gradient with slope | f 2 − f 1 )/b where b << a is the "width" of the grain-void interface.Grey level gradients are then along the radial vector e r and have spherical symmetry.The point X 0 is set to be the center of the grain.The sub-matrix M tt is then clearly diagonal, equal to µ t i with On the other hand the matrices M ωω and M tω do clearly vanish.This means that rotation can not be evaluated, as expected.The noise on the components of the translation vector of the center of the grain can however be evaluated as: Consider for instance CT images with a signal to noise ratio ≈ 20, grains with a radius a = 30 voxels and an interface width of b = 6 voxels.The result would be a standard deviation on grain translation components along the three directions of 0.003 voxel.
Consider now a cubic grain of width 2a and similar grey levels.As a consequence of cubic symmetry, the matrices M tt , M tω and M ωω are again proportional to the identity and take the form M tt = µ t i, M tω = µ tω i and M ωω = µ ω i, with The matrix M is then fully diagonal so that the application of the above general theory is straightforward.The standard deviation on translation component is very similar to that obtained for a sphere (6/π times less for the cubic grains).But the standard deviation of the error on rotation components is now With the above numerical values, one finds σ ω i ≈ 1/3000 radian ≈ 0.02 degrees.These theoretical evaluations will be compared to experimental data on real sand grains, with data consistent with the above numerical assumptions.

Experimental quantification of errors
The experiment that has been run is a simple "zero deformation test" on a sand sample.The material was a Hostun S32 sand, which is a fine grained angular siliceous sand from the south of France.The grain sizes range from 500µm to 1600µm with a mean grain size (D 50 ) of about 900µm.The sand filed a cylindrical rigid container sufficiently transparent to X-rays with an inner diameter of about 13mm and an outer diameter of 16mm.The sand was manually compacted so that grains do not move with respect to each other during the experiment.The sample was imaged with the laboratory CT scanner available at Laboratoire 3S-R at Grenoble.This device is a prototype manufactured by the RX-Solutions company, especially devoted to 3D-imaging of samples during in situ mechanical tests.The X-Ray source is a Hamamatsu 150kV seeled tube ,used at 80 kV for the present experiment, and the detector is a 2520 Varian flat panel with a 127µm pixel size.Geometry was set to get a 15µm voxel size in the reconstructed 3D images.Reconstruction algorithms are those provided by the manufacturer.3D images were produced in 16bits and subsequently converted into 8-bits voxels for DIC analysis.Figure 1 gives a typical cross-section of the sample, together with a grey levels profile.Average grey level in the air is about 3500 while in the grains it was near 10500.Fluctuations of grey level in air provide an evaluation of the noise level, as the attenuation should be constant in air.The standard deviation of grey levels in air was about 300 grey levels.Fluctuations in grains could be related to natural fluctuations of composition.Their levels are however limited, with a standard deviation which is also about 300 grey levels.These fluctuations are thus highly likely to be just noise.The signal to noise ratio of these images can thus be evaluated to (10500 − 3500)/300 ≈ 23.Note that 8 bits images are sufficient to render such a signal.The profile in figure 1 illustrates also the intensity of the grey level gradient at the interface, which is not steep, but spread over a thickness of about 6 voxels.These image characteristics are very close to those assumed in the previous theoretical noise evaluation.Note that these noise level evaluations do not take into account reconstruction artifacts, which were limited but not absent (rings can for instance clearly be seen in figure 1).
Right: typical grey level profile, the interface between grains and void is spread over about 6 voxels.
The sample was first imaged in an upright position.It was then slightly tilted with a little hold put between the sample and its support on the CT rotation stage and imaged again with the same parameters.Although the rotation angle could not be evaluated precisely, these procedure ensures that the rotation axis belongs to the (X, Y) plane of the images.Homologous sub-volumes of 1150 × 1150 × 351 voxels have been extracted from the whole volumes for the purpose of the DIC analysis.
These images have first been analyzed with standard continuous DIC procedures using large correlation windows and regularly spaced reference positions.The average deformation gradient of the whole region of interest has been computed with the finite element code Cast3M by means of the procedures described in [7].Its components are the following: The decomposition of this gradient into rotation and distortion leads to a rotation angle of 2.53 degrees around an axis with unitary direction vector m = 0.726e x + 0.687e y + 0.02288e z which indeed belongs to the (X, Y) plane.The components of the distortion matrix were below 5.10 −4 , which provides an estimate of the error level on this macroscopic evaluation of the transformations gradient.Note that such error levels on the components of F lead to errors on the rotation angle below 0.007 degrees and errors on the direction of the rotation axis below 0.2 degrees.
The reference image has been thresholded and grains segmented by means of the commercial Visilog software and the above described discrete-DIC algorithm was applied.Since the natural con-trast in the grains was too small, segmented grains have been enlarged by 3 voxels in order to take the interface into account for the registration.The so-called zero-mean normalized cross-correlation co- , was used together with a first gradient descent optimization procedure.The iteration was stopped when components of displacement increments dT i were below 2.10 −5 voxels and components of rotation increments dω i below 2.10 −2 radians.Convergence was achieved after 20 to 40 iterations starting from a zero rotation and the translation evaluated from classical DIC procedures.Among the 672 grains which were totally included in both images, 640, i.e. about 95%, could be processed successfully, with correlation coefficient at convergence of the order of 0.02.Other grains were detected by a number of iterations above 50 and a correlation coefficient above 0.08 and were removed from the statistical analysis presented hereafter.The discrete DIC analysis provides for each processed grain g a rotation vector composed of the three last components of the unitary quaternion which described the finite rotation, i.e. the vector n g = sin(θ g /2)m g where θ g is the total (positive) rotation angle and m g is the unit vector pointing towards the direction of the rotation axis.The three components of this rotation vector can be expressed into angles θ g i = 2arcsin(n g i ).The average and standard deviation of these angles are reported in table 1, together with the average and standard deviation of the total angle.The full statistical distributions are given in the right plot in figure 2. The average n of all n g allows to define an average rotation R. The local analysis on each grain g provides in addition the translation vector of its arbitrarly chosen center X g 0 .One can substract from this local translation the translation associated with the average rotation R around some arbitrary reference point X0 .The results should be a uniform translation for which the average and standard deviation of the components can be evaluated.While the average is of poor meaning as it is linked to the arbitrary choice of X0 , the standard deviation provides an evaluation of the standard deviation of the error on the translation components of the rigid body registration of the grains.Results are given in table one, and full distributions are ploted on the left plot of figure 2.
It is observed that the overall transformation is recovered by each successfull individual grain analysis, with an accuracy on rotation angles of the order of 0.1 degree and an error on translation components of the order of 0.1 voxels (i.e.1.5µm).This errors are larger than those estimated from the theoretical model relative to a cubic grain, by a factor 5 for the rotation angles, and about 30 for the translation components.The higher error on rotation angles is not very surprising as first, grain shapes are probably somewhere in between cubes and spheres so that the µ ω coefficient given by equation ( 11) is probably overestimated, and, second, because reconstruction artifacts are not taken into account in the analysis.In addition, other sources or errors, as those already mentionned earlier, may have an effect.
The error levels on translation components might sound more surprising, but it should be recalled that the theoretical evaluation is relative to the translation of a very particular point which is the center of the spherical or cubic grain.Translation errors at other positions are also sensitive to errors on the rotation angles: for instance, at a point 10 voxels away from the center of the grain, the additional error due to an error of 0.1 degree on the rotation angle is of the order of 10 sin(0, 1) ≈ 0, 02 voxels, which is much larger than the 0.003 voxels error evaluated theoretically for the translation of the exact center.In a further error analysis, the rotation centers will thus have to be chosen more precisely.

Conclusions and perspectives
Discrete-DIC is a very promising tool for the investigation of the microdeformation mechanisms in granular type materials.Its feasability has been demonstrated in a earlier work [4].We focused here on the evaluation of its accuracy in terms of rotation and translation vector components.A theoretical model has been proposed, which takes into account white image noise, spatial resolutions of the images and geometry of the grains.It seems to be consistent with experimental data evaluated from the analysis of the rigid body motion of a many-grains sand sample, although experimental errors turn out to be somewhat larger, and are of the order of 0.1 degree for rotation and 0,1 voxel for translation, for grains with an average diameter of about 60 voxels and a signal to noise ratio of the images of about 20.The error model will be further tested by taking into account the real geometry of each grain, and by a better selection of the arbitrary rotation center.

Fig. 2 .
Fig.2.Statistical analysis of rigid body motion evaluation on 640 grains processed separately.Left: distribution of fluctuations of translation vector components of center of grains with respect to average transformation.Right: distribution of components of rotation vector expressed in degrees, and distribution of total rotation angle.

Table 1 .
Statistical analysis of rigid body motion evaluation on 640 grains processed separately.Average values and standard deviations of the rotation vector components and total rotation angle expressed in degrees and standard deviation of the translation residue