Keywords: Projection-based in-situ 4D mechanical testing Tomographie, Mesure de Champs, Corrélation d'images numériques, Identification, Problèmes inverses Projection-based in-situ 4D mechanical testing Tomography, Full field measurements, Digital Image Correlation, Identification, Inverse problems

L'analyse quantitative de volumes 3D obtenus par tomographie permet l'identification et la validation de modèles. La séquence d'analyse consiste en trois problèmes inverses successifs : (i) reconstruction des volumes (ii) mesure cinématique par corrélation d'images volumiques (DVC) et (iii) identification. Les très longs temps d'acquisition nécessaires interdisent de capter des phénomènes rapides. Une méthode de mesures, Projection-based Digital Volume Correlation (P-DVC), raccourcit la séquence précédente en identifiant les quantités clés sur les projections. Cette technique réduit jusqu'à 2 le nombre de radiographies utilisées pour le suivi de l'essai au lieu de 500 à 1000. Cette thèse étend cette approche en réduisant la quantité d'informations acquises, rendant ainsi accessibles des phénomènes de plus en plus rapides et repoussant les limites de la résolution temporelle. Deux axes ont ainsi été développés : -d'une part, l'utilisation de différentes régularisations, spatiales et temporelles des champs 4D

(espace/temps) mesurés généralise la méthode P-DVC (avec volume de référence) à l'exploitation d'une seule radiographie par étape de chargement. L'essai peut désormais être réalisé de façon continue, en quelques minutes au lieu de plusieurs jours; -d'autre part, la mesure du mouvement peut être utilisée pour corriger le volume reconstruit luimême. Cette observation conduit à proposer une nouvelle procédure de co-détermination du volume et de sa cinématique (sans prérequis), ce qui ouvre ainsi de nouvelles perspectives pour l'imagerie des matériaux et médicale où parfois le mouvement ne peut pas être interrompu. Le développement de ces deux axes permet d'envisager de nouvelles façons de réaliser les essais, plus rapides et plus centrés sur l'identification de quantités clés. Ces méthodes sont compatibles avec les récents développements « instrumentaux » de la tomographie rapide en synchrotron ou laboratoire, et permettent de reduire de plusieurs ordres de grandeurs les temps d'acquisition et les doses de rayonnement.

Mes trois années de thèse ont été formidables, en particulier grâce à la qualité des échanges avec de nombreuses personnes que je veux ici remercier.

Avant tout, je voudrais remercier l'ensemble des membres de mon jury pour m'avoir fait l'honneur de participer à ma soutenance de thèse. Mes remerciements vont particulièrement à Cino Viggianni, président du jury ainsi qu'aux rapporteurs Jean-Charles Passieux et Luc Salvo. Ils ont passé une partie de leurs vacances à rapporter cette thèse en un temps record, apportant ainsi corrections et pistes de réflexion qui m'ont permis d'élaborer la version définitive de ce document.

Naturellement, je remercie chaleureusement mon directeur de thèse et mentor Stéphane Roux pour sa disponibilité et son encadrement pendant ces trois années. Ses remarques, témoignant de connaissances dans des domaines scientifiques très variés, ont toujours réussi à m'impressionner par leur pertinence et leur originalité (comme, par exemple, la recette du plâtre au citron). De par son écoute et son ouverture, il m'a permis d'explorer tous les possibles, rendant la thèse encore plus passionnante. Je tiens également à remercier Martin Poncelet, qui, depuis le projet Pyréolophore de Licence, m'a transmis sa passion pour la recherche grâce à nos discussions sur des domaines de sciences variés. Je souhaite exprimer ma gratitude envers Thibault Taillandier-Thomas qui a réalisé la thèse précédant la mienne, en partie sur le même sujet. Je le remercie pour ses explications et pour m'avoir transmis ses travaux, me permettant de commencer ma thèse dans les meilleures conditions.

On ne saurait bien sûr oublier le laboratoire. Je souhaite remercier les membres gravitant autour du tomographe, Amine Bouterf et Benjamin Smaniotto, pour leur disponibilité et leur aide lors des différents essais, même pour les plus surprenants (comme les essais de soufflette sur des plots de gélatine, refroidis au ventilateur/glaçon. . . ). Merci aux membres du centre d'essai, Boubou, Remy, pour leur disponibilité pour la conception des essais. Et enfin merci aux personnels du centre de calcul pour leur aide permanente et qui réussissent à faire en sorte que tout le parc informatique et le réseau fonctionnent de façon optimale. Un énorme MERCI à tous les membres du LMT, thésards actuels et anciens, postdocs et personnels qui font que le LMT est tel qu'il est, un lieu (autour du bar notamment) chaleureux, un lieu d'exigence scientifique, un lieu d'entraide et de travail convivial. Je remercie aussi Faema sans qui cette thèse n'aurait pas eu la même saveur.

Not to mention a special thanks to the Solid Division of Lund University for their warm welcome during the 3 months exchange. I would especially thank Steve Hall for his kindness, enthusiasm and for his invitation in Sweden. This was an very nice experience.

Pour elle,

Chapter 1

Processing data-flow in tomography

1.

Tomography

Tomography is an imaging technique that allows the visualization, inside of a specimen, of the material microstructure as revealed by its local coefficient of absorption, in a nondestructive way. Initially developed in the '70s for medical imaging as it enables distinguishing the different tissues, organs and bones, tomography is now widely used in many other fields (e.g., medical imaging for diagnosis [1], biology, material science [2,3,4], etc.) and performed with many different waves (e.g., X-ray, neutron, electron [5], terahertz, optics [6], ultrasound, etc.) depending on the experiment and material absorption and/or scattering.

A large variety of imaging devices exists. For material sciences, tomography has been first developed using bright synchrotron sources. Nowadays, X-ray Computed Tomography scan (CT-scan) is a commonly available and accessible equipment in laboratories. The development of in-situ testing machine (for mechanical, thermal, environmental purposes) has also opened up new avenues in material sciences allowing the specimen to be imaged during the experiment.

The resolution revolution

Many hardware developments of tomography [7] (e.g., use of sensitive CMOS detector, new optics, very intense and monochromatic sources...) have pushed the limits of temporal and spatial resolution, in synchrotron as well as lab-sources. Figure 1.1, extracted from [3] shows the recent evolution of the tomography space and time resolution. The blue area in the figure gives an order of magnitude of the LMT lab-CT characteristics. Figure 1.1: Evolution of the spatial and temporal resolution for X-ray imaging, extracted from [3]. Red dashed curves denote white beam and dashed black curves correspond to the monochromatic, both at synchrotron facilities. The red dotted curve refers to lab-scale tomographs. The blue area corresponds to the LMT lab-CT resolutions.

In standard X-ray tomography, a basic assumption is that the sample should remain static during the acquisition to enable a correct reconstruction. Decreasing the temporal resolution is thus a crucial objective for the analysis of fast phenomenon. A slight motion or absorption change of the sample during the scanning process and the reconstructed volume becomes blurry and hard to exploit.

For many in-situ and in-vivo imaging applications, the temporal resolution is a limit. Although the acquisition rate using synchrotron sources, especially using X-ray, can be very high (e.g., [8] and [9,10] where full 3D volumes were acquired at 20 Hz with 270,000 frames per second), it remains quite long in lab CT-scan (few minutes up to few hours [11]). Moreover, many experiments require multiple full acquisitions of the same sample, at different times, states, loading steps, in order to follow some evolutions. Such multiple scan acquisitions increase the scanning time (e.g., [12,13]) and are thus highly subjected to artifacts: drifts of the beam, thermal changes, uncontrolled motions, etc.

In addition to the development of the experimental setups, software developments operating as a post-processing on the acquired projections have also favored the reduction of the number of required input data. With different regularizations on the reconstruction algorithms (e.g., exploiting a small number of phases, texture sharpness etc.), the number of required projections can be reduced. For example, with binary materials, the number of required projections to reconstruct a 512×512 pixels image was reduced to 7 in [14] and 6-8 in [15]. The use of dictionary / learning also allows reducing the data input.

Processing data-flow in tomography

Initially used predominantly as a means of acquiring 3D images from which diagnoses could be made based on visual judgment, there has been a steady progression towards extracting key parameters from these images, through quantitative analysis. The experiment is indeed often driven by a precise experimental goal (i.e., specific scientific questions). Those goals can be widely different depending on the experiment, on the scientific community, on the material etc.For example, the identification of material parameters for a given model, the position of cracks, the evolution of voids or the front of a fluid invasion can be such scientific questions. Being able to anticipate the goal allows the experiment to be designed. The acquisitions have thus to be quantitatively analyzed in order to extract the sought answers.

The standard sequence to read tomography data, as performed for example in [16,17,18,19] and presented in figure 1.2 aims to extract the quantities of interest. This sequence first consists of the reconstructions of the stacks of radiographs into 3D volumes for the different states. Those volumes are then analyzed to extract measurements (e.g., displacement fields, thermal fields, gray level variations, etc.). Finally, measurements are reduced to parameters and quantities of interest, generally from a model. During this sequence, three inverse problems are performed successively: reconstruction / kinematic measurement / identification. For a correct and homogeneous analysis, it is required to consider how the acquisition uncertainty circulates through the different elements of the chain. Indeed the acquisition noise on the projections becomes reconstruction noise on the volumes that affects the measurements and identification of parameters. The amount of data on each of the previous element of the sequence can be highlighted. To specify orders of magnitude of an experiment composed of 10 scans:

• Radiographs of an acquisition are generally composed of a stack of 1000 approximately of size 1 MPix (12 MPix for the LMT lab-CT in full resolution). The total number of data is hence 10 10 values for the 10 scans.

• Volumes are composed generally of the same amount of data (depending on the projection sampling) say 10 10 unknowns.

• Measurements exploit redundancy to increase the signal to noise ratio (very low for each individual pixel in imaging compared to other sensors). Described on a mesh, the average of [20,21,22,23,18] corresponds of approximately 10 5 degrees of freedom that have to be identified.

• Parameters that drive the model, as described above are composed of very few unknowns, say 10 unknowns.

The blue background of figure 1.2 shows the evolution of the degrees of freedom for the sequence. It is interesting to ask the question of the relevance and necessity of the massive amount of input quantity compared with the number of output parameters.

To begin answering this question, this first chapter is a state of the art organized as follow: a first focus on the data acquisition (the image acquisition of radiographs to projections and the in-situ testing machines that enable capturing sensitive data) and on the three inverse problems (reconstruction -measurement -identification) are described in the following sections. Then, the introduction of projection-based measurements (section 1.6) allows one to short-circuit the sequence and provides a more direct determination of the kinematics, reducing drastically the acquisition data (time and doses).

Radiographs 1.2.1 Imaging technique

Projections are recorded from the relative beam intensity attenuation for each detector position r = [r, z] (where z is parallel to the specimen rotation axis, and r is perpendicular to it) and rotation angle. The Beer-Lambert law links the line integral of the material absorption along the X-ray path L(r) from source to detector at position r, to the recorded intensity I(r, t) at time t (and rotation rotation angle θ(t)): I(r, t) = I 0 (r, t) exp x∈L(r) µ(x)dx (1.1) with µ(x) the linear attenuation coefficient (this absorption coefficient associated to a location in the sample is generally steady in time. An fluid invasion case where the scanned material evolves in time will be studied in the following) and I 0 (r, t), the intensity recorded without the sample on the beamline called flatfield (or 'white field'). To obtain the projections, one has to normalize the recorded intensity with the flatfields (after possibly subtracting off darkfields acquired in the absence of beam) p(r, θ(t)) = -log[I(r, t)/I 0 (r, t)].

(1.

2)

The collection of N θ projections p(r, θ(t)) for each angle is called the sinogram. A flatfield, a radiograph and a horizontal slice in the raw non-corrected sinogram of a rectangular cast iron sample imaged with a parallel beam is shown figure 1.3. It can be seen that the flatfield is composed of a vertical gradient and white patterns. The slice of the sinogram shows vertical bright lines due to those white patterns. Because they are invariant in the t = θ direction, they may be related to the detector or optics (dust and scratches). If non-corrected, those patterns will generate ring artifacts in the reconstruction. The difficulty is that I and I 0 cannot be acquired simultaneously. A standard assumption is that I 0 does not vary in time. A series of flatfields are hence generally acquired before or after the experiment and averaged (otherwise the sample and testing machine have to be moved outside of the beam for the acquisition of flatfields and repositioned). It can be noticed that this is only an approximation [25] hence the flatfields I 0 should be estimated at each instant t, although they are not available 1 .

In addition to the flatfield evolution, different artifacts can affect the radiographs [26,27]. Those artifacts will pollute all the elements of the sequence thus have to be corrected if possible or at least taken into account. A non-exhaustive list of artifacts encontered in the developed applications and possible correction is • Beam hardening. This artifact can be identified from a non uniform texture in the a phase depending on the edge distance (a darker halo in the bulk of the sample). It corresponds to non-linearity in the beam absorption (often related to the spectral width of the source). If all phases have the same absorption function, a non-linear projection model allows the reconstruction to be corrected. Many different algorithms exist [28,29,30]. A simple polynomial correction of the intensity will be used in the following (as proposed in [31]).

• Acquisition noise are assumed in the following developments as white and Gaussian (after the correction of the beam hardening). During the acquisition procedure, an average of radiographs acquired at the same state increases the signal to noise ratio.

• Local patterns such as scratches and irregularities on the optics (see figure 1.3). Its signatures can evolve and move with the beam drift [32] or intensity variations [25] and generate for example ring artifacts after the reconstruction process 2 .

• Phase contrast. This artifact affects the edges of the reconstruction. It is associated with the microstructure and its evolution, therefore it should not highly affect the developed quantitative procedures that aims to track the microstructure motion.

It is noteworthy that many other artifacts can be enumerated (e.g., dead pixels in the detector, metal artifacts, etc.) but are not described because they were not encountered in the following developments. Motion artifact appearing in after the assumed static reconstruction procedure due to motion or an uncalibrated tomograph will be discussed in Chapter 7.

In-situ tests

The first generation of material and mechanical tests were performed ex-situ. The sample is loaded outside of the tomograph, unloaded and positioned to be scanned. This method enables measurement of the remaining changes (cracks, plastic deformations, phase changes [34,35] etc.). For example, cracks in a carbon composite were observed in [36] with the identification of delamination. However, because it is necessary to unload the sample, it is not possible to measure, for example, elastic deformations or opened cracks. Moreover moving the sample at each loading iteration takes time and often introduces miss-positioning and bias. The development of in-situ testing machines has allowed one to visualize inaccessible phenomena. Developed for synchrotron beamlines [37], in-situ testing machines are also available in lab-synchrotron where the space is generally smaller [38]. Testing machines have been made for many different fields: crack propagation [23,39], thermal or environment variation [40], fluid invasion and pressure driven flow [41], granular materials [42,43], biomechanics [44] etc.

One of the difficulty in designing an in-situ testing machine is to allow the sample to rotate without seeing the frame of the testing machine. Different solutions have been chosen, with transparent frames (for tension compression [45,23], 3 point bending [46], etc.) or with a rotating top and bottom grips. In this latter case, an other difficulty is the correct alignment of the grips in order to not introduce unexpected stress and motion during the rotation.

Supported by an Equipex MATMECA project 3 the LMT has bought a new generation testing machine able to perform tension -compression -torsion of large samples (10 cm height loaded up to 20 kN/0.1 kN.m). With different grips, this machine can be used for many different materials and tests. The stiff frame is fixed and the loaded sample rotates thanks to the top and bottom motors. The development of calibration and alignment tools are currently interesting on-going projects. 2 A correction procedure has been proposed in this thesis [33] and enables to separate the superposition of moving scratches and the fixed low frequency background.

3 'Investissements d'avenir' Program under the reference ANR-10-EQPX-37 MATMECA

From Projections to Volumes

Reconstructing the 3D microstructure from the acquired sinogramm (i.e., an inverse Radon transform) is a huge but very simple inverse problem [47].

First the projection operator has to be defined (i.e., 9 parameters for a conic beam [48,49]). This operator is generally given by the tomograph after a calibration procedure. This procedure can be offline (i.e., with a dedicated calibration system made of known patterns [50]) or online (i.e., during the scanning process; based on the sample itself [51,52] or on a known frame [53]). At the LMT lab-CT, the offline calibration is performed with a series of steel balls rotating before or at the end of the experiment (as in [54,55]). The parameters are generally identified only once, even for a multi-acquisition test in order not to move the sample during repeated scans.

In the following the projections and reconstructions are performed using the ASTRA toolbox [56]. Three main types of beam exist: (i) parallel beams in synchrotron facilities (associated with 180 • scanning rotation), (ii) conic beam in lab-sources and (iii) fan beam in medical scanners. The projection operator Π θ defined, the problem can be written in its discrete form with a simple linear system

p = Π θ f (1.3)
with p the sinogramm composed of N θ acquisitions and f the volume. Despite the fact that this system is simple, it can not be easily inverted (the matrix Π θ is of size 10 9 × 10 9 ). Different reconstruction algorithms have been proposed in the literature and can be separated in two categories: Fourier space and Algebraic methods.

Fourier space reconstruction

A common reconstruction method is the Filtered Back-Projection (FPB) [47]. With this method, each projection, p(r, θ) is first "filtered" with a ramp, or Ram-Lak filter, eventually windowed. Ignoring such windowing, in Fourier space, F[p(r, θ)](k, θ) is multiplied by |k|. This ramp is used to normalize the frequency sampling, increasing the high frequencies. Then the filtered "image" is inverse Fourier transformed, and then back-projected, with B θ the backprojection operator, in real space, thereby producing a field g θ (x) that is invariant along the direction e θ . These fields g θ (x) are simply summed over all visited angles θ, producing the sought initial image, f (x)

f (x) = N θ θ=1 B θ F -1 [F [p(r, θ)] (k, θ) • |k|] .
(1.4)

FPB algorithms have been also developed for conic beam [57] (e.g., FDK algorithms). Fourier space algorithms are efficient and fast but it is not easy to incorporate regularizations and a-priori knowledge (e.g., on the microstructure, number of phases, positivity etc.) that would allow to under-sample the acquisition.

Algebraic methods

Another type of reconstruction method consists in decoupling the inverse problem and solving it iteratively [58]. Those methods, referred to as Algebraic Reconstruction Methods (ART), are generally costly but provide a better quality reconstruction than FBP. They are based on the minimization of the functional, Γ ART [f ], equal to the quadratic norm of the difference between the projected reconstructed volume and the acquired projections

Γ ART [f ] = r,θ Π θ [f (x)] -p(r, θ) 2 (1.5) then f = Argmin f † Γ ART [f † ] (1.6)
Additional prior information may easily be added to this functional through regularization, in order to compensate limited angle range for projections, or coarse sampling for example. A standard SART [59] algorithm (Simultaneous ART) is reviewed. In this technique, each projection is considered successively (in a random order as given by a permutation of random angles π). The reconstruction sequence is in 3 steps: (1) projection of the reconstructed volume, (2) comparison with the recorded projection at this angle and (3) back-projection of the difference with the back-projection operator B θ (i.e., the transpose of the projection operator). A convergence criterion on the functional value can be used to stop the number of iterations (Γ ART [f ] < ), with , a threshold value with respect to noise and artifact acquisition. Generally few iterations (N ART ) are required for convergence (it can be noted that for a faster convergence, the volume initialization can be performed using an FBP algorithm).

The algorithm for this method is detailed in Algorithm 1.

Algorithm 1 Standard algebraic procedure, ART (p)

Initialization n ← 1 Initialization f (n) ← 0 Choose a permutation, π, over N θ indices while Π θ [f (n) (x)] -p(r, θ) > do for k ← 1 to N θ do m ← π(k) ρ(r, θ m ) ← p(r, θ m ) -Π θm [f (n) (x)] ∆f (n+1) (x) ← B θm [ρ(r, θ m )] f (n+1) (x) ← f (n) (x) + ∆f (n+1) (x)
Implement additional constraints on f (n+1) (e.g., positivity) n ← n + 1 end for end while Within this algorithm, different regularizations can be included. Those additional constraints enable the reduction of needed projections and the reconstruction to be more physical or highlight specific features. This may come from prior knowledge on the different phases of the sample (as DART algorithms proposed by [15], reconstruction with binary images [14], Total Variation [60]), dictionary learning [61], etc.

The comparison of different reconstruction approaches and regularizations is presented figure 1.4 and extracted from [9]. It is to be noted that this example is a specific test case (2 phases and very few projections) and the performance of the methods are not the same for other cases.

From Volumes to Kinematic Fields

From the acquired volumes at different states, quantitative analysis can be performed. Between different non-rigid image registration methods, Digital Volume Correlation (DVC) is a full field measurement technique developed in the two past decades [20] that aims to capture the 3D displacement field between two (or more) volumes.

DVC has been highly developed in the literature for a large range of experiments and materials; e.g., in bio-mechanics [20,62], for textured or poorly textured metals [START_REF] Buljac | Evaluation of measurement uncertainties of digital volume correlation applied to laminography data[END_REF], woods [START_REF] Tran | 3d mechanical analysis of low-density wood-based fiberboards by x-ray microcomputed tomography and digital volume correlation[END_REF], granular media [START_REF] Forsberg | 3d deformation and strain analysis in compacted sugar using x-ray microtomography and digital volume correlation[END_REF], polymers [START_REF] Hu | Internal deformation measurement of polymer bonded sugar in compression by digital volume correlation of in-situ tomography[END_REF] etc.

Extended from Digital Image Correlation (DIC) in 2D, DVC algorithms consist in minimizing the weighted difference between the reference image f (x), defined for every voxel x = [x, y, z] corrected by a displacement field u(x) and the deformed volumes g(x) (although it is conventional to correct the deformed image by the displacement field (Lagrangian) so that it matches with the reference one, the reference image will be corrected by the Eulerian displacement to be homogeneous with the manuscript notations). The minimization has to be carefully chosen with respect to the acquisition noise [START_REF] Buljac | Digital volume correlation: Review of progress and challenges[END_REF]. With the assumption of a white and Gaussian noise, the resulting functional is

Γ DVC (u) = 1 2γ 2 f N Ω x∈Ω f (x + u(x)) -g(x) 2 (1.7)
with γ f the standard deviation of the noise and N Ω the number of voxels in the region of interest Ω. The solution of the DVC is given by the displacement field that minimizes the functional

u(x) = Argmin v Γ DVC (v) (1.8)
Digital Volume Correlation is an ill-posed problem. A solution to make it well-posed is to regularize the displacement field: reducing the number of unknowns by introducing a priori knowledge on the expected kinematics. Many different regularization techniques have been developed in the literature. In local DVC, the kinematics is divided in independent subsets driven by few degrees of freedoms (translations, rotations [START_REF] Smith | Digital volume correlation including rotational degrees of freedom during minimization[END_REF], warping [START_REF] Bornert | Mesure tridimensionnelle de champs cinématiques par imagerie volumique pour l'analyse des matériaux et des structures[END_REF] etc.). The kinematics is finally interpolated from the subset motions.

Global DVC has been developed more recently [START_REF] Roux | Three-dimensional image correlation from x-ray computed tomography of solid foam[END_REF]. The displacement field is here expressed on a mesh kinematics. This choice offers continuous displacement fields and constitutes an ideally suited interface with finite element computations if needed for instance for future identification purposes. It is proposed to write the displacement field on the basis of shape functions ψ i (x)

u(x) = Nu i=1 u i ψ i (x)
(1.9)

with u i the amplitude of the shape functions and N u the number of degrees of freedom in the mesh (three times the number of nodes) that is, for an order of magnitude, approximately 10 5 . Different choices of 3D meshes and shape functions have been selected in the literature: cubic C8 elements [START_REF] Roux | Three-dimensional image correlation from x-ray computed tomography of solid foam[END_REF] (8-noded cube with tri-linear interpolation), T4 elements [22] (4-noded tetrahedra with tri-linear interpolation), X-FEM C8 elements [START_REF] Réthoré | Extended threedimensional digital image correlation (x3d-dic)[END_REF] etc.

It is crucial to reduce the number of parameters to describe the displacement fields, but as this number is reduced, the risk of excluding the actual displacement from the trial space increases, and hence it may induce a model error. Being able to quantify the error is of utmost importance. The residual field (i.e., the error field at the end of the procedure, with a correction by a displacement field ũ) allows one to (in)validate the procedure and/or the kinematic model ρ(x, ũ) = f (x + ũ(x)) -g(x).

(1.10)

This field contains what was not captured by the kinematics correction: noise and artifacts, model error (e.g., unexpected discontinuities in the kinematics), convergence problems, gray level variations, etc. Looking at those 3D residual fields allows eventually the operator to modify the kinematics model and enrich the basis. An example of application and result is presented figure 1.5 and extracted from [START_REF] Buljac | Digital volume correlation: Review of progress and challenges[END_REF]. A 4-point bending test on plasterboard samples was scanned in-situ (scanning time of 60 min). A continuous C8 mesh was used to regularize the kinematics. Because of a discrepancy between the chosen kinematic basis and the real displacement field (composed of cracks), errors appear in the gray level residual field. In this volume, in addition to an homogeneous noise, a dozen of cracks are visible. Otherwise, no spatial patterns are visible meaning that the global kinematics was wall captured. The amplitude of the displacement field (approximately 2 vox) is very low and highlights the accuracy of the method.

To solve the problem numerically, a gradient decent algorithm is employed. The functional is linearized at iteration l around the current position u l = u l-1 + δu, with f N Ω x∈Ω ρ(x, u l-1 ) + ∇ f (x)δu(x) 2 (1.11) thus the linear system that has to be solved is

M ij δu j = b i (1.12)
with [M ] the hessian of Γ DVC-lin and b the second member

M ij = x ∇ f (x)ψ i (x) ∇ f (x)ψ j (x) (1.13)
and

b i = x ρ(x, u l-1 ) ∇ f (x)ψ i (x) . (1.14)
In order to avoid updating the gradient of the advected volume, one may use the gradient at theoretical convergence (i.e., ∇g(x)).

The regularized functional is generally quite smooth and convex around the global minimum, allowing gradient-based methods to work efficiently. However, in case of large displacements (large compared to the microstructure correlation length), the convergence may be more difficult. Different methods can be applied such as rigid body motion (RBM) initializations (performed in Fourier space for example) or multi-scale approach [START_REF] Helfen | Synchrotron-radiation computed laminography for high-resolution three-dimensional imaging of flat devices[END_REF][START_REF] Mendoza | Extraction of "topological differences" on woven composite materials[END_REF].

The latter consists in low-pass filtering the texture (e.g., convolution with a Gaussian kernel) in order to depress the high-frequency patterns and thus smooth out the functional. Downsizing (sub-sampling) can also be performed after the filtering to reduce the computation time. Used recursively, this procedure is called a pyramidal algorithm, starting from coarse filtered scales to finer scales, revealed very efficient. The final computation is finally performed without filtering.

Different regularizations or penalizations can be added to the kinematics identification for a better convergence or a more physical fields. An easily implementable penalization method is the Tikhonov regularization [START_REF] Bell | Solutions of ill-posed problems[END_REF] or Levenberg Marquardt algorithms [START_REF] Levenberg | A method for the solution of certain non-linear problems in least squares[END_REF] although the amplitude factor may be difficult to tune. A more physical way to regularize the displacement is to introduce a second functional based on the local elastic equilibrium [START_REF] Leclerc | Voxel-scale digital volume correlation[END_REF]. The weight of this functional (i.e., corresponding to a regularization length) allows selecting an internal elastic length. Developed in material science, such regularizations are also used in the image processing field and registration in medical imaging [START_REF] Bouter | A novel model-based evolutionary algorithm for multi-objective deformable image registration with content mismatch and large deformations: benchmarking efficiency and quality[END_REF].

The above procedure can be performed in a 4D (i.e., 3D space and time) framework [23,18,62,[START_REF] Buljac | Digital volume correlation: Review of progress and challenges[END_REF]. The displacement field is hence written u(x, t) = i u i (t)ψ i (x) and the summation in the functional is in space x and time t. When studying space and time fields, the time evolution can be regularized as it is the case for the space field. The chosen time regularization can be weak (e.g., penalization on the time derivative to impose smoothness) or strong (e.g., from a model, described in the next section).

Being able to quantify noise is essential when it comes to challenge the measured results. The impact of the noise in the DVC procedure can be measured. A standard procedure consists in acquiring few volumes in the same reference state (or after having applied a known rigid body motion) and performing a DVC procedure between them. This also measures the ability of the DVC code to register images (inter-voxel interpolation, evaluation of gradients, ... )

From Fields to Parameters 1.5.1 Weighted identification

For the final step of the sequence, N a key parameters called a are to be extracted from the measured fields. The model containing the parameters has to be adapted so that the numerically generated field matches the measured one. Finite Element Updating Method (FEMU) [START_REF] Kavanagh | Finite element applications in the characterization of elastic solids[END_REF][START_REF] Kavanagh | Extension of classical experimental techniques for characterizing composite-material behavior[END_REF][START_REF] Collins | Statistical identification of structures[END_REF][START_REF] Gras | Identification of a set of macroscopic elastic parameters in a 3d woven composite: Uncertainty analysis and regularization[END_REF] is an identification procedure based on the minimization of the weighted quadratic difference between the experimental measured displacement field u(x, t) and the computed field from the mechanical model u c (x, t, a). The functional is performed for every voxels and every time steps.

Γ FEMU (a) = 1 N Ω x,t u(x, t) -u c (x, t, a) 2 Cov -1 . (1.15)
The • Cov -1 denotes the inverse Covariance norm, weighted by the measurement uncertainty.

When the uncertainty matrix is defined for each degree of freedom of a mesh (it is for global DIC the correlation matrix obtained in the previous section [START_REF] Mathieu | Estimation of elastoplastic parameters via weighted femu and integrated-dic[END_REF][START_REF] Neggers | Big data in experimental mechanics and model order reduction: Today's challenges and tomorrow's opportunities[END_REF]), the minimization is performed for each displacement amplitude u i . Each measured nodal displacement is hence weighted by the gray level gradient of the associated elements volume; a node in a non-texture zone will not be considered in the identification as it is associated with high uncertainty. Thus the functional is written:

Γ FEMU (a) = 1 2γ r N n i,t u i (t) -u c i (t, a) 2 M . (1.16)
The convergence is performed with a simple gradient descent algorithm. Hence, with H the Hessian of Γ FEMU (a) and h the second member

H ij a j = h i (1.17)
with S the displacement sensitivity field such that

S i = ∂u c (t, a) ∂a i (1.18) then H ij = S i M S j (1.19) and h j = S j M (u(t) -u c (t, a)). (1.20)
An application of the procedure is presented figure 1.6, extracted from [START_REF] Pottier | Identification paramétrique par recalage de modèles éléments finis couplée à des mesures de champs cinématiques et thermiques[END_REF][START_REF] Pottier | Contribution of heterogeneous strain field measurements and boundary conditions modelling in inverse identification of material parameters[END_REF]. The experimental (top line) and numerical (middle line) are compared showing the residual fields (bottom line). In identification procedures, each available data is useful and may be (suitably) added as constrained to the primal as it helps converging to an accurate optimal solution. First, one may cumulate many time steps to enhance redundancy and increase signal to noise ratio (a typical 4D tomography experiment is composed of 5 to 15 scans). Then, every other modalities can be added in the framework provided that they are weighted with their own uncertainty. With this condition, different functional can be added easily to the FEMU one. Because it is a standard available measurement in mechanical testing machine, the force can be added in the functional [START_REF] Mathieu | Estimation of elastoplastic parameters via weighted femu and integrated-dic[END_REF][START_REF] Bertin | Optimization of a cruciform specimen geometry for the identification of constitutive parameters based upon full-field measurements[END_REF] (called FEMU-F, then FEMU-UF when coupled with the displacement-based functional). With N f the number of force measurements, C f the correlation matrix of the measured loads and F m the measured and F c the computed force

Γ FEMU-F (a) = 1 N f t F m (t) -F c (t, a) 2 C f . (1.21)
Many other modalities can be included in the identification process (strain gauges, camera / X-ray, thermal, X-ray / Neutron DVC that could be performed in a single step [START_REF] Tudisco | Full-field measurements of strain localisation in sandstone by neutron tomography and 3d-volumetric digital image correlation[END_REF]).

Because the model is not always differentiable analytically, it is a common practice to evaluate the finite difference

S i ≈ u c (t, a + δa i .E i ) -u c (t, a) δa i (1.22)
with E i a vector of N a values with 1 at index i and 0 elsewhere. The computation of the sensitivity fields can be very costly. Other optimization procedures have been developed to accelerate the convergence.

In most FEMU procedures, the numerical model is driven by the boundary conditions, extracted from the full field measurement or force sensors. The minimization is hence performed in the bulk of the images as if the boundaries did not evolve. However, boundary conditions can also be considered as parameters to be identified. This solution is attractive as boundary condition uncertainties are not often taken into account and are generally assumed perfect. The procedure hence is able to update both model parameters and boundary conditions.

Integrated methods

A first short-circuit of the above general scheme of data flow is to perform an integrated identification [START_REF] Mathieu | Estimation of elastoplastic parameters via weighted femu and integrated-dic[END_REF]. Instead of performing successively the two last inverse problems, kinematic measurement and identification, the model can be used to warp the gray level volumes and estimate the parameters values directly from the images. If the full field measurement has well converged, the integrated solution should be rigorously equivalent to the weighted FEMU measurement. It is noteworthy that this technique is natural as the frontier between measurement and model is very fuzzy (should the amplitudes of a polynomial fit be considered as measurement or identification?).

The functional of the integrated identification Γ Int is similar to the DVC functional (i.e., based on images)

Γ Integrated (a) = 1 2γ 2 f N Ω N t x,t f (x + u c (x, t, a)) -g(x, t) 2 . (1.23)
However the sought displacement field is not defined as the most general a mesh kinematics but from the model which can exploit the very same mesh. At iteration l: u(x, t, a l ) = u(x, t, a l-1 ) + S l (x, t) δa (1.24) with S l (x, t) the sensitivity field defined for each voxel space and time. Integrated approaches have been developed for crack identification using the analytic and differentiable Williams' series [START_REF] Mathieu | Analyse de la tenue mécanique d'un liner en titane: apport des mesures de champs cinématiques[END_REF]. It provides the position of the crack tip and evaluation of stress intensity factors, in 2D and 3D [START_REF] Roux | Digital image correlation and fracture: an advanced technique for estimating stress intensity factors of 2d and 3d cracks[END_REF][START_REF] Mathieu | Analyse de la tenue mécanique d'un liner en titane: apport des mesures de champs cinématiques[END_REF][START_REF] Roux-Langlois | Dic identification and x-fem simulation of fatigue crack growth based on the williams' series[END_REF][START_REF] Réthoré | Automatic crack tip detection and stress intensity factors estimation of curved cracks from digital images[END_REF][START_REF] Lachambre | Extraction of stress intensity factors for 3d small fatigue cracks using digital volume correlation and x-ray tomography[END_REF].

The residual field (in gray level for the integrated approach) allows (in)validating the procedure. Model errors can thus be identified and the model be enriched until the result is judged satisfactory [START_REF] Neggers | Improving full-field identification using progressive model enrichments[END_REF].

1.6 P-DVC short-circuit -from Radiographs to Parameters P-DVC for measurement A recently developed Projection-based Digital Volume Correlation method [START_REF] Leclerc | Projection savings in ct-based digital volume correlation[END_REF] allows one to short-circuit the sequence by identifying the kinematics directly on the projections and "helped" by an initial acquisition. Figure 1.7 shows the projection-based measurement short-circuit. First, a reference volume is acquired once, before, during or after the experiment. This volume has then to be positioned in the 3D space such that its projections match the recorded deformed projections. This procedure is not performed through the minimization of 3D residual fields but a stack of projected residuals at different angles. With Π θ the projection operator [56] at angle θ, p(r, θ) the recorded deformed projections, N θ the P-DVC functional can be written

Γ P-DVC (u) = 1 2γ 2 f N r N θ r,θ Π θ [f (x + u(x))] -p(r, θ) 2 (1.25)
The linearized functional around the solution, at iteration l, with f (x) = f (x + u l-1 (x))

Γ P-DVC-lin (δu) = 1 2γ 2 f N r N θ r,θ Π θ [ f (x)] -p(r, θ) + Π θ [∇ f (x)δu] 2 (1.26) CHAPTER 1. PROCESSING DATA-FLOW IN TOMOGRAPHY
The quantity ρ P-DVC (r, θ, u) = Π θ [ f (x)] -p(r, θ) defines the projected residual fields for each angle, quantity to be minimized. Exactly like the 2D or 3D residual field defined in DIC or DVC, the projected residuals highlight what was not captured by the 3D correction, here projected in the 2D space detector. The displacement field can be regularized using finite element mesh as in global DVC or using integrated approaches as it will be the case in the following development of this thesis.

P-DVC has been developed, during the PhD of Thibault Taillandier-Thomas [START_REF] Taillandier-Thomas | Développement de la tomographie spatio-temporelle pour le suivi d'essais mécaniques[END_REF][START_REF] Taillandier-Thomas | Soft route to 4d tomography[END_REF][START_REF] Taillandier-Thomas | Measurement of 3d displacement fields from few tomographic projections[END_REF] for the study of a fatigue crack in a cast iron sample, imaged at the ESRF synchrotron at different stages of the test. In [START_REF] Leclerc | Projection savings in ct-based digital volume correlation[END_REF], the comparison of P-DVC with a standard DVC approach has shown that (i) P-DVC provided the same displacement field as DVC and (ii) a large part of the projections could be unexploited without prejudice on the result. With a mesh composed of 303 T4 elements, the displacement field comparison (see figure 1.8) was performed using a P-DVC procedure with all the 600 acquired projections, 48 and finally only 2 orthogonal ones. The displacement results are very close and prove that 600 projections are not necessary, and hence are highly redundant. They can be reduced to only two projections which leads to a potential gain in acquisition time (thus data and doses) of about 300. The two last projections are selected to be orthogonal such that the nodal displacements are sensitive to the 3 directions of the space. Indeed, the projection operator does not offer sensitivity in the X-ray projection axis (e.g., a displacement along the direction of the beam can not be seen). In the P-DVC procedure, because the minimization is performed simultaneously with the two views, the solution can be found. It can be noted that 90 • angle between two projections is not a requirement in the procedure. Even a small angle can be enough to provide sensitivity in all space directions. An initial optimization, prior to the experiment can be performed to select projection angles. This optimization will be discussed in the manuscript.

The projections of the volume and residual fields are shown figure 1.9. The two projections of the sample at 0 and 90 degrees and two residual fields are shown on the first and second line respectively. The top parts of the residual images are composed of very high positive and negative values (signature of motion) that have not been corrected. Those parts were masked in the procedure thus correspond to the initial residual fields. The huge decrease of the residual fields when motion is correctly accounted for can be seen clearly. Alternative methods: It is interesting to note that this "multi-view" P-DVC framework is similar to stereo-correlation methods that aims at measuring 3D surface fields from multi-view cameras. In fact only the projection operator has to be changed to switch from multi-view P-DVC to stereo. Projection operator Π θ 3D reference (e.g., CAD [START_REF] Dufour | Cad-based displacement measurements with stereo-dic[END_REF][START_REF] Pierré | Finite element stereo digital image correlation: framework and mechanical regularization[END_REF]) Reference image f (x) Mesh Mesh Φ i Set of images X-ray projections p(r, θ)

A similar method has been developed by Khalili et al. in 2016 for the measurement of the rigid body motion of few individual grains in a compression test on a granular media from a discrete projected DVC [START_REF] Khalili | Tracking and modelling small motions at grain scale in granular materials under compression by X-Ray microtomography and discrete simulations[END_REF][START_REF] Khalili | Discrete digital projections correlation: a reconstruction-free method to quantify local kinematics in granular media by x-ray tomography[END_REF] called Discrete Digital Projection Correlation method (D-DPC). As in P-DVC, the motion is estimated from 2 orthogonal projections.

In a similar spirit, a method has been developed in the medical imaging namely 3D/2D registration [START_REF] Markelj | A review of 3d/2d registration methods for image-guided interventions[END_REF]. Developed in 2005 [START_REF] Zeng | Respiratory motion estimation from slowly rotating x-ray projections: Theory and simulation[END_REF][START_REF] Zeng | Estimating 3-d respiratory motion from orbiting views by tomographic image registration[END_REF], this technique has been used for example for the calibration of scanners (often using affine transformation [START_REF] Ouadah | Selfcalibration of cone-beam ct geometry using 3d-2d image registration[END_REF]), to track specific moving patterns [START_REF] Otake | Automatic localization of vertebral levels in xray fluoroscopy using 3d-2d registration: a tool to reduce wrong-site surgery[END_REF]. For the measurement of displacement fields (and sometimes the correction of the reconstruction itself [START_REF] Li | Enhanced 4d cone-beam ct with inter-phase motion model[END_REF]), this method has been applied for the measurement of moving thorax or head [START_REF] Delmon | Recalage déformable de projections de scanner X à faisceau conique[END_REF][START_REF] Delmon | In-room breathing motion estimation from limited projection views using a sliding deformation model[END_REF][START_REF] Suzuki | Projection-based motion estimation for cardiac functional analysis with high temporal resolution: a proof-of-concept study with digital phantom experiment[END_REF].

Different methods have been developed to speed up the computation time of Projection based approaches (e.g., GPU in [START_REF] Leclerc | Projection savings in ct-based digital volume correlation[END_REF], Convolutional Neural Network in [START_REF] Miao | Real-time 2d/3d registration via cnn regression[END_REF][START_REF] Toth | 3d/2d model-to-image registration by imitation learning for cardiac procedures[END_REF]).

Structure of the dissertation

This manuscript aims to progressively introduce new developments in projection-based identification allowing to reduce the number of acquired data. As for the space where it is natural to couple the independent degrees of freedom (e.g., with Global procedures), working on the time regularization is a key point as it unifies all time states in a single identification framework. By relating all instants of time together, redundancies are not only observed in space but also in time. The space data reduction performed in [START_REF] Leclerc | Projection savings in ct-based digital volume correlation[END_REF] can be extended in time.

The document is structured in two parts. The first part (Chapter 2-3-4-5-6) extend the P-DVC development using a reference image to identification and single-view ultra fast procedures. The second part (Chapter 7-8-9) is dedicated to dynamic tomography. test on a tensile compression testing machine. A single reference state is acquired on the virgin unloaded sample. The entire procedure takes 2.5 hours for the reference acquisition and 13 loading states that took 35 second each. This chapter is extracted from [START_REF] Jailin | In situµ ct-scan mechanical tests: Fast 4d mechanical identification[END_REF].

• Chapter 3: P-DVC is extended to a single projection per time step. The acquisition of projections is continuous while the sample is loaded and rotates. The displacement field is regularized in space and time such that each loading step contributes to the global identification. The lack of sensitivity at some steps is hence retrieved at other steps. Because of a large number of degrees of freedom, model reduction techniques (PGD) are developed. The solution is identified mode by mode successively. The application is performed with a tensile test on a slender cast iron sample until failure and imaged in-situ in the LMT lab-CT. One reference state is acquired in 20 minutes and the experiment composed of 127 loading steps is carried out in 6 minutes. This chapter is a Part1 of a two part article extracted from [START_REF] Jailin | Fast 4d tensile test monitored via x-ct: Single projection based digital volume correlation dedicated to slender samples[END_REF].

• Chapter 4: P-DVC is extended to integrated approaches. This chapter is related to Chapter 3 and is applied on the same data set. The sought displacement field is described by a plastic beam model. Two different methods are described: FEMU identification using the displacement fields of Chapter 3 and integrated methods. This chapter is the Part2 of a two part article extracted from [START_REF] Jailin | Fast 4d tensile test monitored via x-ct: Elastoplastic identification from radiographs[END_REF].

• Chapter 5: The projection based measurement is here extended to "nonmechanical" identification with gray level variation. The fluid invasion on a porous sandstone sample has been imaged using neutron tomography. The moving 3D front of the fluid is tracked using projections. The fluid front is parametrized in space (low order polynomials), time and saturation profile. Because of numerous unknowns, the 3 coupled problems are solved using PGD methods. The result is a front tracked at 5 Hz time resolution (300 times faster than standard 3D scans). This chapter is extracted from the work performed at Lund University, Sweden [START_REF] Jailin | Fast tracking of fluid invasion using time-resolved neutron tomography[END_REF].

• Chapter 6: P-DVC is extended to ultra fast in-situ vibration measurements. The measurement of the modal basis of a vibrating sample is demonstrated. With a random time sampling, and very few projection angles, the microstructure and a small number of spatial modes of vibrations shared by all projections and hence these common modes can be reconstructed since the microstructure is known. Two numerical examples are presented and show the accuracy of the method. With this method, the only speed limitation for vibrating phenomenon is no longer the acquisition time but the exposure time (that is (or can be) extremely small).

The second part of the manuscript is dedicated to dynamical tomography. The reference volume is here imaged with motion and the identified motion is used to correct the reconstruction procedure (with an update of the calibration or by changing the calibration parameters), as shown figure 1.12.

Figure 1.12: Projection-based DVC without reference image used for dynamic tomography.

• Chapter 7: The online calibration of the tomograph using P-DVC techniques is performed. Because of approximate calibration parameters of the tomograph, the reconstruction of the scanned sample is very blurry. This low quality reconstruction is used in a P-DVC framework with its own projections. The measured displacement field allows calibration to be perform and good and sharp volumes to be reconstructed. The spatial displacements are regularized with rigid body motions and their time evolution is regularized with polynomials and sinusoidal functions. A PGD procedure is implemented to capture the first motion modes. This work is extracted from [START_REF] Jailin | Self-calibration for lab-µct using space-time regularized projection-based dvc and model reduction[END_REF].

• Chapter 8: A dynamic tomography method is developed in 2D with numerical phantoms. From an initially blurry reconstruction, a coarse displacement can be measured. Then, with a motion compensated reconstruction procedure based an a SART algorithm, the volume can be progressively updated until it it reaches a satisfactory quality. The procedure is applied on two numerical examples: a moving Shepp Logan phantom and a pulsating checkerboard. This work is extracted from [START_REF] Jailin | Dynamic tomography of moving samples[END_REF].

• Chapter 9: The dynamic tomography method is developed in 3D with a real test case. The dynamic correction of a respiratory motion based on 10 full CT of a patient lungs is performed in a multi-scale framework. The 10 scans have been acquired at the CREATIS laboratory. The 4D motion identification allows the reconstruction to be much shaper around the diaphragm.

Supplementary materials

Articles or preprints issued from published or submitted work, are sometimes augmented by some additional comments. In such cases, the new text is added in a red box, not to confuse it with the original version. A recently proposed "Projection-based Digital Volume Correlation" (P-DVC) method is here extended to a cone-beam lab-tomograph in which a mechanical test is performed. It consists in a crack propagation test for an elastic-brittle gypsum specimen. The kinematic analysis is done based on a reduced finite element modeling for which the appropriate boundary conditions, and crack propagation stage, are determined from the radiographs. Considering only two projections per loading step, an integrated model-based analysis of the entire test provides a full space and time identification of the kinematics, including the crack position and the determination of two material parameters. This is achieved with a drastic reduction in acquisition time as compared to classical Digital Volume Correlation analysis. In the considered example, acquisition time is cut down by a factor of 350.

Introduction

Identification and validation of mechanical models used to predict the behavior of materials and structures has been and still is the central focus of experimental mechanics. However, the ever increasing sophistication of mechanical models, and the multiplicity of scales required to assess and quantify the microscopic phenomena at play present every day more challenging demands to mechanical tests. During these last decades, this trend is balanced by the accessibility to more and more powerful measuring and imaging techniques, (and in a lesser extent, richer and more accurate loading setups), from which complex loadings, specimen geometry and full-field measurements can benefit.

The development of Computed Tomography (CT) has been a revolution in materials science [1,2,3]. Accessing the intimate micro-structure of solids in a non destructive way has opened new horizons. The recent evolutions of these imaging techniques give access to higher spatial and temporal resolution. First developed on synchrotron facilities with high brightness monochromatic coherent X-ray beam, X-ray CT scanners are now an accessible equipment in laboratories. Moreover, their state-of-the-art performance may in favorable cases compare well with those of large scale synchrotron tomography.

The development of material tests coupled with tomographic images has been studied in the recent years. After ex-situ testing where the materials are deformed outside of the tomograph, the recent evolution of CT makes in-situ tests possible [4,5,6]. In the latter case, material specimens are deformed inside of the tomograph. This method allows for characterizing new mechanisms (i.e. crack opening [7,8], shear banding, fast transformation [9,10]). Recent works developed in synchrotron facilities reported 20 Hz scans for the study of crack propagation [11]. Combined with imaging analysis techniques such as full field measurement, X-ray CT becomes a powerful tool for experimental mechanics.

Digital Volume Correlation (DVC) [12,13] is a full field measurement method that aims at capturing the way a solid has been deformed between two states captured in 3D images. Extended from 2D Digital Image Correlation, DVC permits to measure accurately 3D displacement fields based on the micro-structure of registered volumes. As in standard mechanical tests, an experiment is generally composed of several loading steps. DVC can be performed with all deformed states of the same specimen. This space-time analysis of the displacement field, called 4D-DVC [14,15] permits for example to identify a constitutive law exploiting all loading steps globally.

Nevertheless, the major limitation of CT imaging especially in lab-tomograph is the acquisition time (around one hour).Demanding application for very high resolution can reach more than 20 hours [16,17]. This limit does not allow for the visualization of time-dependent behaviors that may be considered as artifacts, blurring the reconstruction [18]. In [19,7], 0.5 to 1 hour were spent to wait for relaxation or creep behavior at each loading steps. Yet another limitation is the study of biological materials where a low dose of x-rays is to be used to avoid radiation induced damage [20]. One possible procedure to circumvent these difficulties is the Projection-based DVC (P-DVC) [21]. Instead of working with full 3D images for every loading step, P-DVC consists in operating directly on radiographs (i.e., projection of the deformed volume). A very small number on radiographs turns out to be sufficient to measure a displacement after a first complete 3D image of the reference state has been (classically) reconstructed. In [22,23], it has been shown that the measurement of the 3D displacement field of a specimen of cast iron with a fatigue crack was obtained from no more than two radiographs (instead of 600). The exploited experiment had been performed in situ in the European Synchrotron Research Facility (ESRF) in Grenoble. Ideas in a similar spirit for rigid grain tracking are recently proposed in [24].

The present paper aims at extending this P-DVC methodology on lab-tomograph images, addressing the identification of material properties from a complex geometry. In Section 2.2 a brief overview of DVC and P-DVC procedures is given. This image processing technique is regularized with a coupling to a mechanical model. In Section 2.3, the Double Cleavage Drilled Compression (DCDC) test case is introduced. The proposed methodology is applied in Section 2.4 to the test case, and the obtained displacement fields and projection residuals are presented. The identification of model parameters is discussed. Finally, Section 2.6 recapitulates the main results of the paper, and proposes a discussion on possible ways to improve the proposed methodology.

Method

Digital Volume Correlation

Digital Volume Correlation (DVC) [12,13] is a full field measurement technique for the 3D displacement field that relates two 3D images, one for the reference state and one for the deformed one. It consists in the registration of an image f (x) in the reference configuration and a series of 3D images g(x, t) in deformed configurations indexed by time t. The term x denotes the Cartesian coordinates of the reconstruction. The DVC procedure (written here with the Eulerian transformation to unify notations considering the next section) is the minimization of the quadratic difference between the reference image corrected by the measured displacement u(x, t) and the image of the deformed state (called "deformed image" for conciseness)

χ 2 u (t) = 1 2|Ω|γ 2 f x∈Ω (f (x -u(x, t)) -g(x, t)) 2 (2.1)
where γ 2 f is an estimate of the noise variance over the image, and |Ω| the volume (number of voxels) of the region of interest, Ω. A first kinematic regularization of the displacement field can be introduced by global DVC. In global DVC [25], the displacement field is expressed on a reduced basis composed of a set of fields Φ i (x) such that

u(x, t) = i u i (t)Φ i (x).
(2.2)

A general framework for kinematic bases well-suited to mechanical modeling is those used in the framework of the finite element method. This choice offers continuous displacement fields and constitutes an ideally suited interface with finite element computations if needed say for future identification purposes.

The displacement field is finally obtained from the minimization of the functional with respect to the degrees of freedom u i (t) (i.e., nodal displacements)

u i (t) = Argmin v i (χ 2 v (t)). (2.3)

Global P-DVC procedure

In order to reduce acquisition time (especially the acquisition of the entire sinogram for every deformed states), it is proposed to work with few radiographs instead of reconstructed 3D volumes. The proposed method is Projection-based Digital Volume Correlation (P-DVC) [21,22,23]. The first step of the procedure is the acquisition and reconstruction of a reference image, f (x), (generally without loading) from a complete set of radiographs with a fine sampling of all rotation angles. Then the deformed states are characterized through a much reduced set of radiographs, sampled at a few, N θ , selected angles. Because the number N θ of needed projections may be very small, the deformed image cannot be reconstructed solely based on those radiographs. The idea is on the contrary to obtain it from the reference image, by fitting the displacement field so that the projections of the deformed volume match the few available radiographs.

To clarify our notations, it is recalled that reconstruction is based on images that are computed as the cologarithm of intensities of radiographs normalized by a flat field (i.e., radiograph captured without specimen). The latter are denoted as s(r, θ, t), with r the coordinates of the detector, θ the rotation angle and time t. In order to avoid the confusion with the raw radiographs, they are referred to as "sinogram" in the following, irrespective of the number of angles θ.

In the same spirit as DVC, the displacement field is obtained from the minimization of the quadratic difference between the (re-)projected reference 3D image corrected by the displacement field Π k [f (xu(x, t))] and the sinogram in the deformed states at time t, captured at just a few angles s(r, θ k , t)

χ 2 u (t) = 1 N θ |Ξ|γ 2 s k,r (Π k [f (x -u(x, t))] -s(r, θ k , t)) 2 (2.4)
where the double sum over (k, r) is the discrete integration over all pixels r ∈ Ξ at all projection angles θ k . γ 2 s is the variance of the sinogram noise, and |Ξ| the area (number of pixels) of the detector (or its utilized part). At this step, the displacement field may be expressed on the reduced basis as previously defined. Minimization of the functional leads to the displacement field from the nodal values

u i (t) = Argmin u i (χ 2 u (t)). ( 2

.5)

A Newton algorithm may be used to solve this problem iteratively by a progressive correction of the displacement field, correction obtained from the tangent linear problem about the current point. Because the correction displacement field, δu i , is in the range of small perturbations, the integrand in Eq. 2.4 reference image corrected by this displacement field can be expressed as

Π k [f (x -u i (t)Φ i (x))] ≈ Π k [f (x -u(t)) -δu i (t)Φ i (x)∇f (x -u(t))].
(2.6) Finally the linear system that has to be solved is

δu i = N -1 ij n j (2.7)
with N the Hessian matrix and n the second member vector based on the residual field

N ij = θ k ,r Π k [Φ i (x)∇f (x -u(t))]Π k [Φ j (x)∇f (x -u(t))] (2.8
)

n i = θ k ,r (s(r, θ k , t) -Π k [f (x -u(t))])Π k [Φ i (x)∇f (x -u(t))] (2.9)
As in DVC, the residual field at convergence, ρ(r, θ k , t)

≡ (s(r, θ k , t) -Π k [f (x -u(t))])
gives very precious information on the quality of the solution. However, in P-DVC the residual fields (N θ 2D images) are in the projected domain (parameterized by r). They are also affected by acquisition noise and artifacts, reconstruction and projection inaccuracies and mesh discretization errors (on a crack path for example).

Let us emphasize here a very attractive property of this approach. Because computations operate in the projection space, acquisition artefacts can be corrected and noise taken into account. Noise for instance is close to being white, and Gaussian, what legitimates the use of a simple L2 norm of the projection residuals. Some additional artefacts can also be characterized (dead detector pixels, beam hardening, ...), and hence accounting for these features is feasible. It is fair to add that some artifacts can be corrected in a pre-processing stage, leading to enhanced quality of the reconstruction. However, it should also be emphasized that for the present application, more that the reconstructed volume, each individual projection is to be corrected, and this is very demanding. This contrasts considerably with classical DVC, where, even after a pre-processing of the radiographs to erase some artefacts, noise as well as uncorrected or partially corrected artefacts have been processed in the reconstruction and hence spatial correlation have been built, and spatial stationarity is broken. In this case, a complex metric (inverse of the covariance matrix) should be used in DVC in order to compute the theoretically correct cost function. This is so complex than no-one ever considers such a comparison metric in practice, and hence registration is sub-optimal.

Integrated P-DVC

Aiming to reduce the number of degrees of freedom, another mechanical regularization can be introduced. The displacement field can be expressed as a linear combination of fields obtained from a finite element computation based on the same meshed domain and exploiting a constitutive law. The kinematics is then controlled by fewer degrees of freedom, p i , for example those parameterizing the boundary conditions of the test or a specimen geometry including for example a crack front position.

Depending on how the modeling (including the mechanical behavior) can be trusted either a small weight can be given to the regularization, so that it acts as a mechanical low pass filter [26,23], or a large weight can be ascribed enforcing the mechanical behavior over the entire domain, and thereby reducing drastically the number of effective degrees of freedom, and hence the computation time. In the application, only the limit of an infinite weight will be considered, so that a homogeneous elastic behavior is strictly prescribed. It is noteworthy that boundary conditions are understood here at large, and may contain additional parameters such as the position of the crack tip. Moreover, because the displacement fields are directly expressed on mechanically admissible fields, the present formulation allows for the identification of model parameters. In this framework, the chosen unknown parameters, be they boundary conditions, geometry or constitutive parameters, can be considered altogether as the unknown p i with 1 ≤ i ≤ N .

A similar Newton's descent method is chosen as previously. However, as the expression of the new degrees of freedom (or their incremental change) can be expressed from the nodal displacements, the present formulation can be deduced from the previous one. It is needed to compute the sensitivity fields

S ij = ∂u i ({p}) ∂p j (2.10)
The corrections δp i to the current estimate of the parameters p reduces to

δp i = H -1 ij h j (2.11)
with a Hessian matrix H and second member h

H ij = S mi N mn S nj or H = S † N S (2.

12)

and

h j = S ij n i or h = S † .n (2.13)
After convergence, the residual fields are instrumental in evaluating the validity of the mechanical model and the choice of boundary conditions. If deemed necessary, this choice can be revisited in order to better capture the kinematics.

Use of load measurement

The previous procedure only takes into account images. However, purely kinematic information cannot set a stress scale and hence cannot provide access to parameters such as the Young modulus for an elastic law. Therefore, it is useful to also take into account any available measurement data such as that of a load cell. With the previous integrated procedure, the load F c (t; p) can be computed with any chosen set of parameters p and compared to the measured ones F m (t).

An additional term to the functional can be added [27,28] in order to perform the minimization in displacement and force

χ 2 (t) = χ 2 u (t) + χ 2 F (t) (2.14)
with the force functional

χ 2 F (t) = 1 γ 2 F (F c (t) -F m (t)) 2 (2.15)
where γ 2 F is the noise variance of the load cell, (this formula assumes that only one measurement is being used).

Let us stress that the extensive quantities are the cost functions χ 2 times the number of measurements they involve (provided they are statistically independent). Thus the joint minimization of two sets of measurements A and B should be based on

χ 2 A∪B = N A χ 2 A + N B χ 2 B N A + N B (2.16)
where N A (resp. N B ) is the number of measurements of A (resp. B). This ensures that the expectation value of each χ 2 is unity. Thus when combining different loading states, individual χ 2 u and χ 2 F should be weighted by N θ |Ξ| and 1 respectively. When several (N t ) loading states are considered, because each instant involves the same expression, both N A and N B are multiplied by N t , and hence the above weighted means remains invariant, and N t can be ignored.

In the following, two strategies will be followed: in the first one, the measurements of displacement are first performed at each loading step independently, and then from them and the force values, the Young's modulus will be estimated; A second strategy consists in determining forces and displacement jointly with the above functional, that is coupling different kinematic measurement when no damage is expected.

In the following test case, an elastic law will be used. The Young modulus, E, can be easily identified considering all N t loading steps in a spatio-temporal framework [29]. It consists in the minimization of the total functional (summed over time) with respect to the unknown parameter. Because the force is proportional to E, the force can be computed with a unitary Young modulus, F u (t), and hence F c (t) = E.F u (t). Minimizing the total cost function with respect to E thus leads simply to and hence

t (F m (t) -EF u (t))F u (t) = 0 (2.
E = t F u (t)F m (t) t F u (t) 2
(2.18)

Case study 2.3.1 DCDC test

The Double Cleavage Drilled Compression (DCDC) test refers to a parallelepipedic shaped specimen with a square cross section normal to the sample length axis and a centered cylindrical hole drilled in the center of one of its lateral face. The specimen is subjected to a compressive load in the length direction (perpendicular to the hole axis). This type of test permits to study the fracture of brittle materials (originally used for studying crack propagation in glass [30,31]). Under compression, due to the geometry, a tensile stress concentration is generated around the central hole and initiates two symmetric mode I cracks that propagate on the mid-plane of the sample while the load increases [32,33,34].

Test case presentation

The DCDC specimen used for the present study is a 8.40×8.40×21.79 mm 3 plaster sample with 1.4 vol.% of copper powder (approximative size of 150 µm) in order to provide markers well distributed in the bulk with a sharp contrast for X-Ray imaging (the X-ray absorption of copper (Cu) is much higher than that of gypsum (CaSO ○, the X-ray source and 3 ○, the detector. Note that the distance between the source and detector is not the one used in this article

The sample (figure 2.1(a)) is scanned at the LMT lab-tomograph (cone beam, 149.6 kV, 108.5 µA, W target) inside of an in-situ testing machine similar to the one designed by Buffière and Maire [35] (figure 2.1(b)). The voxel size at full resolution is set to 15 µm. The complete scan of the reference state consists of 1500 radiographs captured at equally spaced angles ranging over a 360 • revolution. One flat-field is acquired after conditioning and before the experiment in order to perform a standard flat-field correction. One darkfield is also acquired before the experiment. Because the reference volume is just to be acquired once, in order to enhance its quality and reduce the acquisition noise, each radiograph is an average over 50 projections. (This high number was found a posteriori to be unnecessary because the procedure is not sensitive to uncorrelated noise). Finally the entire reference scan (including storage time) requires 2.5 h to be completed.

Supplementary materials

One of the key quantity in the P-DVC measurement procedure is the projected gray level sensitivity per degree of freedom (for the global procedure:

σ i (r) = Π k [Φ i (x)∇f (x)]
). From these fields are computed the global Hessian matrix: N and the integrated one H as well as the covariance matrix og the measurement when a white Gaussian noise affects the images: C = 2γ 2 s H -1 . The gray level sensitivity is composed of two terms: (i) the regularization choice (i.e., choice of a kinematics model) that has to be carrefully selected before the experiment and (ii) the projected gray level texture. The latter one depends on the scanned materials. Some materials have an very good X-ray intrinsic texture because they are composed of particles (e.g., cast iron, reactive materials, granular materials), fibers (e.g., wood, carbon composites) or because of their shapes (e.g., foams). However DVC or P-DVC on poorly textured materials can be difficult (e.g., iron, aluminum, ...).

In this application, it was chosen to enhance the texture by adding a pinch of copper powder (1.4 vol.%) before molding it. Because the sample was composed of plaster (light elements: 1 H, 8 O, 16 S, 20 Ca), adding heavy (absorbing) particles would create a texture. Different sizes and heavy elements have been tested ( 26 Fe, 29 Cu, 50 Sn, 20 Pb, 83 Bi). Projections of two volumes loaded with Pb and Cu powder as acquired by the lab-CT NSI software are presented figure 2.2(a,b). It can be seen that in (a), the Pb powder are extremely absorbent. The reconstruction of the Pb loaded sample is polluted with metal artifacts. Concerning the size, particules should be reconstructed at least with 4-5 voxellength. Two sizes of particles are shown figure 2.2. The Fe particles were in this case too small. Cu powder was finally selected because the reconstruction could be performed without too much artifacts.

Other techniques could have been used to enhance the texture: modifying the geometry, apply the powder or other tracers on the surface, etc.

It is noteworthy that enhancing the texture was performed for convenience. Because the plaster is a material with an inner texture (particles, air bubbles, etc.) and kinematics is highly regularized, the natural texture would have worked as well (but with less sensitivity thus a lower signal to noise ratio).

Mechanical test

The in-situ testing machine [35] is a tensile/compression machine controlled by the axial displacement of the lower part. The sample is placed on a composite tube that balances the load exerted on the specimen. The load cell capacity is 1 kN with an uncertainty of 7 N.

The experiment is composed of loading steps using a P-DVC approach (i.e. without acquisition and reconstruction of the deformed states). Previous works [21,22,23] have shown that the procedure with a synchrotron beamline gave accurate displacement field with 2 angles (chosen orthogonal). Similarly, angles for the P-DVC procedure are to be selected. A high sensitivity to a degree of freedom allows for an accurate identification of its amplitude. The sensitivity to the crack position is here chosen to be maximized. Prior to the experiment, an elastic FE simulation with simple boundary conditions was performed on a synthetic gray level volume with a microstructure representative of that of the real sample, and the norm of the sensitivity b j was estimated for all integer angles values so as to select the most sensitive one. First, the angle 0 • (i.e. in the crack front direction) gave the highest sensitivity due to an important orthogonal displacement around the crack. However, because of projection, some information (i.e. displacement along the projection axis) was lost. Therefore, the second angle was selected at 90 • to provide the lacking sensitivity.

Supplementary materials

One key quantity in the measurement procedure is the projected sensitivity field Π k [f ] The pre-computation of the sensitivity has been performed on a synthetic volume composed of a light matrix and heavy absorbent particles (3D Gaussian). The volume was generated from a matlab function due to Dr. Jan Neggers. The sensitivities were computed using finite differences with a matlab-FEM script and were used to weight the 3D gradient before projecting it with a guess of the conic projection operator. The projection of the sensitivity at 0 • with respect to the symmetrical advance of the cracks is shown figure 2.3(a). The most important part of the residual is located in the central part of the sample and is the signature of an horizontal motion. The norm of this sensitivity with respect to the 360 degrees angle is plotted figure 2.3(b). It was chosen to optimize the angles from the crack sensitivities thus 0 • and its complementary 90 • were selected. Other criterion on the sensitivity could have been introduced. For example maximizing the minimal eigenvalue of the Hessian matrix, or its determinant, would also have been a possible criteria.

Eight radiographs (acquired with an average over 50 frames) at each loading state were captured at every multiple of 45 • (8 were chosen beforehand for caution but only two are used in the following procedure, namely at 0 • and 90 • ). The acquisition time for each of these steps is 35 s. The test is composed of two parts:

• 7 loading steps without crack (i.e. in an elastic regime), from 0 to 600 N, spaced by approximately 100 N

• 6 loading steps with crack initiation and propagation.

The fast 8 radiograph acquisitions are performed after approximately 15 minutes hold time to avoid possible relaxation phenomena. This procedure (again motivated by caution) revealed unnecessary. At the end of the experiment, a final 3D scan (i.e. with a full acquisition of 1500 radiographs) was performed. This reconstruction was designed to track possible non-planar crack propagation, and for possible validation of the P-DVC results. Table 2

Sinograms

In order to control the reconstruction as well as the projection processes, the reconstruction of the reference state is performed with the ASTRA toolbox [36] and cone-beam projections. Because of computation time and GPU memory requirement for ASTRA, the reconstruction from an FDK algorithm is carried out with different sinogram scales:

• Scale 2 : Gaussian filtering of the sinogram with a characteristic length of 2 pixels followed by a down-sampling of 1 out of 2 pixels in both directions, so that at each angle, the sinogram is a 751 × 972 pixels image, and hence after reconstruction the effective voxel size is doubled to 30 µm. This scale is used for the reconstruction shown in fig. 2.1(a).

• Scale 4 : Gaussian filtering of the sinogram with a characteristic length of 4 pixels followed by a down-sampling of 1 out of 4 pixels in both directions, resulting in 372 × 486 pixels images). This coarse scale is used for development and initialization of the finer scale.

The multi-scale procedure is applied on the deformed sinogram too.

As can be seen in figure 2.4, the specimen contains a large air bubble on one of the top corner. This bubble is meshed in the following finite element computation in order to account for its possible influence on the crack propagation.

A bounding box around the sample (e.g. white dashed rectangle in the projections fig. 2.5) on the residual is chosen to not bias the analysis with the still composite tube and avoid the top and bottom artifacts on the loading plates. Finally the central hole and notches are hidden by a mask (delineated by full white lines) on the projections because of important phase contrast. The studied area, inside of the bounding box and excluding the central mask is called region of interest (ROI). It is defined in 3D, and its projection is considered on the sinogram (the exterior of the bounding box is considered to be void, but the excluded inner parts mask their entire shadow). The sinogram and ROI for step 8 are shown in figure 2.5. The crack may be visible at θ = 0 • due to a phase contrast built along the optical path. This is to be contrasted with the following P-DVC procedure that estimates the position of the crack from the displacement field.

Residual fields

Artifacts residual fields

As shown in section 2.2.2, the P-DVC code minimizes the quadratic norm of the 2D residual fields ρ(r, θ k , t) ≡ s(r, θ k , t) -Π k [f (xu)]. It is interesting to extend this procedure to the initial state t = 1 where u is identically null per definition. Due to measurement noise and slight inaccuracies of the reconstruction and re-projection algorithm, the artifact residual ρ(r, θ k , 1) ≡ s(r, θ k , 0) -Π k [f (x)] is not 0. The initial artifact residual is shown in fig. 2.6 with a divergent color map to highlight differently positive and negative values. All residual fields are presented in percent of the initial projection amplitude, after cutting out the highest and lowest 5% of gray levels. In this field, some features could be clearly interpreted as detector artifacts and some due to phenomena neglected in the reconstruction (beam hardening, phase contrast, metal artifacts). It was observed that this artefact residual field could still be read in the following residuals ρ(r, θ k , t). Because ρ(r, θ k , 1) appears to be a limit than cannot be overstepped, in the following, we will focus on the difference ∆ρ(r, θ k , t) = ρ(r, θ k , t) -ρ(r, θ k , 1), in order to minimize the effects of this initial systematic bias. where the standard deviation is computed over the clipped and unmasked region of interest in projection space.

Mesh and regularization

As previously explained, the displacement field is expressed as a combination of fields related to the sensitivity to the selected parameters. These fields are computed based on a mechanical model, here implemented with a finite element code. Because gypsum is an elastic brittle material, the chosen mechanical behavior is an elastic law (2 parameters: the unitary Young modulus, E = 1 MPa, conventional value chosen for identification purpose and the Poisson ratio, ν = 0.2) and a possible crack propagation that "splits" the volume into two parts.

According to the reconstruction of the deformed state at step 13, the crack path is observed to remain flat in the mid-plane of the sample. Hence, as shown in figure 2.8 the mesh is composed of two conforming parts that can be tightened (if there is no crack) with a common displacement on facing nodes at the interface. The mesh is defined on the reference and deformed images. It has 1343 nodes and 4248 T4 elements. The large bubble observed in the corner is meshed as it may affect the specimen stiffness. In order not to leave 'floating voxels' in the advection of the reference volume, the mesh is very slightly larger than the specimen volume. The effective cross-section of the mesh is about 4 % larger than that of the volume itself.

CHAPTER 2. CRACK IDENTIFICATION FROM MULTI-VIEW TOMOGRAPHY

The chosen parameters to control the kinematics are • 12 degrees of freedom that represent the applied loading (i.e. 6 dofs on each of the top and bottom faces). They describe a rigid body motion of each of these faces. Equivalently, they can be reorganized to describe the 6 degrees of freedom of the mean rigid body motion of the specimen, and 6 others that correspond to different elementary loadings: one tension and one torsion along the z axis, two shear and two bending modes (x and y axes). The sensitivity to the latter six boundary conditions are computed with a unitary displacement because of the linearity of the elastic law. These modes are shown in figure 2.8.

• 2 degrees of freedom for the crack front position. The crack is assumed to be plane. The sensitivity for this parameter is the difference in displacement field due to a crack propagating from a crack length of p i (with i equal to 13 and 14 for the two crack fronts) to a slightly longer one p i + dp i , normalized by the increment dp i .

Figure 2.8: T4 Mesh adapted to the reconstruction for the six degrees of freedom applied on the upper face (while the bottom one is held fixed). On the first row: two shear and one tension modes, on the second row: two bending and one torsion modes. The color corresponds to the displacement amplitude

Results

Methodology

The results of the P-DVC procedure are now presented. First, in Section 2.4.2, the loading steps are considered independently in order to identify the displacement field boundary conditions and crack position. In this part, the first loading step where the crack has propagated, step 8, has been selected to display the results in more details. In Section 2.4.3, once the different displacement fields have been independently evaluated for each loading step, it is possible to identify the Young's modulus. Section 2.4.4 shows a refined evaluation of the displacement fields and crack position taking into account the previously determined Young's modulus and the measured load. Finally in Section 2.4.5, a minimum bound for the toughness is evaluated using the 7th load step, just before crack propagation.

Independent analysis of load steps

The same analysis was performed for the 12 loading steps, and the entire series will be presented for some global results, but when focusing on one loading step, the 8 th one is here systematically chosen, since it is the one of crack inception. Convergence criteria are defined on both variations of the SNR and displacement field (∆SNR < 10 -3 and ∆u(x, t) 2 < 0.3 µm or 0.01 voxels at scale 2) and the computation stops when both conditions are satisfied. The P-DVC procedure is observed to converge in 4 to 6 iterations depending on the loading step, as exemplified in Figure 2.9 for step 8. It can be seen that the residual field values are much smaller and smoother at the end of the process than in the initial state shown in Figure 2.7. The previous alternating positive and negative patterns (signaling a displacement) have been erased, confirming thereby that the kinematics has been well captured. The crack path is much more apparent than on the initial residuals. Even the central part of the sample, masked in the ROI, has been well corrected because of regularization. Some horizontal and vertical lines are clearly visible. They correspond to detector artifacts that are not corrected by our procedure.

The displacement field in the three directions is shown figure 2.11. This displacement is composed of a large vertical rigid body motion of 4 voxels (due the positioning of the sample and to the low rigidity of the testing machine), a compressive state and a bending The two crack front positions were two degrees of freedom of our identification. At step 8, they indicate that both top and bottom cracks have propagated up to the top and bottom faces of the sample. This was unexpected as the DCDC geometry should lead to a stable propagation for a displacement control, but the low rigidity testing machine may be responsible for releasing an additional energy at the onset of propagation.

P-DVC computed at step 7 leads to no crack inception, (i.e., no further than the notch), and indeed, in the residual before and after convergence as shown in Figure 2.12, no trace of a crack is apparent in contrast to step 8. This observation is interesting because it provides a lower bound of the material toughness (see section 2.4.5). Moreover step 7 load is only 19 N smaller than that of step 8. Other loading steps give similar results for the improvement of the residual fields. The SNR ratio before and after the procedure for each step is shown figure 2.13. As expected, the initial SNR decreases with the loading step as the load and displacement increase. Even the initial residual at step 2 can be well improved although it mostly corresponds to a rigid body motion.

It can be seen that after the crack propagation, at step 11 to 13, the SNR is lower than before. With the reconstruction of the final 3D volume (step 13), one may see many cracks due to damage at the top and bottom ends of the sample, as well as flexure cracks transverse to the sample axis at the position of the central hole. Such damages are not described in our modeling and hence they induce a degradation of the registration quality and the SNR decreases.

Identification of elastic properties

The Young's modulus identification has been performed with a space-time regularization considering all previous analyzed loading steps. The measured forces are compared to the computed ones at convergence of the 14 kinematic parameters. As shown previously, the first 7 steps are considered to be in the elastic regime (no crack propagation), hence they should give accurate estimate of E without significant model error. The following steps (including the through crack) should be amenable to a similar elastic analysis until a significant further damage appears. Exploring different ranges for such an analysis lead to the conclusion that beyond step 11, the elastic model was not appropriate, and hence the Young's modulus was estimated based on the first 10 loading steps (including the computed crack extension). Considering the force uncertainty, the uncertainty on the Young modulus is 1.5 %. From Eq. 2.18, the estimate is E = 1.12 GPa. This value is low compared to the literature [37,38] (2 to 5 GPa) but corresponds to ex-situ pre-tests on the same plaster and geometry (three tests with a Young modulus of [998, 931, 940] MPa. A high water/plaster ratio was chosen for an easy moulding of the specimen. Inserting back this value into the elastic model provides a computed force that can be compared to the measured one as shown in Figure 2.14.

Further regularization using identified elastic property and force signal

Once the Young's modulus has been identified from all loading steps, it is possible to update the displacement field considering in the P-DVC procedure the additional constraint given by the force measurement (equation 2.14). This new 4D measurement gives a revised displacement field that can be compared to the previous one. The comparison of the norm of the nodal displacements for the 10 first steps are shown figure 2.15. The maximum difference, 0.3 %, is small meaning that the kinematics was well captured even prior to accounting for the force minimization. A small difference is visible in the 2 last considered steps due probably to less accurate model assumptions.

Toughness evaluation

Because step 8 corresponds to crack inception, the analysis of step 7 gives a lower bound for the sample toughness. In the same spirit, the loading of step 8 applied on a non-cracked specimen gives a higher bound of the toughness. Moreover, step 7 is only 19N less than step 8 and the crack geometry is known (notches) hence the bounds should be accurate. The G f parameter is obtained with an elastic computation performed in Abaqus ® using the J-integral method. The two bounds are G f = [10.4, 11.1] J/m 2 . This result can be compared with standard G f measurements performed without full field analysis. In these standard procedures, only the compressive force is considered (hence the displacement field is assumed to be a uniform compressive state). In our experiment, due to an important bending of the upper face, it is important to consider full field measurements. When the standard formula (ignoring bending) [39] is applied to the present case, it leads to G f = [9.5, 10.1] J/m 2 , corresponding to 9% difference with These values of G f are high compared to the literature [37,38]. However, the low value of the identified Young modulus gives K Ic = [0.108, 0.111] MPa.m 1/2 that is in the low range but compatible with literature values ranging from K Ic ≈ 0.1 to 0.5 MPa.m 1/2 .

Discussion

This study is a first attempt to exploit a complete in-situ mechanical test using P-DVC, in a lab-tomograph. This method is 4D, in the sense that one has access to the entire displacement field, u(x, t), in space and time at each loading step. The fact that only two projections (per step) are needed endows the methodology with a much enhanced time resolution, by a factor of several hundreds.

This opens new horizons for biological tissues where a low dose of x-ray radiation is needed to limit radiation damage. However, one should emphasize that it is necessary that a reference volume has been reconstructed. Saving on x-ray dose only comes from the subsequent tracking of the motion in time where the additional dose is very small. Such a property may be used to add, at low additional cost, It also allows considering time dependent phenomena that could not be captured without P-DVC, prosaically because of the evolution of the sample during the scan time (several hours or more). This may prevent a proper tomographic reconstruction and hence preclude the usage of volume correlation (classical DVC) to track the motion, or any other evolution of the sample. Hence, on this ground, the superiority of P-DVC over DVC is obvious, as it renders possible the study of situations that were simply out of reach with classical means.

Complex specimen shapes or microstructures can be handled by this technique, without restrictions. It is also to be emphasized that when the geometry of the sample evolves during the test (as in the present case where the crack was not present in the initial geometry but appeared along the test), one may resort to the final state to extract its very precise geometry. More generally, possible non-planar cracks could be considered in their final geometry and cohesive element implemented in the model in order to progressively open the crack as called by the studied projections.

In the present study, a simple elastic model revealed appropriate (as judged from the residuals). In the case of a more complex mechanical behavior measurement and identification, a more sophisticated model could be used. In the same spirit as the crack propagation that is a non-linear parameter, the finite element model as well as projected sensitivities may have to be re-computed when needed and possibly at each iteration of the procedure, that would increase the computation cost of the analysis.

Numerous artifacts are known to be detrimental to tomographic imaging. When working with raw projections, a number of possible corrections can be envisaged and tailored to the specific conditions of acquisition of any given projections. In fact reading projection residuals is typically easy to indicate if motion, or intensity modulation, or any other phenomenon, be it localized in space or uniformly distributed, is needed to reduce their level. Exploiting the above flexibility allows for a very precise adjustment of artifact corrections that can hardly be incorporated prior to reconstruction as this would increase considerably reconstruction time.

Noise is also a key feature often limiting the quality of reconstruction. Long-range correlations, present in the latter, are the result of the reconstruction that processes the noise. However, in the projections, noise is much less complex, and hence, taking it into account by a relevant weighting of the input data is easily accessible. This is expected to enhance the reliability and lower the uncertainty on projection registration.

The limits of the proposed technique mostly comes from our ability to model the test with a reduced set of parameters. It is essential that no (or very small) model error is introduced, otherwise one may not be in a position to approach the actual projections, from the deformation of the reference volume. The appropriateness of the model can be read in the fact that the residuals may be reduced to a low level that can be ascribed to noise (or to small amplitude artifacts that would not warrant corrections). Another potential limitation is the presence of large displacements that would violate the relevance of small perturbation analysis. Although a multiscale approach may help correcting large amplitude motion, it is wise to perform acquisitions of projections at small increments of loadings.

In terms of acquisition time, one may consider a much faster time lapse than the one used in the present study. First 8 angles were captured but only 2 were used. Moreover a very long sequence (average over 50 projections) was acquired and is today not the limiting factor for the proposed analysis. 10 projections are certainly sufficient, so that the two projections could be acquired within about one second. One limit here may be the time needed to rotate the sample in between the two views.

An appealing perspective to achieve a higher temporal resolution is to use the time regularization to acquire a single radiograph at each loading step and perform a continuous rotating test. With such a 4D procedure, the unknowns contribute for the mechanical response of several steps, and with a suitable design they can be determined with the entire time series, challenging the computational aspects of the inversion.

Conclusion

The analysis of a complete DCDC test with 13 loading steps was performed successfully allowing for the analysis of the entire kinematics, and for the identification of the Young's modulus as well as a bound on the critical elastic energy release rate. Identification of the Young modulus has been performed considering all loading steps at once, a specific feature made possible with the 4D space and time proposed procedure.

This analysis was performed taking the actual boundary conditions of the test rather than assuming an ideal one which would have led to 9% error on the toughness. In this measurement of a quasi-brittle sample, a simple elastic model with a crack propagation has been chosen. This kinematic model has been validated because the low residuals showed that the kinematic was well captured. The kinematic regularization (based on only 14 parameters) has been chosen considering the studied sample and test. A scan of the cracked specimen at the end of the experiment permits to design the mesh according to the real geometry.

This basis could be enriched by additional freedom such as shrinkage of the bottom and top part due to a Poisson effect, or a more complex crack path and front. However, phase contrast may induce specific features in the projection of the crack faces that are not modeled at present and that may limit the resolution of a much finer description.

This developed P-DVC procedure in a lab-tomograph permits, from only one reference volume, to measure the displacement field of a "moving" specimen from two radiographs. P-DVC required 350 times less acquisition than standard DVC procedures. Further elimination of extraneous (unused) projections and performing fewer frame averages for radiographs, would further increase this ratio well above three orders of magnitude.

Bibliography

[1] E. Maire, J.Y. Buffière, L. Salvo, J.J. Blandin, W. Ludwig, J.M. Létang. On the application of X-ray microtomography in the field of materials science. Advanced Engineering Materials, 3(8):539-546, 2001.

[2] L. The measurement of 4D (i.e., 3D space and time) displacement fields of in-situ tests within X-ray Computed Tomography scanners (i.e., lab-scale X-CT) is considered herein using projection-based Digital Volume Correlation. With one single projection per loading (i.e. time) step, the developed method allows for loading not to be interrupted and to vary continuously during the scan rotation. As a result, huge gains in acquisition time (i.e., more than two orders of magnitudes) to be reached. The kinematic analysis is carried out using predefined space and time bases combined with model reduction techniques (i.e., Proper Generalized Decomposition with space-time decomposition). The accuracy of the measured kinematic basis is assessed via gray level residual fields. An application to an insitu tensile test composed of 127 time steps is performed. Because of the slender geometry of the sample, a specific beam space regularization is used, which is composed of a stack of rigid sections. Large improvements on the residual, whose SNR evolves from 9.9 dB to 26.7 dB, validate the procedure.

Introduction

The development of Computed Tomography (CT) has been a major breakthrough in materials science [1,2,3,4], providing nondestructive measurements of the 3D microstructure of imaged samples. Initially developed using synchrotron sources, it is now a common equipment accessible in laboratories. In addition to the broad range of spatial resolutions (from nanometer to few meter scans), significant progress has been achieved in ultra-fast X-ray tomography as illustrated with the 20 Hz full scan acquisition for the study of crack propagation [5].

Coupled with mechanical tests, these acquisitions become an attractive tool in experimental mechanics. First used for ex-situ tests where the materials are deformed outside of the tomograph [6], the recent evolution of testing machines and CT makes in-situ tests possible [7,8,9]. In the latter cases, the sample is imaged inside a tomograph, either with interrupted mechanical load or with a continuously evolving loading and on-the-fly acquisitions. Visualization of fast transformation, crack opening, or unsteady behavior become accessible. Combined with full-field measurements, in-situ tests coupled with X-ray CT offer a quantitative basis for identifying a broad range of mechanical behavior.

A robust and accurate method to quantitatively measure kinematic data from the acquired images is Digital Image Correlation (DIC) in 2D and its 3D extension, Digital Volume Correlation (DVC). The latter aims at capturing the way a solid deforms between two states from the analysis of the corresponding 3D images. The measured displacement field is then used to calibrate model parameters from inverse problem procedures (e.g., finite element model updating [10], virtual fields method [11]). The more numerous the acquisitions, the more accurate and sensitive the identification procedure. DVC methods have hence been developed in a 4D space-time framework [12,13,14] using all available volumes globally. In 4D analyses, an acquisition of 5 to 15 steps is usually performed. However, the major limitation of CT imaging, especially in lab-tomographs is the acquisition time. To give orders of magnitude, each reconstructed CT volume, being composed of about one thousand 2D radiographs, takes approximately one hour to be acquired. This duration limits the number of possible acquired scans and hence restricts such identifications to time-independent behavior. Moreover, spurious motions may occur during acquisitions, blur the projections and reconstructions. In some cases it is required to wait 0.5 to 1 hour for relaxation or creep stabilization at each loading step, before each scan [15,16].

The recently developed Projection-based DVC (P-DVC) [17] is an interesting method to circumvent these difficulties. Instead of working with a series of 3D volumes, it is proposed to directly measure the 3D displacement field from few of those radiographs. This procedure exploits the property that two projections at a given orientation of the sample under different loads contain a partial information about the full 3D kinematics, and sampling few angles may be sufficient to extract the entire displacement field. In a similar spirit, and very recently, the 3D tracking of rigid grains was proposed using few radiographs [18]. Previous works have developed the P-DVC method for the analysis of a cracked cast iron sample imaged with a synchrotron X-ray source [19,20]. It was shown that the measurement was possible with only two orthogonal projection angles, thereby allowing for a huge gain in acquisition time (i.e., a factor 300). In Ref. [21], a spatiotemporal framework was proposed in order to analyze the kinematics with a single projection per angle, provided a suited temporal regularization was used. This methodology was applied to a lab-tomograph in-situ experiment where crack propagation in plaster was monitored and quantified (via double cleavage drilled compression test).

In order to deal with space and time displacement field, a Proper Generalized Decomposition (PGD) framework [22,23,24] will be used. The displacement field is decomposed over a basis of separated functions, namely, spatial modes and time (or projection angle) modes. The spirit of the approach, which is inspired from the PGD technique, consists of a progressive enrichment of the space-time modes for displacement corrections. Modes are progressively added until a convergence criterion based on the residual field is reached. First introduced in DIC and DVC [25,26] to provide a reduced basis for the displacement field (separating different spatial directions), PGD has been extended to the P-DVC framework for the calibration of material parameters in a lab-tomograph using a single projection instead of two [27]. In that work, deviations from a perfect rotation were taken into account over one scan in order to reach a good quality tomographic reconstruction. The spurious motions were described as small amplitude rigid body motions over time in addition to the ideally expected rotation. In the same spirit, a developed projection-based approach allows the imbibition process to be followed in sandstone at a frequency of 5 Hz using Neutron tomography [28].

In the sequel it is proposed to measure 4D (space and time) displacement fields based on a series of projections acquired at different angles. Using only one projection per angle allows the sample to be loaded and continuously rotated without hold and/or dwell time. This reduction in the number of radiographs, which is offered by the mathematical formulation, together with the now tolerable continuous loading leads to half the acquisition time of the previous method. The analysis of a complex kinematics, up to localization, with such an approach is a novelty. In Section 3.2, a highly regularized PGD framework coupled with a P-DVC procedure is introduced in order to capture the kinematics of the studied sample.

Section 3.3 focuses on the specific formulation of the approach to slender samples, and the space and time bases are discussed. Section 3.4 is devoted to the analysis of a tensile in-situ test on nodular graphite cast iron with the proposed methodology. The experiment, which is composed of 127 time steps, is performed in 6 minutes with a continuous loading until failure. The measurement of the displacement field for each time step is presented and shown to provide very low registration residuals. A gain in acquisition time of more than two orders of magnitude compared with standard methods is obtained.

Full-field measurement

Projection-based DVC

The proposed approach to fast 4D (space and time) measurements is called Projectionbased Digital Volume Correlation (P-DVC) [17,20]. Instead of working with reconstructed volumes as in standard DVC (whose acquisition time is one of the major limitations of CT in laboratory facilities [29]), it aims to measure the 4D displacement field from a series of 2D projections acquired at different angles θ(t) and loadings.

The registration operation consists of minimizing the sum of squared differences between N θ 2D projections g(r, t) of the deformed configuration at different times t, or angles θ(t) and loading steps. The procedure makes use of the 3D reference image, F (x), which is reconstructed using classical means. It provides for all voxels of the 3D space x = (x, y, z) the local X-ray absorption coefficient. This reference volume, corrected by the displacement, u(x, t), and projected with the orientation θ(t) should coincide with the acquired projections, g(r, t) when the displacement, u(x, t), is correctly measured. In other words, introducing the reference volume deformed by any trial displacement field, v,

F v (x) ≡ F (x -v(x, t)) (3.1)
and the so-called residual field, ρ(r, t; v),

ρ(r, t; v) = Π θ(t) [ F v ] (r, t) -g(r, t) (3.2)
the 3D displacement field, u(x, t), is sought as the minimizer of the following cost function

χ 2 [v] = r,t ρ(r, t; v) 2 (3.3)
where Π θ(t) is the projection operator in the θ(t) angular direction, and r = (r, z) the coordinates in detector space. The integrand of the previous functional can be linearized considering small displacement field corrections δu compared to the microstructure correlation length

χ 2 [u + δu] = r,t ρ(r, t; u) -Π θ(t) [δu • ∇ F u ] (r, t) 2 (3.4)
with ∇ the 3D gradient operator. It is noteworthy that after each evaluation of the displacement corrections δu from a known displacement u (n-1) such that u (n) = u (n-1) + δu, a correction of the volume F u (x) is performed so that the previous equation is used without approximation. The P-DVC framework requires the acquisition of one reference volume in order to compute the correction terms. The latter is to be performed ideally in the same conditions as for the experiment, but without load (or with a very modest one). This only one classical tomographic scan is generally not challenging.

In order to validate the proper evaluation of the displacement, one should consider the magnitude of the residual field that highlights all projection differences that are not captured by the measured displacement field (e.g., noise, artifacts of the detector, illconvergence, model error). Ideally, it should be statistically indistinguishable from noise. The Signal to Noise Ratio (SNR) can also be defined to evaluate the residual quantitatively. The higher the SNR, the better the solution. It is commonly defined as 20 times the decimal logarithm of the ratio of the standard deviation of the projections over that of the residual fields [30].

Regularization procedure

CT images are defined by a huge number of data (usually billions of voxels treated independently). However even if the displacement field is to be defined for all voxels, its regularity legitimates the use of much less degrees of freedom so that it is sought in a vector space generated by a reduced kinematic basis. In the following, it is proposed to choose as reduced basis the product of separated space Φ j (x) and time σ i (t) fields for which the sought amplitudes are

u ij u(x, t) = Nt i=1 Ns j=1 u ij σ i (t)Φ j (x) (3.5)
with the space and time dimensionality respectively N s and N t . Such separated expression is standard practice and implies no restrictions (provided N m = N s × N t is large enough), but it will reveal convenient for the following model reduction technique. The previous displacement is written in a different way

u(x, t) = Nm l=1 Nt i=1 a l i σ i (t)   Ns j=1 b l j Φ j (x)   (3.6)
Hence the previous amplitude matrix is

u ij = l a l i b l j .
Different regularization procedures of the displacement field have been introduced in the literature for global DVC when the kinematics is discretized over finite element meshes [31]. Spatially, local (elastic) equilibrium-gap penalty was proposed [32,33,34]. Other model-based regularizations can be used to calibrate model parameters [35,13,14,29]. These regularizations lead to drastic reductions in the number of unknowns and enable for seamless experimental/numerical interfaces. They will not be considered herein. A more generic approach will be followed hereafter.

Greedy approach to P-DVC

The entire problem composed of a large number of unknowns (i.e., 630 degrees of freedom in the present application) may be costly. The method proposed to solve the minimization problem is the Proper Generalized Decomposition (PGD) approach, where modes are successively determined and added as long as the residual level is considered too high to be explained by noise, which corresponds to detector artifacts or if the noise is assumed to be white and Gaussian: χ[u] > 2γ f N pix , with γ f the standard deviation of the noise and N pix the total number of pixels on which the procedure is applied).

Proper Generalized Decomposition techniques [22,36,23,24] consist of successive enrichments of the displacement field u(x, t) adding a new contribution at each iteration, each term of the sum being sought a priori in a separate representation. PGD-DIC and PGD-DVC [25,26] with one-dimensional space functions are here extended to 3D space and time (i.e., 4D) analyses.

In the following progressive PGD procedure, the time and spatial modes are identified successively, one per iteration, with a greedy approach [37]

u l (x, t) = u l-1 (x, t) + Nn i=1 a l i σ i (t)   Ns j=1 b l j Φ j (x)   (3.7)
Let us note that only the product a l n b l j matters so that a convention such that {a l } = 1 or {b l } = 1 can be freely chosen without consequences. A fixed point algorithm is used to get the solution. Alternate minimizations of the two unknown column vectors {a l } and {b l } are proposed. The minimization of the functional leads to the determination of the unknowns with two coupled equations

{a l } = Argmin {a} (χ 2 u ({a}, {b})) (3.8)
{b l } = Argmin {b} (χ 2 u ({a}, {b})) (3.9)
i.e., the minimization of χ 2 u with respect to the additional mode is considered. It is proposed to write the projected sensitivity for each degree of freedom of each section as

S i (r, t) = Π θ(t) [Φ i (x ) • ∇F (x )] (3.10)
and the associated matrix

B ij (t) = r S i (r, t)S j (r, t) (3.11)
where x is the corrected position of any voxel x with the previously identified modes such that x = xu (l-1) . The sensitivity S i (r, t) is composed of degrees of freedom i, pixels of the detector r and angles

θ(t) hence it is of size [N s × N p × N θ ].
Computing and storing this matrix is the longest operation in the procedure. The other quantities are easily obtained from combinations of these projected 3D sensitivity fields. It can be noted that S i (r, t) should be recomputed for each identified mode while x is updated. Because the degrees of freedom per sections are quite independent, the sensitivity matrix can be stored as a sparse matrix (and could be highly parallelized if the beam is non-diverging).

The two parts are obtained from the above linearized integrand using Newton's scheme. The derivative with respect to {b l } leads to

{b l } = [N l ] -1 {n l } (3.12)
where [N l ] is the spatial Hessian matrix of χ 2 u with respect to {b l } (i.e.,

N l ij = ∂ b i ∂ b j χ 2 u
) and {n l } the second member vector based on the residual field, which is written as a weighted sum over time of the sensitivities with α l (t) = Nn i=1

a l i σ i (t) N l ij = t α l (t)B ij (t)α l (t) (3.13)
and

n l j = r,t ρ(r, t, u l-1 )α l (t)S j (r, t) (3.14)
Similarly, the derivative with respect to {a l } leads to

{a l } = [M ] -1 {m} (3.15)
where, as previously, [M l ] is the temporal Hessian matrix (i.e., M l ij = ∂ a i ∂ a j χ 2 u ) and {m l } the second member vector based on the residual fields

M l ij = t σ i (t)b l n B nm (t)b l m σ j (t) (3.16)
and

m l j = r,t ρ(r, t, u l-1 )σ j (t)b l n S n (r, t) (3.17)
A general overview of the 4D PGD P-DVC procedure is shown in the algorithm 2. Even though a maximum value of iterations or convergence criteria, p and α , can be enforced to stop the fixed-point algorithm, this revealed unnecessary as the maximum number of iterations to reach stagnation is usually quite low (i.e., 3-5).

Application to slender samples

It is proposed to analyze the kinematics of a nodular graphite cast iron dog bone sample during in-situ tension until failure (Figure 3.1(a)). As many of the standard in-situ test geometries, the proposed sample is slender and may be considered as a beam with one dimension much higher than the other two. An appropriate kinematic regularization based on a very small number of degrees of freedom is thus proposed.

Slender-shaped specimen kinematics

The slender geometry is axially divided into N s undeformable beam sections, normal to the beam axis, with trilinear interpolations. The displacement field is hence written in this reduced basis

Φ i (x) = Ns j=1 6 k=1 p j (z)ψ k (x) (3.18)
with i = (j, k), ψ k (x) the 6 rigid body motions of the entire sample, p j (z) the shape functions that allow specific sections to be selected and interpolations to be performed between them. For a linear interpolation, the shape function will be triangular functions whose maxima are located at the section positions and its length will correspond to the size of two sections. It can be noted that the shape functions are the same for all sections (thus do not depend on mode j). The consequence is some coupling between rotation and translation for the section far from the center of rotation. In the present application, the beam is composed of 15 cross-sections (see Figure 3.1(c) for the depiction of few sections on the slender sample used in the application), each animated with 6 Rigid Body Motions (RBMs, i.e., 3 translations and 3 rotations).

Time dependence

Generically, in the time dimension, σ i (t) are chosen functions that also introduce temporal regularizations [38,39]. The availability of additional measurements, such as tensile loads, T (t), also offers the possibility to define temporal variations that may depend on load rather than time, or any combination that may be physically motivated. This is a very convenient and non-intrusive way of incorporating some information, or some possible relationship, in the measurement parameterization. Ultimately, no time regularization can be introduced if wanted by choosing as elementary functions σ i (t), functions valued 1 only for time step t = t i and 0 otherwise. Other basic examples are provided by tent shaped (piecewise linear) functions, polynomials, Fourier modes and splines.

In the treated application, the time basis does not change with mode identification thus does not depend on i, but such cases could be designed (if the goal was for example to extract first the elastic part of the test).

Test case

Tomography acquisitions

The application case for this study is an in-situ tensile test on a nodular graphite cast iron sample (similar sample geometry can be found in Refs. [40,13,14]). The geometry of the sample is described in Figures 3.1(b) and 3.4.1. The central part was thinned with a radius of 20.5 mm in order to ensure that the specimen would break in the ligament area and not in the grips. The sample was mounted in an in-situ tensile testing machine similar to that used by Buffière et al. [41] (see Figure 3.4.1) and was scanned in the LMT equipment (X-View X50-CT, North Star Imaging, 180 kV, 130 µA, W target). The voxel size (using 4 × 4 binning at the acquisition stage) was set to 10.7 µm.

Two flat-fields and one dark-field were acquired after conditioning and before the experiment in order to perform flat-field and dark-field corrections. Each radiograph was averaged with 5 frames in order to reduce acquisition noise without loosing too much time (it is noteworthy that a single frame would have been sufficient inasmuch as noise were white and Gaussian). The radiographs had a definition of 954 × 768 pixels and an initial crop of the edges allowed its size to be reduced to 954 × 432 pixels (Figure 3.1(a)). All the projections used herein have been normalized to 0 for air, and 1 for the maximum log-attenuation, after dark-field and white field corrections as well as beam-hardening corrections. In order to perform a multiscale approach, one lower scale was used in the following procedure. Coarse graining of 2 × 2 elementary pixels into one superpixel was carried out with the convolution of the projections by a Gaussian kernel with a characteristic width of 2 pixels, and downsampling to a coarse 2 × 2 regular square grid to create smaller images (called images at scale 2). The projections were obtained after flat field normalization and standard (i.e., third order polynomial) beam hardening corrections [42] due to the high absorption of the ferritic matrix. Reconstructions and projections were performed with the ASTRA toolbox [43], using the Feldkamp-Davis-Kress (FDK) procedure suited for cone beams [44]. The initial projection f (r, θ) were compared with

In-situ tensile test

The fast space-time in-situ experiment is composed of three phases:

1. loading to T = 250 N in order to remove the backlash that would introduce rigid body motions, 2. a complete scan of the reference state (at 250 N) that consisted of 600 radiographs captured at equally spaced angles ranging over a full 360°rotation. This scan took 22 min to be acquired, 3. continuous rotation of the sample with 50 projections per full rotation at a rate of one projection every 2 s. One hundred twenty seven projections were acquired during 2.5 full rotations ( as illustrated in Figure 3.3). The first full rotation (i.e., 50 time steps or 100 s) was performed at constant load and was used to quantify the uncertainty. The remaining rotation (starting after 100 s) was carried out with a continuous load change (from 250 to 750 N), as shown in Figure 3.4, controlled at a constant stroke velocity of 2 µm/s. The particular choice of these angles will be discussed after the analysis of the results. The 60-pixel wide right and left edges of the radiographs and 100-pixel long top and bottom parts were discarded. This operation avoided the top and bottom parts whose quality was low because of the divergent X-ray beam and reconstruction process. The temporal basis used to measure the rigid body motions is composed of N n = 7 time functions σ n (t) (Table 3.1). The sine and cosine functions are introduced as they are expected to occur for a slight misalignment of the sample with respect to the rotation axis [27]. Table 3.1: Temporal basis chosen for the kinematic measurement.

σ 1 (t) σ 2 (t) σ 3 (t) σ 4 (t) σ 5 (t) σ 6 (t) σ 7 (t) T (t) sin(θ(t)) cos(θ(t)) t 0 = 1 t 1 max((t -50), 0) max((t -50), 0) 2
The initial residual field is shown in Figure 3.5 for the three selected angles at different time steps and loads of the procedure with a divergent color map that highlights the positive and negative patterns (so as to ease the visual interpretation of a residual displacement). Let us stress that here and in all subsequent figures showing residuals, the gray level normalizations of the initial projections between 0 and 1 have been preserved so that the color bar values can be compared. Large levels of the residual field are observed at the edges of the sample corresponding to large rigid body motions. The initial SNR is 9.94 dB.

Rigid body motion measurement

It is first convenient to start erasing the mean RBMs that occurred during the mechanical test, which are due to compensated backlash or the testing machine compliance. Instead of studying each section independently, it is proposed to use p j (z) = 1. Six spatial degrees of freedom (i.e., the 6 RBMs) are hence measured for each mode (i.e., the three translations and three rotations). This makes a total of 42 unknown amplitudes. The PGD approach consists in sorting out all those 42 degrees of freedom in "modes," such that the first is the "major" one, in the sense that it allows the residual to decrease by the largest amount.

Eight rigid body modes are measured at scale 2 until the root-mean-square average residual no longer decreases. Because of the successive updates of the non-linear problem (image registration) after each mode acquisition, the identified modes are not expected to be strictly orthogonal to each other. Hence, there is no reason to converge after the 6 modes that would result from the singular value decomposition of a full identification procedure. The choice of scale 2 is a matter of convenience as long as it is much cheaper computationally, and sufficient for RBM evaluations, as this is only a pre-correction. The results are then applied to scale 1. The residual fields for the three selected angles are shown in Figure 3.6. Note that the color scale is magnified by a factor of 4 since a large part of the initial residual has been erased. It is observed that the edges of the projection are composed of large vertical and horizontal scratches due to the registration process. Those areas are not taken into account in the SNR measurement. The SNR at this step of the procedure is 25.4 dB confirming this observation. What is visible from these images is that the central part has been properly corrected by the rigid body motions. The top and bottom parts are composed of moving features that will be corrected by the next procedure. The mean translation values are shown in Figure 3.7, and designated as u x , u y and u z . The vertical u z component evolves with the applied force and corresponds to the compliance of the testing machine, and the motion of the bottom part of the sample by the grip. The u x component is a transverse motion between the rotation axis and the sample, which accidentally almost coincides with the x-axis as can be seen from the fact that u y is much smaller. 2) and ( 3) respectively corresponding to u x , u y and u z .

Tensile deformations

From the previous corrected residual fields where the rigid body motions have been erased, it is now possible to measure motions corresponding to a tensile deformation. Each section of the beam is now considered as being independent. 4 deformation modes are measured (for a problem composed of 630 degrees of freedom). The residual field at the end of the procedure is shown in Figure 3.8 and reaches a mean SNR value of 26.7 dB. It is observed that the previous alternating positive and negative features have disappeared. The captured displacement field correctly reduces the residual and thus is deemed trustworthy. Some residuals are still visible in the center of the sample and may be a consequence of the shrinkage of the section expected from plastic incompressibility, but not included in the chosen kinematics. The dark residual located in the central part for the last steps (see Figure 3.8) is due to localized necking, a precursor to ductile fracture that will break the sample shortly thereafter, at step 127.

The initial SNR for each step is shown in Figure 3.9. The time periodicity that is seen on the graph, with a period of 25, is due to specific angles at which the residuals are more sensitive to the displacement field (e.g., angles where the faces of the sample are parallel to the X-ray beam). The significant improvement from the initial images shows the accuracy of the proposed 4D approach. Step 50 has an accidentally low SNR value and corresponds to the beginning of the load variation. The vertical displacement field and the vertical strain ε zz = ∂ z u z (t) for each section is shown in Figure 3.10. The mean tensile extension between the top and bottom part reaches about 15 voxels. A large plastic strain concentration (i.e., not only proportional to the force measurement) is visible in the central part as expected from the chosen geometry.

The very small displacement field in the first part [0;50] should be null. However it cannot be considered as uncertainty because some real displacements due to creep for example could appear. It can also be noted that the strain field displays some heterogeneities (see e.g., Section 4 at the end of the experiment) that was not anticipated from the geometry. Adding some regularizations (e.g., local mechanical regularization or directly from models) would provide a smoother field.

Discussion

Kinematics

The kinematics was simplified in the spirit of that used in beam theory (i.e., slender body). The displacement field was assumed to be well approximated by a rigid body motion for each section of the sample, but these rigid body motions were slowly varying along the longitudinal axis. This assumption is quite generic and believed to be applicable to many uniaxial mechanical tests. In the present study, six parameters (i.e., three translations and three rotations) for each rigid body motion defined at 15 cross-sections were linearly interpolated along the sample axis. Hence 90 kinematic parameters defined the motion at each instant of time. Seven time functions were introduced to account for the loading history, namely, a discretization involving a total of 630 kinematic degrees of freedom to be determined.

The choice of such a discretization is, by itself, a regularization. It was here designed to suit a slender body loaded along its longitudinal axis. However, for other cases, additional degrees of freedom may be introduced. Any information about the experiment may be used to regularize the displacement field and find a 'smart' 4D-kinematic basis, e.g., crack propagation, sudden motion (because of a shock), temperature variations. The number of unknowns may not be the major limitation of the method when a PGD procedure is used, yet including part of prior knowledge in a relevant kinematics is always beneficial. With such a suited parameterization of the kinematics, the approach followed in the present study can be generalized without notable restrictions. The only requirement is to keep the evolution between consecutive radiographs progressive and that during the experiment, the rotation has been large enough (a few full 360°rotations is appropriate), otherwise uncertainty may be large in some directions.

PGD

The proposed methodology relied on model order reduction, namely, PGD. By itself, PGD is not restrictive (i.e., separated form for the time and space variations is not limiting) nor does it rely on any smoothness assumption. It is meant to exploit the intrinsic "simplicity" of the problem, which is however difficult to formulate in other terms than saying that only few "modes" are needed for describing the kinematics. Additionally, this number of modes is not a priori defined. Rather, being driven to lower the residuals, one may argue that PGD is a smart technique able to pick up only those degrees of freedom, or modes, that are relevant, hence making the algorithm very efficient but not intrinsically based on smoothness. It has been shown that smoothness could independently be tuned through regularization [45].

In the present case, instead of an exhaustive analysis of all 630 degrees of freedom, PGD was used to select the relevant modes. Convergence was considered satisfactory based on residuals after 8 rigid body modes and 4 deformation modes. Thus the efficiency of this model reduction technique is very significant. Let us stress that a space and time separation accounting for the deformation of samples is quite generic and is expected to have a very broad range of applicability with a similar high level of efficiency.

Residuals

The residual field (i.e., differences at all angles between the projected and corrected volume and the deformed radiographs) that were minimized are the key information to validate the method. It gives a global and quantitative evaluation of the trustworthiness of the results and -if and where relevant -allows error areas (due to noise, imperfect convergence, sensor artifacts or model errors) to be visually localized and interpreted. In the present case, it was observed that the residuals mostly showed imperfections of the X-ray detector. From this observation, it was earlier suggested that the choice made prior to the experiment of averaging over 5 frames each radiography could have been reduced presumably to a single acquisition. Let us also note that the last residuals show localized features than can be interpreted as the inception of a crack. This observation stresses once more the value of such residual fields.

Additionally, residuals can also be quantified in terms of Signal-to-Noise Ratio (SNR). In the reported example, the mean SNR of the residual fields increased from 9 dB to 26.6 dB along with the correction by the 4D displacement field. This residual change highlights the accuracy of the results. The measured displacement field was presented for the vertical component (although all other components were captured). A tensile displacement amplitude of 15 voxels was measured at the end of the experiment and high strain levels were concentrated in the central part of the specimen, as expected from its geometry, until localization leads to ductile failure.

Conclusion

The full 4D (i.e., 3D space and time) kinematics of the 6 min experiment was captured with an extension of P-DVC that uses a model order reduction technique (PGD). Based on regularized fields relying on the slender sample geometry as well as a dense sampling in time, this method measures displacement fields from single projections at each time (or load) step of the experiment instead of reconstructed volumes in standard DVC methods.

The procedure was tested with an in-situ tensile test on nodular graphite cast iron composed of 127 radiographs with continuous load changes and rotation of the sample until failure. The experiment was carried out in a lab tomograph with an X-ray cone beam source. The entire experiment was carried out in 300 s, which is more than two orders of magnitude faster than standard methods. This performance goes together with the benefit of having a continuous (i.e., uninterrupted) loading so that load and rotation can be varied simultaneously. PGD was used to only focus on important space and time separated modes, thereby reducing the number of effective kinematic degrees of freedom, not from prior judgment, but as called by the experiment itself.

The measurement of the kinematics is a first step toward mechanical identification, which is one of the major goals of experimental mechanics [46]. Being able to measure the 4D motions of the studied sample at fast rates is a major asset. From the estimated displacement fields at each space and time steps, it is possible to use an inverse method in order to calibrate the mechanical parameters of, say, an elastoplastic model. Another variant called "integrated measurement" directly identifies the model parameters from the images, here radiographs, and is presented in a companion paper [45].

Supplementary materials

Using the same experiment, another short application was performed on a review article. The test was focused on the central plastic part. After the same RBM removing procedure, the kinematics was not chosen as a single C8 element. Because the behavior is expected to be that of a plastic hinge (i.e., uniform displacement in the top and bottom part of the sample and large strains in the center), the chosen space regularization for the vertical displacement is to use one single element (with 8 nodes of size 200×200×200 voxels), each composed of a single degree of freedom, vertical displacement, with inside a trilinear interpolation, i.e., a reduced version of a C8 element where transverse displacements are neglected (since the mean transverse translation has already been corrected, this assumption neglects the transverse strains). In the top and bottom parts, the (uniform) displacement is chosen as a constant extension from the central cube face by continuity. The element is located in the central part as shown in Figure 3.11. The time regularization is based on second order polynomials of the force measurement. A linear interpolation starting at time step 0 is added and another linear function valued 0 during time steps [0, 49] and linear during the load increase part:

• σ 1 (t) = 1 • σ 2 (t) = T (t) • σ 3 (t) = T (t) 2 • σ 4 (t) = t • σ 5 (t) = max(0, (t -49))
The 4D problem is composed of 40 degrees of freedom (i.e., 5 time and 8 space shape functions, σ(t) and Φ(x) respectively) and focuses on the region of interest corresponding to the projected cube element. Figure 3.11 shows the residual field change before and after registration. The residuals are rescaled to the dynamic range of the reference projection. It can be seen that a large part of the residual has been erased meaning that the kinematics has been well captured. The measured displacement for the 8 nodes at each loading step is shown in Figure 3.12. At the end of the experiment, the cube has been stretched in tension by about 11 voxels, or about 110 µm. The gray level residual is quantified, in the region of interest (see Figure 3.11), by the signal to noise ratio (SNR). The raw difference of projections leads to an SNR of 11.0 dB. The periodicity of about 25 time steps, which is seen in the initial residual, corresponds to the angles where the edges of the sample are aligned with the x-ray beam, namely, where the sensitivity to the radial displacement field is high. After the transverse translation correction, it increases to 23.4 dB. Finally, when the axial strain is accounted for as above described, the SNR reaches 26.8 dB. The 3D displacement for the 127 loading steps has been captured. Run in 300 s, this 4D procedure based on radiographs offers a fast measurement of the experiment with a gain of three orders of magnitude as compared to standard techniques. The entire procedure, including the experiment itself (after the initial tomography of the reference volume) and its analysis, is performed in approximately 10-15 minutes. A more complex and refined mesh could be used in order to characterize the kinematics more thoroughly. However, in order to avoid the increase of the number of degrees of freedom, regularization based on a mechanical model or other experimental observations would be welcome. A Projection-based Digital Image Correlation (P-DVC) method (presented in a companion paper [1]) is extended to an integrated approach for an elastoplastic law identification based on a single radiograph per loading step. Instead of following a two-step sequential procedure (i.e., first, measurement of the displacement field; second, identification), the integrated method aims at identifying few model parameters directly from the gray-level projections. The analysis of an in-situ tensile test composed of 127 loading steps performed in 6 minutes is presented. An isotropic elastoplastic constitutive law with free form hardening behavior (i.e., controlled by only 8 parameters) is identified and shows a very ductile behavior (up to 6.3 % strain before failure). A large improvement on the residual quality is shown and validates the proposed model and procedure. The obtained displacement fields revealed to be similar to those measured with no mechanical integration. A different parametrization of the constitutive law provides a very similar result, assessing the robustness of the procedure.

Introduction

The identification and validation of increasingly complex mechanical models is a major concern in experimental solid mechanics [2]. Model refinements require ever more precise and complex experimental protocols, from experimental setup and sample (e.g., complex multiaxial loading, optimized geometry) to sensor and acquisition devices (such as infrared, X-ray, or neutron imaging coupled with full field measurements). Especially in 4D situations (i.e., 3D space and time), the acquisition, processing and analysis of such a big amount of data (for the most part redundant) is challenging [3].

Computed Tomography (CT) allows materials to be imaged in 3D, thereby revealing their 3D microstructure in a non-destructive way [4,5,6]. Coupled with ex-situ or in-situ mechanical tests [7,8,9,10], it is possible to measure the deformation of a specimen with full field measurement techniques (Digital Volume Correlation, or DVC [11,12]) that provide a link with mechanical models and simulations. Intimate experiment-model dialogue has been developed to study crack propagation [13,14], plastic behavior [15], damage model [16], and allows for model identification and validation [17].

The measured displacement fields are used to identify model parameters with inverse methods (e.g., finite element model updating [18], virtual fields method [19]). Because model identification is the final goal of the entire experiment (contrary to displacement fields that are an intermediate step), one may identify and optimize model parameters directly on the acquired images. This fusion of two steps into one, has been called "integrated approach" [20,15] and leads to a drastic reduction of the number of unknowns. The more loading steps are acquired in the experiment, the more accurate and sensitive the identification procedure. DVC methods have hence been developed in a 4D space-time framework [15,21] using all available volumes globally. In lab-CT analysis, no more than 10 to 20 acquisitions are usually performed because of the acquisition time (approximately 1 hour is necessary to acquire the 1000 radiographs needed for one full 3D scan).

Projection-based DVC (P-DVC [22]) is a measurement method that aims to measure the 3D displacement fields directly from few radiographs instead of 3D volumes. The outputs of this method are the same as DVC procedures but require two to three orders of magnitude less data, and thus acquisition time. The comparison between two projections of the sample at the same angle with different loadings can be read as due to 3D motion. Previous works have developed this method for the analysis of a cracked cast iron sample imaged with a synchrotron X-ray source [23,24]. It was shown that the measurement was possible with only two orthogonal projection angles leading to a huge gain in acquisition time (of a factor 300). The identification of a brittle-elastic behavior of a plaster sample with crack propagation inside of a lab tomograph was performed within such a 4D framework [25].

In a companion paper [1], a measurement of the kinematics from a series of single projections per loading step was presented. Coupling P-DVC and model reduction methods provided a global measurement of the full 4D (i.e., space and time) displacement field. These results will be used as comparison means of the present approach. It is proposed in this paper to identify an elastoplastic model based on a fast 4D measurement. A standard approach would consist of performing the identification from the previously measured displacement fields, using for instance a Finite Element Model Updating (FEMU) procedure. In comparison with this pathway, it is shown that the fast 4D kinematic measurement and the identification can be performed in a unique "integrated" procedure, leading to more accurate determination of the constitutive parameters. A free form representation of the hardening behavior is chosen and controlled by 8 parameters. A variant with 7 parameters leads to a very similar hardening behavior showing the robustness of the procedure. The test case that was chosen in the companion paper for the kinematic measurement [1] is further used herein. A nodular graphite cast iron dog bone sample is tested in-situ with 127 loading steps performed in 6 minutes through uninterrupted loading until failure (leading to a gain in acquisition time of more than 2 orders of magnitude compared to standard methods).

In-situ mechanical test

The main features of the test are summarized here. The interested reader will find additional details in Ref. [1]. The tested material is a cast iron sample (similar to those studied in Refs. [26,15,17]). It is tested in uniaxial tension with a slender dog-bone geometry for the sample. An in-situ tensile testing machine similar to that of Buffière et al. [27] was used. The sample was scanned in a lab-tomograph (X-View X50-CT, North Star Imaging, 180 kV, 130 µA, W target). The voxel size (using 4 × 4 binning at the acquisition stage) was set to 10.7 µm. The radiographs of initial definition 954×768 pixels are cropped to a region of interest of 954×432 pixels. They are further coarsened by a factor of 2 in all directions (referred to as "scale 2") for speed up. Reconstructions and projections are performed with the ASTRA toolbox [28] with Feldkamp-Davis-Kress (FDK) procedure [29]. The initial projections f (r, θ) can be compared with the re-projection Π θ [F (x)](r) of the reconstructed volume F (x). The SNR of these projection residuals is 21.70 dB.

A 3D scan of the reference state (subjected to a small load of 250 N) is first acquired consisting of 600 radiographs captured at equally spaced angles ranging over a full 360°r evolution. The initial force is measured during the procedure (hence considered in the following stress measurement) but the initial elastic deformation of the sample, assumed to be very small are not considered in the displacement measurement. This scan is acquired in 22 min. Then, a continuous rotation of the sample is set with 50 acquisitions per turn. A total of 127 loading steps are acquired over 2.5 turns (see figure 4.1). A first series (50 projections, up to 100 s) is performed at constant load and is used to quantify the uncertainty. A second series (77 projections starting at 100 s) are carried out with a continuous load increase (from 250 to 750 N), controlled at a constant stroke velocity of 2 µm/s. 90) is representative of most of the test where the kinematics will be very well captured. The second one [START_REF] Suzuki | Projection-based motion estimation for cardiac functional analysis with high temporal resolution: a proof-of-concept study with digital phantom experiment[END_REF] will correspond to a more heterogeneous displacement field (strains not yet localized but concentrated in the central section), it is a more difficult case but will be seen to be very well captured. The last one (123), with a very similar load as the second, will be at the onset of localized strain (actually a crack has initiated). The axial force measurement κ during the procedure for each radiograph acquisition is shown in Figure 4.2. The herein proposed "integrated identification" is performed within the P-DVC framework [22,24]). The latter aims to evaluate the quantities of interest (be they kinematic or mechanical parameters) directly from the radiographs (i.e., the projection at an angle θ(t) of the 3D microstructure) acquired at different angles and loadings and not from the reconstructed volumes as in standard DVC methods [11,12,30].

The registration consists of finding the appropriate 3D displacement field, u(x, t) by minimizing the residual quadratic norm

χ 2 [u] = r,t ρ(r, t; u) 2 (4.1)
where ρ(r, t; u) is the residual field defined as the differences between N θ 2D projections g(r, t) of the deformed volume at different angles θ(t). The procedure makes use of the 3D reference image reconstructed using classical means, that provides for all voxels of the 3D space x = (x, y, z) the local x-ray absorption coefficient, F (x). This reference volume, corrected by the displacement, u(x, t), and projected at orientation θ(t) should coincide with the acquired projections when the solution is reached. In other words, introducing the reference volume deformed with any arbitrary trial displacement field, u,

F u (x) ≡ F (x -u(x, t)) , (4.2)
the residual field defined as

ρ(r, t; u) = Π θ(t) [ F u ] (r, t) -g(r, t) (4.3)
with Π θ(t) the projection operator in the θ(t) direction angle and r = (r, z) the coordinates in detector space.

The previous summand can be linearized considering small displacement field correction δu compared to the microstructure correlation length

χ 2 [u + δu] = r,t ρ(r, t; u) -Π θ(t) [δu • ∇ F u ] (r, t) 2 (4.4)
with ∇ the 3D gradient operator. It is noteworthy that after each evaluation of the displacement corrections δu from a known displacement u (n-1) such as u (n) = u (n-1) +δu, a correction of the volume F u (x) is performed so that the previous equation is used without approximation. The P-DVC framework requires the acquisition of one reference 3D volume in order to compute the correction term. The latter is to be performed ideally in the same conditions as the experiment itself, with or without load. This single (classical) tomographic scan is not challenging.

In order to validate the proper evaluation of the displacement, one should consider the magnitude of the residual field that highlights all projection differences that are not explained by the displacement field (e.g., noise, artifacts of the detector, ill convergence, model error). Ideally, it should be statistically indistinguishable from noise.

Beam regularization of the kinematics

Different regularization procedures of the displacement field have been introduced in the literature for global DIC methods where the kinematics is expressed on a finite element mesh [31]. Model-based regularizations, called strong regularizations or integrated methods with a reduced basis composed of elementary fields issued from a mechanical model [20,15] can be introduced. It is proposed, in the same spirit as global DVC, to express the displacement field in a reduced basis where the entire kinematics is controlled by the amplitude,

u i (t), of a few N f fields u(x, t) = N f i=1 u i (t)Φ i (x) (4.5)
The chosen reduced basis exploits the slender geometry of the sample that makes it amenable to a beam kinematics description [32]. Because of the slowly varying cross section, it is chosen to parametrize the displacement field by rigid body motions of N s equally spaced cross-sections, with a linear interpolation of the displacement in between. Because of the uniaxial tension imposed during the test, a single component of the motion (i.e., axial displacement along the z tensile axis). Hence the spatial displacements of the beam regularization, Φ i (x), is defined as

Φ i (x) = p i (z)e z (4.6)
where e z accounts for a uniform translation along the tensile axis and p i (z) the shape functions selecting a specific section and interpolating between them. If required by a poor quality of the final residual fields, additional degrees of freedom could easily be added.

It can be noted that no restriction holds on the number of control sections, nor on the type of axial interpolation, and hence the proposed discretization is very generic and is suited to any slender sample geometry.

Mechanical identification

Because the goal is to identify a mechanical behavior with an integrated method, the unknowns of the displacement field are written considering a constitutive model depending on N γ constitutive parameters, gathered in the column vector {γ}. Knowing the constitutive law, the local cross section and the loading at any time t (and possibly its complete history if needed as for elastoplasticity with unloading), it is straightforward to compute the mean axial displacement at control cross section, and its incremental variation for a variation δγ of the constitutive parameters about a reference set γ 0 for which the section displacements are u 0 i (t),

u i (t) = u 0 i (t) + Γ ij (t)δγ j (4.7)
where Γ ij denotes the kinematic sensitivity associated with each parameter γ j .

FEMU identification

Finite Element Updating Method (FEMU) [33,34,35] is an identification procedure based on the minimization of the quadratic difference between the experimental measured displacement field u exp i (t) for each beam section i (identified for each analyzed section [36]) and the computed field from the mechanical model u c i (t, γ)

χ 2 FEMU = i,t (u exp i (t) -u c i (t, γ)) 2 (4.8)
Here the most standard FEMU method is presented although a metric different from L2 (rather based on the covariance matrix of the measurement data) would be more appropriate [3].

The computed model has to be controlled by boundary conditions (generally the displacement field measured on the edges). In this identification approach, the measured load during the experiment is imposed and controls the kinematics of the beam. The minimization of the linearized functional with respect to the parameters γ leads to

δγ i = t Γ ij (t)Γ jk (t) -1 t Γ kl (t)(u exp l (t) -u c l (t, γ 0 )) (4.9)
This FEMU method will be compared with the results of the integrated approach based on the gray level images, and is considered as a baseline reference.

Integrated identification

The integrated approach aims at identifying the model parameters directly on the images hence minimizing the following functional

γ = Argmin γ (χ 2 u (x, T (t), {γ})) (4.10)
It can be noted that this cost function is similar to the one used in Ref. [1] for the P-DVC procedure, at the exception that the displacement field is now parameterized with the constitutive parameters directly. The derivative of this functional with respect to {γ} gives an equation very close to FEMU-U but weighted by the gray level uncertainty of each degree of freedom

{γ} = [H] -1 {h} (4.11)
where [H] is the "gray level" Hessian of χ 2 u with respect to {γ} and {h} the residual vector. It is convenient to write the projected gray level sensitivity for each degree of freedom of each section (i.e., the projected signature of a small motion for each degree of freedom)

S i (r, t) = Π θ(t) [Φ i (x) • ∇F (x)] (4.12)
and the associated matrix components

B ij (t) = r S i (r, t)S j (r, t) (4.13)
Computing and storing this sparse matrix is one of the longest operation in the procedure. However all the sensitivity fields for all parameters of the model are easily obtained from combinations of these projected 3D gray level sensitivity fields. It follows that

H ij = t Γ im (t)B mn (t)Γ jn (t) (4.14)
and

h j = r,t ρ(r, t)Γ jn (t)S n (r, t) (4.15)
Before applying this integrated approach, a first measurement and correction of the Rigid Body Motions (RBMs) is convenient. A 4D analysis based on a regularized space (including RBM) and time evolution is first performed using a Proper Generalized Decomposition approach [1]. Then a first estimate of the displacement field is based on scale 2 measurement, which allows large displacements to be coarsely captured (this scale is particularly focused on the edge alignment of the specimen). When scale 2 has converged, the measured displacement field is doubled and used to initialize scale 1 computations.

Application

Elastoplastic behavior

This paper aims at identifying an isotropic elastoplastic constitutive law. A "free form" isotropic hardening is proposed with N γ = 8 degrees of freedom to describe a piece-wise linear stress strain curve. Such discrete model can be the input behavior for industrial FE softwares. Contrary to some parametric models (e.g., Ludwik [37], Voce [38] or more complex models [39]) whose parameters may lack of sensitivity in the measurement procedure, the chosen free form behavior can be designed for identification according to the analyzed experiment.

The degrees of freedom are chosen to have a fixed stress level (imposed by the applied load and sample geometry) while the corresponding strain is unknown. The parameters and initial values are reported in Table 4.1. The last parameter is chosen as the ultimate strength of all sections during the test. Therefore it controls all sections with high strain levels (i.e., central part of the sample at the highest load). A cubic interpolation, in the [σ, ε] space is performed between the control points. Other interpolations have been tested in the procedure (e.g., linear) but show no appreciable difference. Because of a low sensitivity on the elastic part, the behavior is focused on the plastic part (i.e., the first slope of the stress-strain curve will not be interpreted as the Young modulus).

In addition to the 8 unknown parameters that control the elastoplastic behavior, a vertical translation is introduced as an additional degree of freedom that may vary arbitrarily in time. Although the vertical motion had been corrected (with respect to gray level residuals) prior to the identification procedure, this extra freedom allows a fine tuning of the vertical position without any interference with mechanical identification.

In order to study the robustness of the identified plastic law with respect to the discretization, another identification based on 7 parameters, at different stress levels, is also performed. The same overall behavior is expected.

Initially corrected residual fields

An initial measurement and correction of all rigid body motions is performed based on the method described in Ref. [1]. The SNR increases very significantly from 9.9 to 25.4. The residual field for the three selected angles is shown in Figure 4.3 with a diverging color map to highlight positive and negative patterns. With this first procedure, the mean displacement field is corrected and hence the motion appears to be well corrected in the center of the specimen, i.e., the top (resp. bottom) part is moving upward (resp. downward). Alternating positive and negative spots are visible on the top and bottom parts, which are signature of a vertical motion of the sample. Those residuals will be corrected by the following procedure. Note that the color bar amplitudes differ in the figures as they are adapted to the residual amplitude.

Results of the integrated approach

From the initial guess of parameters (shown in Table 4.1), the procedure aims to identify the set of parameters that leads to the residual field having the lowest norm. The procedure converges in 5 to 7 iterations. A new set of parameters is found (Table 4.2). The values of the parameters have significantly changed with respect to their initial guess (Table 4.1). This result shows that even though the initialization was very far from the converged solution, it did not prevent the algorithm to find a good solution (i.e., the gray level residuals were significantly reduced). This new set of parameters leads to a mean value of the SNR of 26.4 (the SNR for each radiograph is shown in Figure 4.4, curve (b)). A decrease of the residuals (i.e., increase of the SNR) is observed from step 80 on, where the displacement correction is large because of the plastic yield. The SNR decreases at the end of the loading (120) step where the model may not be adapted (i.e., high strains including localized features such as necking and shear bands). Low SNRs are located every 25 time steps and correspond to angles where the sample edges are "aligned" with the X-ray beam hence where the residual field is very sensitive to radial positioning. A very small mis-positioning of the sample generates large residuals. A small SNR at step 50 is due to the first load increase that may have generated spurious vibrations. and Integrated P-DVC with the 7 parameter variant). The positions of the control points in stress were selected at will, which allowed the hardening curve to be refined where needed. The three approaches lead to very similar shapes especially for the higher strain part (that also corresponds to a high displacement sensitivity after integration of the strain). The difference between 7 and 8 control points (blue and black curves) is minimal, which means that the procedure is not much sensitive to the chosen discretizations. These three results are very different from the initial guess, which proves the robustness of the proposed framework. The stress-strain curve shows a ductile behavior with a total strain reaching approximately 6.3 % in the central part. The longitudinal space-time displacement field obtained from the identified model is shown in The residual fields at the end of the analysis are presented in Figure 4.6. Most of the previous patterns have been erased. Some residuals are still visible in the bottom part of the sample, especially at the end of the experiment and may be due to the localization of plastic strain and/or damage due to a softening behavior. This interpretation is supported by a thin transverse mark at mid-height of the sample (z ≈ 500 pixels) at the end of the experiment that is visible in Figure 4.6(c) and may be interpreted as damage localization. This is also consistent with the fact that failure took place at the next scheduled radiograph. The onset of this localized strain band could be traced back to time step 115 where a degradation of the SNR could be detected. Although such a phenomenon may occur in reality, as it leads to a breakdown of uniqueness in the solution [40,41,42,43], this situation was excluded from the scope of identification, and the absence of softening was enforced in the parametrization of the constitutive law.

Comparison with kinematic measurements

A comparison with the results obtained in Ref. [1] can be performed. A first point is that the SNR whose mean value is 26.4 for the integrated approach and 26.3 with FEMU, is barely smaller than that obtained from the kinematic measurement, namely, 26.7. The temporal changes of the SNR with the two approaches are shown in Figure 4.4, curves (bc). Even though it is obtained with significantly larger amount of degrees of freedom the kinematic approach has only slightly lower residuals. The accuracy of the plastic model based on no more than 7 or 8 parameters is remarkable.

A second comparison has already been performed with the FEMU procedure using the measured displacement field of Ref. [1]. The stress-strain curve is extremely close to the integrated identification (the first two degrees of freedom are different but they correspond to low displacement sensitivities). The two approaches lead to results that can hardly be distinguished.

The relevance of the chosen model can also be judged from the difference in the identified displacement field. The latter is reported in Figure 4.7. The mean of the absolute values of the difference is 0.37 voxel and the standard deviation is 0.56 voxel. The low spatial frequency differences are due to model differences (i.e., the identified integrated displacement cannot reproduce the localized behavior). , obtained from the difference between the displacement field obtained from the integrated approach and that measured with the kinematic approach of Ref. [1] The comparison of the two displacement fields is presented here to show a visual appreciation of the correspondence. It cannot be considered as a judgment criterion of the methods as none of the two fields can be considered as a reference since they are built on different kinematic bases. The only objective comparison metric is given by the gray level residual as discussed previously [17,3].

Conclusion and Perspectives

An integrated identification method is developed within the Projection Based Digital Volume Correlation framework. Relying on a mechanical model with 8 unknown parameters, the identification is carried out from radiographs (i.e., projections of the 3D microstructure) instead of 3D (reconstructed) volumes. As in the previous purely kinematic study [1], this approach enables huge gains in acquisition time (about 2.5 orders of magnitude as compared to standard approaches). Additionally, rather than providing a space-time displacement field, the output of the integrated approach is the numerical estimate of the parameters of the constitutive law, together with a mechanically admissible displacement field.

An illustration is shown with an in-situ tensile test on a cast iron sample. The acquisition is performed on the fly (i.e., with continuous rotation of the sample inside of the tomograph and a continuous load history until failure). The 127 radiographs at different loading steps are acquired in 300 s. The aim of this test is to identify a 1D elastoplastic model of a sample treated as a beam. The (free form) constitutive law is controlled with 8 strain parameters associated to fixed stress values.

An initial analysis of the test is devoted to the measurement (and correction) of the entire set of rigid body motions. Then the identification procedure shows an important decrease of the residual norm. A large part of the residual fields is erased with the kinematic correction, meaning that motions have been well captured by the model and this validates the proposed constitutive law. The onset of strain localization can be guessed from the late residuals showing the benefit of using residuals to evaluate the success and limitation of the identification procedure.

One requirement of P-DVC is that the 2D contrast obtained from the projections of the imaged volume can be safely exploited. The gradient of this radiograph weighs the sensitivity fields and is thus a key element of the identification procedure. In the present application, nodular graphite cast iron is a very good candidate for P-DVC as it is composed of graphite nodules (whose X-ray absorption is much lower than that of iron and provides the needed contrast). The application of the method to weakly contrasted samples is a challenging issue for future applications.

Acknowledgement

This work has benefited from the support of the French "Agence Nationale de la Recherche" through the "Investissements d'avenir" Program under the reference ANR-10-EQPX-37 MATMECA and ANR-14-CE07-0034-02 COMINSIDE project.

Chapter 5 Water flow in a sandstone sample is studied during an experiment in-situ in a neutron tomography set-up. In this paper, a projection-based methodology for fast tracking of the imbibition front in 3D is presented. The procedure exploits each individual neutron 2D radiograph, instead of the tomographic-reconstructed images, to identify the 4D (space and time) saturation field, offering a much higher time resolution than more standard reconstruction based methods. Based on strong space and time regularizations of the fluid flow, with an a priori defined space and time shape functions, the front shape is identified at each projection time step. This procedure aiming at a fast tracking the fluid advance is explored through two examples. The first one shows that the fluid motion that occurs during one single 180°scan can be resolved at 5 Hz with a sub-pixel accuracy whereas it cannot be unraveled with plain tomographic reconstruction. The second example is composed of 42 radiographs acquired all along a complete fluid invasion in the sample. This experiment uses the very same approach with the additional difficulty of large fluid displacement in between two projections. As compared to the classical approach based on full reconstructions at each invasion stage, the proposed methodology in the studied examples is roughly 300 times faster offering an enhanced time resolution.

Fluid invasion tracking with neutron tomography

Introduction

Studies of fluid flow and transport in porous rock have received increasing interest related to a number of different resource engineering challenges, e.g., hydrocarbon production and CO 2 sequestration, amongst other applications. Standard laboratory approaches to study such phenomena are based on "global" flow evolution (e.g., total sample saturation vs. time), which only provide spatial average measurements over the whole specimen. Local measurements of the saturation field during a fluid invasion process would however be highly desirable to validate and calibrate pressure driven flow and imbibition models, especially in samples where a heterogeneous microstructure affects porosity, permeability and wetting properties in an intimately coupled way, meaning that it is difficult to disentangle the influence of each one. Heterogeneity is particularly important, for example, in stratified media or mechanically deformed samples that may contain cracks or shear bands. In such cases, the fluid flow process is intrinsically a 4D phenomenon that requires dedicated 4D measurement methods [1,2]. Being able to image this fast time dependent fluid advance at high speed as in real cases is a challenge in the direction of the study of pressure driven flow.

The recent developments of fast synchrotron tomography [3] allow one to image 3D evolutions of fluid invasion inside samples at a fast rate (20 Hz in [4,5]). Both X-ray [6,7,8,9] and neutron [10,11,12] imaging are used to image absorbing fluid (e.g., zinc iodide for X-rays [13,14,5] or water for neutrons). These two tomographic modalities can even be combined to benefit from their different contrast sensitivities as shown by Hall [15].

Many methods have been developed to extract the position of the fluid from 3D volumes or slices. From simple segmentation of the gray level images to more sophisticated method based on model driven evolution such as history matching [16,13,17]. History matching is an iterative inverse problem, commonly employed in petroleum reservoir engineering, whereby a reservoir simulation model is tuned to be consistent with experiment and accurately describes the fluid advance. However, since these identification methods are based on 3D volumes, they require quasi-static flows (although this constraint becomes less stringent with the progress towards increasingly fast 3D image acquisition). Fluid flow during a tomographic scan leads to blurring in the reconstructed volume around the fluidfront. Eventually, some acquired scans can not be reconstructed because of the discrepancy between the initial and final saturation states over the scan duration. One way to circumvent such difficulties is to resort to samples prepared for being invariant in one direction, and where the fluid advance is expected to be read easily in the radiographs [18,19,20]. Armstronget al.proposed in [21] a method based on radiograph differences to extract the fluid advance with a higher time resolution.

Similar temporal issues to those described above can be found during mechanical testing in-situ with tomography where the sample may evolve (due to viscous effects, creep, relaxation, thermal variations, etc.) during the tomographic acquisition. These effects lead to blurring that affects quantitative measurements based on the image contrast (i.e., full-field space-time measurement [22]). To address this problem, recent methods called Projection-based Digital Volume Correlation (P-DVC) [23,24] and Discrete Digital Pro-jections Correlation [25] have been developed for the measurement of 3D displacement fields based on radiographs rather than volumes. The P-DVC method, after a first (classical) tomography performed before the mechanical test, does not require more than two orthogonal radiographs per loading step of the mechanical test to identify the full 3D displacement field. This method, enriched with a model-based regularization, and a model reduction strategy to extract principal modes, makes fast-rate 4D imaging -even in laboratory tomographs -possible [26].

The principle of the P-DVC method is, at first order (i.e., unless iterative methods are used), the correction of a reference 3D volume by the weighted gradient of this volume so that its projections match with radiographs of the deformed volume. The method proposed in this paper, being a projection matching optimization, is in a similar vein. The final shape of a fluid volume is obtained from an acquisition of the saturated sample. This volume is then corrected iteratively to match, in the projected space, with radiographs at all invasion states. This optimization method is based on the minimization of the quadratic difference between the projected model of the fluid column and the real fluid advance.

The resolution of space-time evolution can be carried out with a separation of space description and time evolution using Proper Generalized Decomposition (PGD) techniques [27,28]. This method consists of successive enrichments of a field adding a new space-time direction at each iteration, each term of the sum being sought a priori in separate representations in space and time. This approach captures the most important modes first and allows one to stop when enough fine details have been incorporated. The spacetime decomposition for measurement has been studied with PGD-DIC and PGD-DVC developed in [29,30] with one dimensional space functions.

The aim of this study is to capture the 4D evolution of fluid invasion (4D space-time evolution of the saturation front) at each neutron radiography time, that is at a rate of 5 fps, for the particular example of the study. First, a description of the method based on a space-time separation is proposed. This is followed by applications of the method to two test cases, which are extracted from an experiment that was initially designed for serial 3D measurements. This experiment involved imaging of water flow into a sandstone specimen by neutron tomography. Because the studied porous sandstone is reasonably homogeneous at the scale of the observations, the fluid flow can be pictured as a rather steep variation from a fully saturated lower part of the sample to a dry upper part. Hence, the fluid (or saturation) front is modeled as a moving 2D surface about which a saturation gradient may be added as a further refinement. The first test aims at tracking the flow at the radiography rate, i.e., 5 Hz, during a single 3D scan, while the sample was being rotated. In this case, the front traveled by a short distance of about 8 voxels. This evolution cannot be resolved by a classical reconstruction procedure with the assumption that the front remains perfectly still during the full scan duration. The second test case addresses the entire invasion process, where the front sweeps through the entire sample height. The saturation profile and velocity field are described and compared with a front detection method from 3D reconstructed volumes.

Method

Projection based method

Tomography, be it obtained with X-rays or neutrons, consists of recording (twodimensional) radiographs ρ(r, t) where ρ is the intensity, r the position in the detector, acquired at the rotation angle β(t) of the sample about the z axis at the time of capture t.

To reconstruct a 3D volume from the radiographs, it is necessary to use a so-called "flat field", ρ 0 (r), that corresponds to the beam intensity at the same position in the detector. According to standard reconstruction procedure, the ratio ρ/ρ 0 is the dimensionless beam attenuation, whose logarithm p(r, t) = log(ρ(r, t)) -log(ρ 0 (r)) is the line sum of local absorption coefficient along the ray hitting the detector at position r. In the following p is called a "projection". In the present case of neutron tomography, the beam is approximated as being parallel and the volume A(x), (where x = (x, y, z) designates the 3D coordinates) can be obtained from an inverse Radon transform of the projections, p(r, t).

For different reasons, the simple Beer-Lambert law used above to relate the length of neutron ray through the material to the log of the attenuation may be corrected by non-linear (so called "beam hardening", BH) corrections [31]. This can be seen from fully homogeneously saturated or dried sections of the sample that are not reconstructed with a homogeneous gray level. This effect can however be exploited to compute the projection correction. A first reconstruction and segmentation of the volume allows this solid and fluid length to be measured. At this stage, the collected data should allow for the determination of the intensity reduction after a length traveled through the pure solid phase and in the mixed dry-fluid phase. It can be noticed that the dry part has a very low absorption and does not exhibit significant beam hardening effects compare with the fluid-filled part. Hence an intensity correction of the high absorption would correct the beam hardening of the fluid part without changing the dry part. A simple beam hardening correction proposed by Herman et al. [32] where the linear projection is used as part of a higher order polynomial projection has been used in this work. The correction of the projections as a function of a 5th order polynomial allows the beam hardening of the fluid-invaded part to be removed without changing the intensity of the dry part.

The proposed approach for the identification of the fluid advance in a porous sample is inspired by the Projection-based Digital Volume Correlation (P-DVC) technique [23]. It is important to note that the challenge here, in the fluid-tracking, is to achieve a 3D description of the advancing fluid front at the same temporal rate as the 2D projection acquisition, in contrast to the usual practice of performing a 3D reconstruction assuming a still state during the acquisition.

The starting point of the method is a reference ("standard") tomography scan, performed at the beginning of the experiment and, ideally, with a dry sample. This reference is used to extract a 3D image composed only of the shape of the projected fluid (as in [20]). (Note however that the initial state in the experiment presented later was not perfectly dry and a small layer of fluid was present at the bottom.) The aim is thus to determine a description of a scalar saturation field S(x, t) in 3D space.

Introducing the porosity φ(x) of the porous medium, the projection can be written

p(r, t) -p 0 (r, t) = Π β(t) [φ(x)S(x, t)] (5.1)
where Π β(t) is the projection operator at angle β(t). This decomposition is illustrated in figure 5.1. The rotation axis of the sample, denoted as z, is perpendicular to the neutron beam, and is one principal axis of the detector. A gray level color coding has been chosen for radiographs and a diverging color palette (where 0 is encoded as white, while positive and negative values have a different dominant color) is used for all residual images to facilitate reading. Based on the above description of the projections for a given saturation state, the determination of the 4D saturation field is performed through the minimization of a cost function that is chosen as the quadratic norm, χ 2 , of the residuals between the predicted and the actual projection data η(r, t) = p(r, t) -p 0 (r, t) -Π β(t) [φ(x)S(x, t)]

(5.2)

thus χ 2 = r,t η(r, t) 2 (5.3)
The minimized quantity is the norm of the residual field that gives a precious information on the quality of the solution because it highlights all patterns that have not been correctly captured (because of noise, artifacts, model error, ill-convergence, etc.) and it is based on the raw acquired data, with a minimal influence of prior assumptions that could bias the judgment.

Fluid field model and parametrization

Different descriptions of fluid flow in porous materials can be found in the literature, from detailed models accounting for presence or not of the fluid in the pore space, or more homogenized versions where the fraction of the pore space invaded by the fluid, i.e., the saturation S, is used. The choice of the description used is essential, as it enables tuning of the number of needed parameters, which has a drastic impact on our ability to invert the problem. The ideal choice is to restrict the saturation description to the smallest possible number of parameters yet maintaining a sufficient representation of the reality.

To have an accurate model for the projections, it is necessary to know the amount of fluid intersected by the neutron beam, which is given by the product of the saturation S(x, t) and the (connected) porosity φ(x) (independent of time).

The saturation itself can be determined by exploiting the following assumptions: -First, the fluid "front" is assumed to be an abrupt transition in the saturation field where it drops to 0 in one pixel distance. It can be described by a single valued surface, z = h(x, y, t), and hence S = 0 for z > h(x, y, t). This front moves upward as invasion takes place -Second, behind the front, z < h(x, y, t), saturation increases progressively in time and finally reaches a full saturation. The saturation profile seems to be well described by a function of the distance to the front h(x, y, t) -z, this property is not a universal property of pressure-driven flow, but it is observed to be well obeyed in the studied case.

The surface of the fluid front can be expressed, with a separation of variables, as a sum of modes, composed by the products of a spatial field and a temporal function, respectively F m and τ m , h(x, y, t) = Nm m=1 τ m (t)F m (x, y).

(5.4)

From this general expression of the front, (that is not restrictive) the spatial field can be specialized to a reduced space generated by N s shape functions, ϕ i (x, y) with i = 1, ..., N s , weighted by amplitudes f m i . Different levels of description are possible to regularize the space front. One is "free-form" from a fine finite-element mesh description of the surface. This description requires many degrees of freedom. A local "smoothing" can be added to constrain all independent nodal degrees of freedom, and obtain a better conditioned problem. Alternatively, a "strong" regularization can be defined extracting the saturation parametrization from a pressure driven flow model. The latter is generically the best choice, as it leads to very few degrees of freedom, although it does require a physicallyrelevant model of the saturation consistent with the scale of observation.

It is desirable to reduce the number of degrees of freedom, to couple them together for better conditioning and to ease dialogue with models. The fluid front can thus be written as

h(x, y, t) = Nm m=1 τ m (t) Ns i=1 f m i ϕ i (x, y).
(5.5)

The space shape functions, f m i , can be tuned to adjust their spatial resolution: from voxels (Dirac) representation, to finite element shape function for a 2D mesh of the cross section, or to polynomials of low order, for example. In the same spirit, the time evolution can be decomposed over temporal shape functions, θ l (t),

h(x, y, t) = Nm m=1 Nt l=1 g m l θ l (t) Ns i=1 f m i ϕ i (x, y) . (5.6)
It can be noted that only the norm of the product g m f m matters so that an additional convention such as: max Ns i=1 f m i ϕ i (x, y) = 1, has to be chosen without consequence.

Parametrization of bulk saturation

Detailed study of the images of the invasion process indicates that (at least in the regime relevant for our experiment) the saturation shows a rapid first increase from 0 to an intermediate value followed by a slower second stage increase up to 1. This secondary stage can be seen as being ruled by a unique 1D function of the vertical distance to the front. Following this observation, the saturation can be written as

S(x, t) = ζ p (h(x, y, t) -z), (5.7)
where ζ p is a smooth function of its argument, null for negative values, providing the abrupt cut-off at the front. This function ζ p can be parametrized by a set of parameters, collectively denoted as a vector p.

As a final step, to soften the very steep vertical variation occurring at the front, the above saturation field is convoluted with a Gaussian along the z direction. An additional component in the vector p will refer to the width of this Gaussian. Since ζ is a slowly varying function in z except at the front where it is discontinuous, the effect of this last smoothing operation only affects the immediate vicinity of the front and has no perceptible influence on the bulk saturation.

Inverse problem resolution

The previous section described how the saturation field is described by the front (with parameters gathered in two vectors f and g) and by the bulk saturation fields with parameters p. In this section, the solution to the inverse problem is described.

Although the best was done to introduce the fewest possible number of unknown parameters to describe the fluid advance, their number is still too large for the problem to be tractable. The PGD framework [27,28] to progressive enrichment of the model description with modes is both a very powerful and convenient approach. At each iteration of the algorithm, the fluid front can be written as the sum of previously identified modes plus a new one. Hence,

h (m) (x, y, t) = h (m-1) (x, y, t) + Nt l=1 g m l θ l (t) Ns i=1
f m i ϕ i (x, y) .

(5.8)

The minimization of the functional leads to the determination of the unknowns from 3 coupled equations,

f m = Argmin f (χ(f , g, p) 2 ), (5.9)
g m = Argmin g (χ(f , g, p) 2 ), (5.10) p = Argmin p (χ(f , g, p) 2 ).

(5.11)

To solve this coupled system, a fixed point procedure is used, solving for the above three minimization problems successively, keeping the other two sets of parameters fixed, as described in Algorithm 3. It was observed that a very low number (3-5) of iterations was needed to reach a stable solution, and an ill-convergence escape condition based on a maximum number of iterations was never reached. Even if the two space and time unknown amplitudes are defined by the norm of their products, the fixed point algorithm, composed of separated equations leads to a correct determination (from the choice of the initialization of g m from which f m is computed) allowing convergence. The choice of a norm affecting each unknown can finally be set (or reset) at the end without consequence. The three parts of the fixed points are obtained from the above functional using a Newton algorithm. The derivative with respect to h leads to

f = N -1 n, (5.12)
with N the spatial Hessian matrix of χ 2 with respect to f (i.e., N ij = ∂ f i ∂ f j χ 2 ), built from the sensitivity field, Ψ m i (r, t), i.e., the variation of the projection for an infinitesimal variation of input parameter,

Ψ m i (r, t) = ∂Π β(t) [φ(x)S(x, t)] ∂f m i , (5.13)
n is the second member vector based on the residual field, which can be written

N ij = t,r Ψ m i (r, t)Ψ m j (r, t), (5.14
)

n i = t,r
η(r, t)Ψ m i (r, t).

(5.15)

The exact same formulation is found for the derivative with respect to g and p. It is noteworthy that p does not depend of m in the fixed point procedure. It is hence an update of this parameter that is performed at each mode identification.

Applications 5.3.1 Presentation of the test

The studied sample is made of a porous sandstone from the Vosges mountains in the east of France [33]. The sample shape is cylindrical (diameter 38 mm and height 80 mm) with two opposing flattened surfaces on which two shallow notches were cut for a related study on the formation of localized deformation under triaxial compression (see [15,34]). The sample geometry is shown in Figure 5.2(a). The fluid advance in the sample was imaged with neutron tomography at Helmholtz-Zentrum facility in Berlin (HZB) [35]. For the experiments, the sample was wrapped with thin PTFE tape then two Teflon inserts were placed on the flattened surfaces to make the overall shape cylindrical and a heat-shrink FEP membrane, which is reasonably transparent to neutrons, was shrunk to fit tightly around the sample and inserts. The PTFE tape was placed between the sample and the inserts to avoid preferential flow paths. The dry sample was placed and sealed with silicon in a cup that could be fed, by a tube from the bottom, with distilled water. The experimental setup, described in details in [36], is shown figure 5.2(b). Water was supplied to the cup at the base of the sample and the fluid advance was captured with 137 neutron tomography scans over a period of 2 hours.

The pressure of the water was given from the height of the water tank, which was adjusted during the experiment to compensate the water loss and keep the pressure constant. The fluid advance in the sample was driven more by the water pressure head than by capillarity.

Neutron tomographies were acquired with an acquisition time of 0.2 sec per radiograph and 307 projections over 180°, for a total scan time of 1 min. The rotation of the table was continuous with a fixed velocity. 130 flat-fields and dark-fields were acquired and used to perform a standard flat-field correction. The images are 541×648 pixels 2 with a pixel size of 110 µm. For a reduced computation time and because the front did not reveal any high spatial frequencies, the images have been down-sampled with a factor of 4 (i.e., after a Gaussian filtering of the images with a characteristic length of 4 voxel size, the smaller images are created sampling 1 pixel out of 4 in all space directions). Hence the images treated in the following are 138×162 pixels 2 with a pixel size of 440 µm.

To estimate the actual porosity field of the sample, φ(x), the final state after steady fluid invasion was reconstructed from a complete 180°scan reconstruction using an algebraic SIRT procedure from the ASTRA Toolbox [37] similar to the initial state reconstruction. Because the first scan was performed with a small quantity of fluid in the bottom of the sample, a mask was applied to the bottom part of the projections. The following surface height is defined from z 0 = 20 voxels.

Choice of parameterization

Spatial parameters

The saturation front is described using a very simple model: a sum of low order bivariate polynomials of order 3, thus 10 degrees of freedom for the surface function are allowed for each mode.

Temporal parameters

A single projection cannot be used to identify all modes and indeterminacies would result if each single time step (and hence projection angle) were treated independently. This can be resolved by describing the time evolution by time shape functions -denoted as θ l (t) in Eq. 5.6 -that span several time steps or angles; these also introduce an additional smoothness. Again, the limit is that this smoothening should be consistent with reality and should not enforce an oversimplified definition of the front.

The time evolution is chosen to be expressed on a reduced basis composed of low order polynomials. Because the first case example has a simple and low frequency fluid advance, the time regularization was set using polynomials of order up to 2. The second test case is more complex and was regularized in time with polynomials up to 4th order. The proposed regularizations can be decided to be valid, or not, based on the residual field at the end of the identification procedure.

Bulk saturation parameters

The bulk saturation of the porous medium behind the invasion front σ(x, t, p) = ζ(h(x, t)z) is described with two parameters, p 1 and p 2 , and an exponential function.

ζ(δh) =    0 δh < 0 1 -(1 -p 1 ) exp - δh p 2 else (5.16)
The saturation jump at the front is considered as a discontinuity from 0 to p 1 . Although abrupt, the saturation evolution at the front is better described if the transition is smeared over a small distance by a convolution of ζ with a vertical Gaussian kernel of width 2 voxels. Note that this smoothing can be performed in the detector space for a faster computation.

The plot of the identified ζ(δh) at the end of the procedure can be used to validate the chosen intensity model (in the example, p 1 = 0.64 and p 2 = 7.7 vox); for the current example this is shown in Figure 5.3. The chosen and identified profile, especially after the Gaussian smoothing are close to the profile extracted from the real data. Note that the reconstructed 3D volume may have been blurred because of fluid motion during the acquisition.

Case 1: during a scan rotation

The first test case involves the measurement of the fluid advance from the radiographs acquired during a standard tomographic scan. The considered data were extracted from the middle of entire experiment and involve 307 radiographs that were acquired at 5 Hz frequency over 60 seconds during a continuous rotation of the sample from β = 0°to 180°. The identification procedure consists of 3 modes, each composed of 2 unknowns for the saturation profile (vector p), 10 unknowns in the spatial domain f m i and 307 (resp. 3) unknowns in the time domain, g m l , for a free identification (resp. for a 2nd order polynomial).

Identification of the first mode for the fluid front

The first spatial and time evolution is shown in figure 5.5. The range of variation of the first spatial mode, about 0.1, is small relative to its mean value (about 0.9), and hence the time evolution, to first order, reveals an invasion of the fluid into the sample that is initially fast and then more slowly saturates the region of the sample by the end of the 60 s acquisition. The variation of the mean front position is a displacement of approximately 2 voxels at scale 4 (hence 8 voxels at scale 1, for a full resolution reconstruction). It is possible to extract the amplitude at each single time step by taking into account only this first spatial field and to use it to match projections. This provides a free-form estimate of amplitudes f 1 i for local time shape functions, f 1 i (t) = δ(t -t i ). Such amplitudes are shown in red in Figure 5.5b. If instead, the time shape functions are chosen as polynomials of order 2, the resulting determination of amplitudes is shown as the black curve. A fair agreement is observed, in spite of the drastic reduction in the number of degrees of freedom (from 307 down to 3). Moreover, this reduction comes with a much more robust (noise resilient) behavior. In addition to describing the mean volume invasion, the non-trivial front shape also captures a persistent non-uniformity of the front. However, higher mode orders are expected to contribute significantly to a finer description of the front shape.

The true metric to evaluate the merit of our procedure is the residual field η(r, t) which shows everything that has not been captured by the procedure (noise, artifact, model error etc.). To evaluate the relevance of the non-trivial front shape of the first mode, it is easy to perform the same analysis with only one mode, but restricting the spatial part of F , to a single degree of freedom, namely a constant height (polynomial order equal to 0). An arbitrary choice of the amplitude of the space and time function has to be chosen. It was decided to set the maximum space amplitude to 1 to help for a visual understanding of the surface evolution. At convergence, the final residual, is shown in figure 5.6 at the beginning and end of the scan. The most critical part of this residual field that has to be scrutinized is the fluid surface zone (interface between empty and saturated pore volume). It can be seen that the mean evolution is well captured, but alternate positive and negative values show that this description is still inaccurate. The value of χ 2 reaches 144 for a flat front fixed in time, and drops down to 103 when a variation in time is considered, showing that the latter choice allows one to improve significantly the solution. Let us note that with a single mode, when higher polynomial orders are considered in space, the residual decreases only slightly to 102. When more modes are added, the spatial enrichment reveals much more rewarding as detailed in the following.

Identification of additional modes

A comparison of the final residual field considering different time regularization is proposed in table 5.1. From these results it can be seen that, a third order polynomial for the time evolution is required to correctly capture the fluid advance.

Most of the residuals are erased with 3 modes, which indicates that a solution for the front shape has been found. The final residual field is shown in figure 5.7. Some structure can still be distinguished in the residuals below the front, which likely corresponds to reconstruction artifacts from the projected final fluid column.

The identification procedure was been stopped after the identification of 4 modes because the norm of the residual fields indicated convergence. The evolution of the norm of the residual fields is presented table 5.1. It can be seen that the residual field does not decrease further after mode 3, indicating that the main trends have been captured. The space time evolution (with a product expressed in voxels) of the different modes is shown in figure 5.8.

The final surface can be reconstructed from the three identified modes. In figure 5.9, the front surface at two time steps are plotted. The evolution from time t = 0 s and t = 60 s is shown with the same color coding. It can be noted that the surface curvature changes progressively during invasion over the scan time.

Comparison with 3D front extraction method

Using classical tomographic approaches, the beam hardening corrected projections can be used to reconstruct a 3D volume (assuming implicitly a steady fluid), from which a segmentation method, based on the value where the front saturation drops to 0 can be used to extract the shape and position of the saturation front. The measurement of the surface height can thus be determined. Obviously, this method allows a single front to be extracted for the time interval, which should be, in some way, an intermediate mean position. If the 3D front is determined by standard segmentation procedures, its height will be given with integer voxel resolution. Figure 5.10 shows the difference between the interface extracted from the tomography and the time-averaged interface as identified with the method presented in this paper. In this figure few steps are visible corresponding to jumps due to voxel rounding. The difference is always less than 1 voxel, corresponding to the 3D segmented uncertainty. The standard deviation of the difference is 0.36 voxel showing a good agreement in those two results. However the time evolution cannot be obtained with classical reconstruction procedures.

Figure 5.10: Difference, in voxels, of the surface height between the mean identified surface with the proposed approach and the surface extracted from the segmented (voxelized) 3D volume.

Case 2: during the entire test

The above described identification procedure is now applied to identify the fluid front over the entire pressure-driven flow experiment. The analyzed test is composed of 42 180°-scan series. To simulate the acquisition of a fast 4D procedure, only a single projection per invasion state has been selected from each scan in the series for a projection angle that evolves constantly with time just as if the experiment would have been conducted in 4 full turns of the sample (approximately 10 projections per turn). Projections at 3 times

Conclusion and discussion

In order to follow a fast invasion process (water invading a porous sandstone), a new methodology has been introduced, based on a time regularization and a model reduction technique. Such a procedure allows to resolve high rate invasion (up to 5 Hz) that cannot be captured using classical tomographic acquisitions and reconstructions. The method relies on a space and time regularization through a choice of a suited basis, and is based on a projection-based identification of the unknown parameters.

To further reduce the sought parameters, a Proper Generalized Decomposition has been designed, to capture "modes" one at a time until the residuals have been reduced down to an acceptable level. The residual field at the end of the procedure permits a validation of the chosen basis.

A first case example, addressed the question of a slow evolution during a 5 Hz scan. The residual field, after convergence of the identification procedure, does not reveal any salient 5.5. ACKNOWLEDGMENT 123 unresolved features at the fluid front, which indicates that the front was well captured by the chosen description and the regularization is clearly well suited.

A second test case aimed at analyzing a fast pressure-driven flow experiment, and 40 radiographs acquired during about 4 full revolution of the sample were extracted from a longer time series. In contrast to the previous example, large displacements were observed and yet the proposed P-DVC algorithm could deal with such large changes.

The example applications of the method were performed at scale 4 (i.e., the volume was reduced in size by a factor of 4), so only small data volumes were considered. For the two test cases, because the surface was smooth, this resolution was enough to capture the fluid shape. The measurement at full resolution would be in the same spirit, with a higher computation time. With the actual procedure, each mode of case 1 is obtained in 55-65 minutes and each mode of case 2 in 4-6 minutes depending on the convergence rate of the fixed point procedure.

The simplicity/complexity of the problem (i.e., the number of degrees of freedom, in space and time) is conditioned by the scientific background of the experimentalist about the observed phenomenon. For example, the knowledge of fluid flow models composed of few parameters allows the problem to be easily solved with a small number of unknowns. In the opposite, no expectation about the experiment requires to leave much more freedom to the system (too many would lead to ill-conditioning). At convergence, the residual fields enable the quality of the parametrization and chosen model to be appreciated. It is then possible to enrich the studied basis (i.e., add an additional curvature of the front, space or time discontinuities, evolution of the front width etc.).

The method could be applied to more complex cases such as imbibition or pressuredriven flow of water, or other fluids, into sample with cracks or shear bands. In such cases, the choice of the regularization basis would require particular attention (for example a finite element mesh may be suited, and a model based approach may also open a promising route).

Acknowledgment

Clément Jailin would like to especially thank the members of the Division of Solid Mechanics of Lund University for their warm welcome and for giving him the opportunity to work on this subject.

Chapter 6

Vibration measurements in a tomograph

This short chapter in an unpublished ongoing project that aims to develop a modal basis identification method for a sample with in-situ measurements.

An in-situ vibration mode measurement method has been developed based on Projection-based Digital Volume Correlation techniques. Three projection angles are selected and a large number of radiographs of the vibrating sample are acquired at random and unknown times with a small exposure time in order to 'freeze' out the displacement, without motion blurs. The measurement, based on an initial reconstruction acquired in a static configuration, is performed using Proper Generalized Decomposition techniques. All projections are related with the known microstructure and the unknown modal basis associated with amplitudes. The key lies in being able to give a (weak) relationship to link those amplitudes together. A statistical a priori was used assuming that the modal participation is equal for each angle. The method was applied on two synthetic test cases where the volume deformation (using chosen mode shapes and random amplitudes) and projections where numerically performed. In both cases, the space modes were recovered as well as the associated time amplitudes allowing the method to be validated.

Introduction

Full field measurement is a key element for model validation and identification. The measurement of the displacement field during an in-situ test in a tomograph is indeed very rich [1]. With the aim to identify or validate a model, the experiment can be designed to be discriminant, i.e., highly sensitive to parameter changes. Various techniques can be used to enhance sensitivity: adding texture in the volumes (e.g., adding inert inclusions to the materials, paint or surface marking, introducing surface roughness, etc.), optimizing the shape of the sample (e.g., with holes or notches to enhance a plastic behavior as in [2] or more sophisticated method based on topology optimization [3]), the loadings and crack paths (e.g., in [4,5]), the choice of more and more developed experiment procedures and test (e.g., multi-axial, coupled loading with thermal, electrical or magnetic actions, cyclic load/unload, vibrations) etc.

When designing an experiment, temporal resolution may be an issue. For example of experimental time constants, a fatigue test is usually performed at say 5 to 30 Hz depending on the material dissipation , a thermal test may have it variations occurring in few seconds, creep in minutes. On the opposite, standard 3D full field measurement techniques (e.g., Digital Volume Correlation [6] (DVC)) require the acquisition of several volumes, each taken, for example in a lab-CT, minutes or hours. Different methods have been developed to be able to acquire tomography volumes in some of those fast applications. In fatigue tests for example, the sample is loaded cyclically and the test is interrupted at some steps for the scanning process.

An extension of DVC where the measured quantities are identified on projection is the Projection-based DVC (P-DVC) [7,8,9,10]. Those developments allows to follow a 4D (space/time) kinematics from a single projection per state. The developed P-DVC techniques could also be applied to follow the cyclic deformations, radiograph after radiograph, while the sample is rotating. Nevertheless, the measurement of vibrations (say from 10 to 100 Hz) is not feasible because the experiment can not be stopped. Some dynamic tomography techniques, detailed in the next part of this manuscript aim at imaging a moving sample by capturing the kinematics and reconstructing the sample at the same time. However the vibrations are too fast and could not be captured by such approaches (especially if the sample has a high damping and hence can not be imaged during a long time) unless constant phase is used to trigger acquisition, and a steady state vibration regime is imposed.

Being able to measure the displacement fields (thus the modal basis) of a sample, in 3D, is an interesting challenge. More than identifying the kinematics, the modal basis of a sample is extremely sensitive to the parameters of models. The identification of 2-3 modes is hence very discriminating and would permit model (un)validation.

In [11], the modal basis of a vibrating plate was measured using a deflectometry setup with a random and unknown (but numerous) sampling in time. The collection of images were acquired using a short exposure time to freeze the displacement without motion blurs. Without time information, frequencies cannot be measured. However, images are related to each others by a kinematic field that is often composed of the first modes. From DIC analysis of those images, a weighted Principal Component Analysis could extract those modes that were compared with numerical results.

A measurement method is proposed for an in-situ vibrating sample based on a P-DVC procedure with reference volume in the same spirit as [11]. Three projection angles for the multi-view P-DVC are selected and a standard static acquisition is performed in a reference configuration. A large series of radiographs ( 150) is acquired at each angle, at different and random vibrating states, using a short exposure time. To mention a realistic order of magnitude, 1/72 second is available in the LMT lab-CT. The kinematics identification is then performed with a PGD procedure mode by mode using statistical assumption(s) on the displacement field. Two numerical applications show the accuracy of the proposed method with a sample that is virtually deformed with 1 or 2 modes. With the proposed approach, the identification enables the modal basis to be retrieved.

Modal measurement method

The method is based on a 'multi-view' P-DVC procedure as described in the above Chapter 2 and in [12]. With f (x) the reference image, Π θ the projection operators [13] at angle θ, p(r, θ) the recorded deformed projections, acquired at random times using the shortest exposure time possible, the functional to minimize is

Γ Vib (u) = 1 β f r,t Π θ(t) [f (x + u(x, t))] -p(r, t) 2 (6.1)
with the implicit convention that the projection Π θ(t) of point x on the detector plane is r.

The displacement field can be written on the modal basis as

u(x, t) = i α i (t)Φ i (x) (6.2)
with Φ i (x) spatial vibration modes, indexed by i, and α(t) the temporal amplitude of each mode. Because the time acquisition is not known in our experiment, this amplitude can not be regularized as sinusoidal functions (as we would do for experiment with a fine sampling). The modal basis can be written on a finite element mesh kinematics, using mesh shape functions ψ j (x) as usually done in global DIC or DVC procedures [14]:

u(x, t) = i α i (t) j u i j ψ j (x) (6.3)
In a PGD framework, each mode is identified successively (see Chapters 3 and 5). Thus

u k (x, t) = u k-1 (x, t) + α k (t) j u k j ψ j (x) (6.4)
The time evolution can not easily be written on a priori time basis as it is acquired at random time. However, every additional sensor information (i.e., laser interferometry, acoustic measurements) could be used as regularization.

Vibration displacements are supposed to be small and centered around the reference configuration (i.e., null mean displacement), hence the volume is not advected for each deformed position. Thus the computed sensitivities do not have to be updated and are computed once. This allows huge gains in computation time as the linearization is only performed once at the reference state, no update are performed. The linearized functional is hence, with β p the noise variance for the entire ROI:

Γ k Vib (u k , α) = 1 β p r,t ρ(r, t, u k-1 ) + Π θ(t) ∇f (x) i u k i ψ i (x) α k (t) 2 (6.5)
with the projected residual fields ρ(r, t, u) = Π θ [f (x + u(x, t))] -p(r, t). It can be noted that because the volume is not advected, the residual fields can be computed simply as the difference between deformed and reference projections. This mere difference allows one to avoid the possible reconstruction and projection artifacts (e.g., inaccurate reconstruction parameters). Because the above space -time functional is coupled, it is proposed to solve it using a fixed point algorithm. The two minimizations are

u k = Argmin v k Γ k Vib (v k , α) (6.6) α = Argmin a Γ k Vib (u k , a) (6.7)
and gives the following linear systems, with l indexing the different N l angles, each composed of N t projections

u k i = N l l=1 H l ij -1 N l l=1 h l j (6.8)
with H l the hessians of Γ Vib minimized with respect to u k for angle θ l , indexed by l

H l ij = r,t l α k (t)Π θ l [f (x)ψ i (x)] α k (t)Π θ l [f (x)ψ j (x)] (6.9)
and the second member

h l i = r,t l ρ(r, t, u k-1 ) α k (t)Π θ l [f (x)ψ i (x)] . (6.10)
The key to solve this problem is an assumption on α(t) because the acquisition times are randomly sampled. It is proposed to resort to a statistical property that, with a large number of snapshot, the 3D motion is similar for any angle hence

t l α k (t) 2 = 1 with
the normality by convention. It is noteworthy that this condition is applied on the 3D displacement and it is not similar to its expression in residual fields. If for example the motion is parallel to a projection axis, i.e., zero sensitivity in a perpendicular projection direction, the 3D amplitude is however the same for the different X-ray views. With this assumption the Hessian matrices can be written

H l ij = N t r (Π θ l [f (x)ψ i (x)]) (Π θ l [f (x)ψ j (x)]) (6.11)
What is interesting is to see that this Hessian is exactly the same as the one described in the Chapter 2 and in [12] with the multi-view P-DVC. Indeed the assumption on the amplitude allows one to treat this measurement as if it were a static case (at least for the Hessian). The system is hence well conditioned and the sensitivity to every node is recovered with the different angles.

The very small problem in time is

α k (t) = N l l=1 M l -1 N l l=1 m l (6.12)
with M l the scalar Hessians of Γ Vib minimized with respect to α for angle l

M l = r Π θ l f (x)Φ k (x) Π θ l f (x)Φ k (x) (6.13)
and the second member

m l (t) = r ρ(r, t, u k-1 ) Π θ l f (x)Φ k (x) (6.14)
An interesting property of this PGD analysis is that it may be intrinsically linked to the modal measurement. As explained in [15,16,17], the Principal Component Analysis (PCA) of an infinite collection of displacement fields of a vibrating system corresponds to the modal basis if the mass is an identity matrix. This can be easily generalized to an arbitrary distribution considering the PCA decomposition of M 1/2 u instead of u so that the L 2 norm of this quantity is equivalent to an energetic norm. The obtained modes are also orthogonal with respect to the mass product. In this tomography application, the minimized quantity treated by PGD is the residual field. This quantity could be multiplied by the image projection itself such that it would be exactly corresponding to a weighted by the gray level intensity, in other words the exact mass of the volume. The measured PGD modes should hence converge to vibration modes. The corresponding functional would thus be

Γ Vib (u) = 1 β f r,t Π θ(t) [f (x + u(x, t))] -p(r, θ(t)) 2 • Π θ(t) [f (x)] (6.15)
Another assumption, optional, not required, is to impose that the mean amplitude of the 3D displacement is null, for each angular sector. A Lagrangian multiplier λ(t) can be added to the functional

Γ Vib-Reg (u, α) = Γ Vib (u, α) + l λ(t) t α(t) (6.16)
This assumption is used in the second application.

The modal measurement method is summarized in the following algorithm (obeying a normal distribution to be realistic).

First test. The modal basis of the first test case is composed of a single horizontal mode, with e x and e y the vector of the x -y plane and z the normalized vertical coordinate (ranging from 0 to 1). The standard deviation amplitude for each x and y direction is 0.1 voxel.

Φ(x) = 0.1 • z • (e x + e y ) (6.17) 150 projections for each angle are generated at random amplitude of the mode with the addition of a white Gaussian noise whose standard deviation is 0.1% of the gray level dynamic of the projections (to be compared with the very small signal due to small displacements). A projection of the volume and two residual fields are shown figure 6.3 and 6.2. It can be seen in the residual fields patterns that the amplitude are different because they have been acquired at random time. The direction is also different (the positive and negative signature of the motion have an opposite sign).

The measurement is performed with a coarse T4 mesh (4 nodes element with tri-linear interpolation) composed of 16 nodes. This corresponds to 48 degrees of freedom in space and 450 in time. The measured displacement field of the first identified mode, after 6 iterations (measured field and α(t)) are compared with ground truth and shown figure 6.4 and 6.5. The displacement field corresponds to an 'in-plane' motion whose amplitude evolves with z. The amplitudes of two top nodes are not perfect as they are related to low sensitivities (tangential to the cylindrical shape). Adding local mechanical regularization or using models should help converging on more mechanically admissible fields. The identification of the amplitude is really close to the reference. The plot has to be read in three parts of 150 projections. Each third corresponds to a different projection angle. It can be seen that the error is more important for the third part as the sensitivity is low (incidence angle of 30 • ). Although the error is high in amplitude, it has no consequences in residuals, which is the minimized quantity. In fact, all angles do not need to be sensitive to the mode kinematics as long as it can be identified with others. This is the reason why three projection angles have been selected. This error is thus expected and inevitable.

Second test. The second test case is composed of two modes: a first one which is a vertical extension mode Φ 1 and a quite similar horizontal mode as in the first test case Φ 2 . The incidence angle for the third projection angle is nearly 0 and should result in an almost (not exactly 0 and the beam is not parallel) null sensitivity.

Φ 1 (x) = 0.5 • z • e z (6.18) Φ 2 (x) = 0.1 • z • (e x + 1.5 • e y ) (6.19)
150 projections are also used and corrupted with the same white Gaussian noise. The assumption of zeros mean amplitude is used and implemented in the standard procedure as a normalization step. Two residual fields acquired for the same projection angle and different times are shown figure 6.6. The first one is mainly composed of a large vertical motion, the second one displaying an essentially horizontal motion. The second mode displacement field and time amplitudes are shown figure 6.9 and 6.10.

With this second correction, the norm of the residuals decreases from 447 to 402. It can be noted that the third angle is, as in the first application, not associated to high sensitivities. The captured amplitudes are therefore not in accordance with the input signal but have a small influence in the gray level minimization. Figure 6.9: Measured displacement field Φ 2 (x) for the second mode of the second test case.

Conclusions and perspectives

An in-situ vibration mode measurement of a sample based on P-DVC has been developed. After the acquisition of a reference volume, a large series of projections are acquired at few angles using a fast exposure time that freezes out the displacements without motion blurs. Because the vibration frequency may be much higher than the acquisition rate, the radiographs can not be related in time and the concept of frequency is out of reach. The time sampling is hence considered random and unknown (although other measurements could help providing this missing information). What links all projections is the microstructure (known) and the modal basis (unknown). Based on a PGD approach, each mode is identified successively using the statistical assumption that the norm of the 3D displacement amplitude is constant for each angle of acquisition (valid for a statistically steady excitation, and a sufficiently abundant sampling).

The method has been applied to two numerical test cases. Each test is generated by deforming a volume with vibration modes associated with random amplitudes. The numerical experiment is performed with 3 angles and 150 projections per angle. A first test is composed of a single horizontal mode. The kinematics of the second application is composed of two different vibration modes. In both cases, the measured displacement modes (i.e., field and amplitude) are compared with the ground truth. Both applications allow us to validate the modal identification procedure. This proposed procedure is very light to compute as the volume is never updated and the sensitivity fields computed only once. A fine mesh could easily be designed and is not a limit (and could be associated with some Tikhonov [18] or local elastic regularizations [19] for example). In real cases, where the displacement would be very large at some stages, a multi-scale procedure filtering out high frequencies in the microstructure could be envisioned, or the corresponding projection library could be pruned out by erasing those cases from projections.

From the Rayleigh-Ritz ratio, a coupling of the obtained results with assumptions on the stiffness [K] and mass [M ] (given by the tomography though it has to be correctly weighted) matrices enables the frequencies of the measured modes to be identified, such that, with H the spatial Hessian matrix as defined previously that refers to the measurement uncertainty (positive definite):

ω k 2 = Φ k H 1/2 [K]H 1/2 Φ k Φ k H 1/2 [M ]H 1/2 Φ k (6.20)
In the two examples, the lack of sensitivity for two projection angles resulted in an identification of the displacement amplitude that was more or less accurate depending on the angle. However although the measurement is not perfect, because of a low sensitivity, it does not affect the residual field thus the minimized quantity. Designing a test case, i.e., choice of angles, surface texture, number of projection per angle, is an interesting perspective. This prior optimization could be performed with respect to an initial numerical model (e.g., an elastic model).

Implementing a multi-scale approach may be crucial when dealing with large displacements, especially because the volume is never advected. The procedure would be the same as the one described in other Chapters: (i) low pass filter to smooth the functional shape (convolution with a Gaussian kernel), (ii) downsizing to reduce computation time. As performed in previous chapters, this pyramidal identification would be performed from the coarsest to the finest scales.

Performing an in-situ measurement in a lab-CT is also an exciting perspective. As a little teaser with order of magnitudes: the exposure time have to be approximately say 5-10 times shorter than the period to ignore motion blur. For the LMT lab-CT, an experiment to measure modes at 5-15 Hz could be designed. The sample could be soft, in gels or silicons to have low frequency modes and the excitation (white if possible) performed with for example acoustic excitations, air puffs, shocks, etc.Moreover, with the use of a low absorbing sample, a X-ray beam chopper and an intense beam (i.e., a W source) one could highly reduce the acquisition time and catch much faster phenomenon.

Because P-DVC has very strong links with stereo-steps, this modal measurement approach could also be implemented with 'camera projection' and surface measurements. A stereo-DIC setup composed of a single camera moved at different places for the acquisition of snapshots could be designed.

Being able to couple this measurement with other modalities or sensor is promising. With acoustic measurements or LASER velocimeter (based on Doppler effect), the frequencies could be captured and would help measuring the displacement field with or without models.

This developed method is insensitive to damping. Indeed if a large damping occurs, the acquisition can be performed very quickly after the excitation. Moreover, repeating the mode determination at different delay time after the excitation could provide modal viscous damping.

Finally, the identification of material parameters could be performed with this modal P-DVC framework. With a very high sensitivity, sample shapes could be designed and tested. An identification procedure, for example based on a simple elastic modal computation, would allow to identify both kinematics and parameters. An online calibration procedure for X-ray lab-CT is developed using projection-based digital volume correlation. An initial reconstruction of the sample is positioned in the 3D space for every angle so that its projection matches the initial one. This procedure allows a space-time displacement field to be estimated for the scanned sample, which is regularized with i) rigid body motions in space and ii) modal time shape functions computed using model reduction techniques (i.e., proper generalized decomposition). The result is an accurate identification of the position of the sample adapted for each angle, which may deviate from the desired perfect rotation required for standard reconstructions. An application of this procedure to a 4D in situ mechanical test is shown. The proposed correction leads to a much improved tomographic reconstruction quality.

Part II P-DVC without reference reconstruction

Introduction

Micro-Computed Tomography (µ-CT) is widely used for non-destructive imaging [1,2]. It consists in reconstructing full 3D volumes of a sample from sets of radiographs and a reconstruction algorithm [3]. From the first developments with parallel X-ray beams in synchrotron facilities, it is now an accessible equipment in laboratories. Recent developments allow for imaging at high spatial and temporal resolutions [4,5].

To deal with the reconstruction procedure (i.e., an inverse Radon transform), different algorithms exist such as filtered back projection (FBP), its (Feldkamp-Davis-Kress) extension [6] for cone beam geometries or algebraic methods [3]. All these reconstruction methods require the geometric parameters of the tomograph (e.g., position and orientation of the source, rotation axis and detector, pixel sizes) to be known. If the geometric parameters are not well calibrated or if they change in time, the reconstruction suffers from artifacts that may create blur or magnification effects. They may forbid the quantitative use of the acquired volumes (e.g., space-time full-field measurements [7]).

Various types of methods exist to estimate the 9 geometric parameters for cone-beam X-ray CT [8]. These methods can be separated into two categories, namely, offline and online calibration.

Offline procedures are calibrations using a perfectly known pattern with radiopaque markers. Virtual projections of this geometry have to match the acquired projections. In Refs. [9,10,11], the known patterns are composed of few steel balls whose orbits are measured during the rotation of the turntable. In Refs. [12,13,14,15], the authors used an iterative procedure to identify geometric parameters from a known complex phantom and a single projection. These offline methods, which are based on a known pattern or phantom, often have to be performed before or after the real test. The parameters are generally identified only once, even for a multi-acquisition test (e.g., a 4D in situ mechanical test) in order not to move the sample during repeated scans. The parameters are therefore assumed to be stable in time. However, in some cases, this assumption may be violated and lead to blurred reconstructions. A variety of causes can be listed, such as motions of the sample or the entire set-up due to unsteady thermal expansion [16], creep (i.e., delayed mechanical response [17]) of the sample under load, motion of a mechanical testing machine due to its compliance, accidental hitting, uncompensated backlash).

Online calibration is the measurement of the geometric parameters directly using the projections by exploiting knowledge of the object being imaged. In the spirit of offline procedures, Andò [18] used a known frame attached to the sample to find the geometric parameters of the scan during acquisition. Pannetta et al. [19] proposed a calibration method based on the minimization of the quadratic difference between projections at a certain angle and the mirror projection at the opposed angle (for parallel beam). Kyriakou et al. [20] developed an optimization of geometric parameters based on the minimization of the information entropy of the reconstructed volume. With 3D-2D image registration [21,22], Ouadah et al. [23] used a pre-scanned 3D image as the reference volume (hence a known very complex pattern). In Ref. [24], a motion tracking system allows for measuring and then compensating the in situ motion of the acquired image. A recent publication of Sun et al. [25] proposed an optimization procedure based on the initial reconstruction of the volume with a first set of geometric parameters. Then the minimization of the "projection residual" fields (i.e., the difference between the projected reconstructed volume and the initial sinogram) with respect to rigid body motions of the volume, provides an estimate of the motion of the scanned sample. This motion is finally used for the correction of projections. For the application aimed by the authors, i.e., patient motion compensation for medical imaging, the measured displacement field is large (5-20 voxels as an order of magnitude), and cannot be smaller due to several assumptions or approximations (on the sensitivities, on the axial displacement supposed to be null, on filtering).

A robust and accurate method that measures sub-pixel or sub-voxel displacements of patterns between different images is provided via Digital Image Correlation (DIC) in 2D [26] and Digital Volume Correlation in 3D [27,28]. A recently developed projectionbased Digital Volume Correlation (P-DVC) procedure, in the same spirit as 3D-2D image registration, lies in between 2D DIC and 3D DVC as it deals with 3D displacement fields of the reconstructed reference sample with respect to deformed radiographs with the minimization of a set of 2D error maps. This method was initially developed for the measurement of uniform displacement fields during different steps of a mechanical test [29], with possible cracks [30,31] and finally for identification purposes in a lab-tomograph [32]. This procedure leads to huge savings in acquisition time. However, as such, P-DVC uses a reconstructed 3D image to evaluate its motion, and is not designed to identify geometrical reconstruction parameters.

In this paper, an online two-step iterative procedure is proposed. First, the rotation axis is identified in the same spirit as Ref. [25], from a pattern that is the loaded sample itself, reconstructed from its projections hence without a perfectly known geometry. Thanks to P-DVC, the sample is repositioned in the laboratory frame so that its projection matches the corresponding acquired projections for every angle. The result is an identification of the sample position at each instant (i.e., projection angle) with a sub-voxel uncertainty and thus projections can be corrected to coincide with nominal geometrical parameters, and allow for reconstructing a new volume and iterate if needed. A major advantage of this procedure is that an imperfect kinematics for the acquisition, such as a precessing rotation axis, can be accounted for, thereby enhancing the reconstruction quality without any additional demand on the acquisition protocol.

The displacement field is decomposed over a basis of separated functions, namely, spatial modes (here chosen to be rigid body motions) and time (or projection angle) modes. Even when restricted over few spatial modes, a large temporal flexibility leads to a poorly conditioned system. Two routes for dealing with such issues are explored. Either few temporal modes are introduced (with the difficulty of choosing a suited temporal basis), or a model reduction technique is used. The spirit of the latter, which is inspired from the Proper Generalized Decomposition (PGD) technique, consists in a progressive enrichment of the space-time modes for displacement corrections. Based on the minimization of the projection residual, a greedy approach extracts one mode per iteration. This model reduction technique reveals very instrumental to determine the displacement modes.

The method based on 4D P-DVC is introduced in Section 7.2. In Section 7.3, the application to an in situ mechanical test with a cone-beam geometry is performed. The rotation axis is determined with a single iteration of the procedure and the correction of the axis shows significant improvement of the reconstructed volumes.

Method

General principle

In the following, raw radiographs acquired by the CT-scanner are first normalized by their flat fields (i.e., radiographs acquired with the same conditions as the actual scan but without the specimen) and their logarithm computed (i.e., Beer Lambert law) to produce "projections" that constitute the sinogram. It is from such preprocessed data, which will be referred to as "projections," that reconstructions are performed using the Filtered Back Projection (FBP) algorithm (or suited extensions for cone-beam geometries).

Starting with an initial reconstructed 3D image, hereafter also called "volume", it is positioned in the 3D space for every projection angle and its projection is computed so that it matches at best the corresponding projection. This difference between initial projection and projected volume is called projection residual field. Let us note that these residuals can be interpreted as the sum of four contributions due to:

1. acquisition noise and artifacts (e.g., beam hardening, phase contrast, approximations during reconstruction and projection);

2. approximate (or erroneous) geometric parameters of the tomograph in the reconstruction and projection procedure;

3. rigid body motions of the scanned sample during acquisition (i.e., deformation of the testing machine or turntable, uncompensated backlash);

4. deformations of the sample itself (e.g., thermal effects, relaxation). This latter effect is assumed to be null (or extremely small) in the examined case.

Projection-based digital volume correlation

The proposed method for the identification of the displacement field based on radiographs is Projection-based Digital Volume Correlation (P-DVC) [29]. The first step of the procedure is the acquisition and reconstruction of a reference volume, f (x), from a complete set of N θ projections s(r, θ) and an initial choice of geometric parameters. In the present notations, x is a 3D vector in the sample frame, r is a 2D vector denoting positions in the detector frame, and θ the rotation angle. The idea is to estimate the 4D displacement field u(x, θ) so that each projection of the displaced volume f (xu(x, θ k )) matches the initial projection, s(r, θ k ). After updating the reconstruction, from the corrected displacement, the procedure is iterated until corrections of the geometrical parameters are less than a predetermined threshold.

In the same spirit as global DVC [33], the displacement field is obtained from the minimization of the quadratic difference, χ 2 u , between the (re-)projected 3D image corrected by the displacement field Π k [f (xu(x, θ k ))] and the initial sinogram for every angle s(r, θ k ), with Π k being the conical projection operator along the θ k direction

χ 2 u = 1 N θ |Ξ|γ 2 s k,r Π k [f (x -u(x, θ k ))] -s(r, θ k ) 2 (7.1)
where the double sum over (k, r) stands for the discrete integration over all pixels r ∈ Ξ of the detector (or its utilized part, which will be referred to as region of interest or ROI) and all N θ projection angles θ k . γ 2 s denotes the variance of the sinogram noise, and |Ξ| the area (number of pixels) of Ξ. The normalizing parameters of Equation 7.1 are chosen in such a way that when the residuals ρ = Π k [f (xu(x, θ k ))] -s(r, θ k ) are only associated with acquisition noise, χ u = 1. Any deviation from 1 (i.e., χ u > 1 is an indication of model error (i.e., the chosen kinematic basis is not fully consistent with the studied experiment).

When small displacement levels are assumed compared to the correlation length of the imaged microstructure, a small perturbation expansion is written about the current estimate ũ of the displacement field, with ∇f (x) the gradient of the 3D image

χ 2 u ≈ 1 N θ |Ξ|γ 2 s k,r (Π k [f (x) -∇f (x -ũ(x, θ k )) • δu(x, θ k )] -s(r, θ k )) 2 (7.2)
so that the minimization operates on a quadratic function of the unknowns, thereby leading to easily accessible search directions. It is noteworthy that after each evaluation of the displacement corrections δu, a correction of the volume is performed so that Equation 7.1 is used without approximation.

4D Regularization

Different regularization procedures of the displacement field have been introduced in the literature for global DIC methods where the kinematics is expressed on a finite element mesh. Spatially, a local mechanical elastic constraint, in 2D [34,35] or 3D [36,37], strong regularizations or integrated methods with a reduced basis composed of elementary fields from mechanical computation [38,7]. These regularizations lead to a drastic reduction in the number of unknowns and enable for seamless experimental/numerical procedures. The displacement field is expressed as a combination of N m spatial modes Φ i (x) that are weighted for each angle by α i (θ k ). Such a separated expression is standard practice and as such bears no consequence, but it will reveal convenient for the following model reduction technique

u(x, θ k ) = Nm i=1 α i (θ k )Φ i (x). ( 7.3)
The spatial modes are expressed as a reduced basis composed of the 6 Rigid Body Motions (RBMs) ψ j (x)

u(x, θ k ) = Nm i=1 6 j=1 α i j (θ k )ψ j (x). (7.4)
The time changes α i j (θ) can similarly be constrained by additional temporal regularizations [39,40]. This constraint may come from previous knowledge or anticipation of the motions. The N n temporal modes are expressed as elementary functions σ n (e.g., Dirac distribution if no regularization (free) is considered, polynomials, sinusoidal functions)

α i j (θ k ) = Nn n=1 a i jn σ n (θ k ) (7.5)
where a i jn are the amplitudes of elementary basis functions being the product of the N n temporal and N m spatial modes. In the treated application, the time basis does not change with mode identification thus does not depend on i, but such cases could be designed. The displacement field is finally written as

u(x, θ k ) = Nn n=1 6 j=1 a jn σ n (θ k )ψ j (x). (7.6)
Two methods are proposed to solve the minimization problem, depending on the complexity and number of degrees of freedom:

• A full determination of the displacement field, in Section 7.2.3. This method is the most complete but requires an appropriate prior knowledge of time and space basis functions, which may involve many degrees of freedom, and hence possibly poor conditioning. It is referred to as "full identification" in the following.

• A PGD approach, where additional modes are successively determined as long as the residual level is considered too high to be explained by noise (i.e., χ u > 1). This method is presented in Section 7.2.3.

Full measurement

Because the x referential frame is linked to the rotating sample, a choice of time functions can, for example, be low order Fourier modes, plus possibly linear motions to account for a slow and steady drift that would break periodicity. The full determination of all space and time amplitudes, a i n , (the N n time modes times the N m = 6 spatial modes) may be accessible; the risk being here that the system become poorly conditioned. The following subsection will describe a method to reduce the number of modes to only those that are needed.

The problem is rewritten as

u(x, θ k ) = n,j a jn σ n (θ k )ψ j (x) (7.7)
so that its minimization is achieved with Newton's descent method. A single uncoupled system is obtained in terms of the corrections δa in

δa in = T -1 injm t jm (7.8)
with

T injm = k,r σ n (θ k )Π k [ψ i (x ) • ∇f (x )]σ m (θ k )Π k [ψ j (x ) • ∇f (x )] (7.9)
and

t jm = k,r (s(r, θ k ) -Π k [f (x )])σ m (θ k )Π k [ψ j (x ) • ∇f (x )] (7.10)

Free modal measurement

An integrated approach is proposed along the line of Proper Generalized Decomposition (PGD) techniques [41,42,43], which consists in successive enrichments of the displacement field u(x, θ k ) summing an additional contribution at each iteration, each term of the sum being sought a priori in a separate representation. PGD-DIC and PGD-DVC [44,45] with one-dimensional space functions is here extended to 3D space-time (angles) analyses.

In the following progressive PGD procedure, the spatial modes, which are defined as Φ i (x) = p i j ψ j (x), will be identified successively, one per iteration, with a greedy approach [46].

u l (x, θ k ) = u l-1 (x, θ k ) + Nn n=1 a l n σ n (θ k )   6 j=1 p l j ψ j (x)   . (7.11)
Let us note that only the product a l n p l j matters so that an additional convention (but just 1) such as a l = 1 or p l = 1, could be freely chosen without consequence.

A fixed point algorithm is used to get the solution. Alternate minimizations of the two unknown vectors a l and p l are proposed. The minimization of the functional leads to the determination of the unknowns with two coupled equations p l = Argmin p l (χ u (α, p) 2 ) (7.12)

a l = Argmin a l (χ u (α, p) 2 ) (7.13)
i.e., minimization of χ 2 u with respect to the additional mode is considered. A general overview of the 4D PGD P-DVC procedure is shown in the algorithm 5. Even though a maximum value of iterations or convergence criteria p , and α can be enforced to stop the fixed-point algorithm, this revealed unnecessary as the maximum number of iterations to reach stagnation is usually quite low (i.e., 3-5).

Algorithm 5 General 4D-P-DVC fixed-point procedure while High residual norm do Initialize a l and p l Correction f (x) ← f (xu l-1 ) while ∆Φ i < p and ∆α l < α do Compute spatial mode p l , Equation 7.14 Compute temporal amplitude a l , Equation 7.17 end while Update displacement field u l , Equation 7.11 l = l + 1 end while

The two parts of the fixed point algorithm are obtained from the above linearized functional using Newton's scheme. The derivative with respect to p leads to

p l = N -1 n (7.14)
with N the spatial Hessian matrix of χ 2 u with respect to p (i.e., N ij = ∂ p i ∂ p j χ 2 u ) and n the second member vector based on the residual field

N ij = k,r α l (θ k ) 2 Π k [ψ i (x ) • ∇f (x )]Π k [ψ j (x ) • ∇f (x )] (7.15)
and

n i = k,r α l (θ k )(s(r, θ k ) -Π k [f (x )])Π k [ψ i (x ) • ∇f (x )] (7.16)
where x is the advected position of the volume with the previously identified modes such that x = xu l-1 . Similarly, the derivative with respect to a l leads to

a l = M -1 m (7.17)
where, as previously, M is the temporal Hessian matrix (i.e., M ij = ∂ a i ∂ a j χ 2 u ) and m the second member vector based on the residual field

M ij = k,r σ i (θ k )Π k [Φ l (x ) • ∇f (x )]σ j (θ k )Π k [Φ l (x ) • ∇f (x )] (7.18)
and

m j = k,r (s(r, θ k ) -Π k [f (x )])σ j (θ k )Π k [Φ l (x )∇f (x )] (7.19)
The residual field at convergence, ρ(r, θ k ) ≡ (s(r, θ k ) -Π k [f (x)]) for the N θ 2D projections, gives a very precious information on the quality of the solution in the projected domain, which is parameterized by r. Its norm is the minimized quantity and defines the convergence criterion. Ideally, it consists of acquisition noise and artifacts. It is thus a natural place to judge whether some spurious displacement is to be corrected for.

Comparison metric

In order to show the improvement provided by the method, comparisons at 2 steps are performed. The first one is the measurement of the norm of the residual between the projections of the final reconstructed volume and the acquired radiographs. This norm is the minimized quantity of interest hence this is the true metric of the present procedure. The Signal to Noise Ratio (SNR) can also be defined to evaluate the residual quantitatively. The higher the SNR, the better the solution. It is defined as 20 times the decimal logarithm of the ratio of the standard deviation of the projections σ(s) over that of the residual fields

σ(ρ) SN R(t) = 20 log 10 σ(s) σ(ρ) (7.20)
One may also be interested in the quality of the reconstructed volume itself. A visual estimation of the sharpness of the image can be supported by the value of the Shannon information entropy that has to be minimized,

S = - f p(f ) log 10 [p(f )] (7.21)
where p(f ) is the gray level distribution (i.e., probability of observing a gray level equal to f ).

Application

Test case

The application case where the set of radiographs has been extracted is an in situ tensile test on a cast iron sample. The sample (Figure 7.1(a)), which was mounted in an in situ testing machine similar to that used by Buffière et al. [47], was scanned at LMT lab-tomograph (X-View X50-CT, North Star Imaging, 129 kV, 95 µA, W target). The voxel size at full resolution was set to 2.8 µm. The complete scan of the reference state consists of 1,000 radiographs captured at equally spaced angles ranging over a full 360°r evolution. Two flat-fields are acquired after conditioning and before the experiment in order to perform flat-field corrections. One dark-field has also been acquired before the experiment. Each radiograph is averaged over 30 frames in order to reduce acquisition noise. The acquired radiographs have a definition of 1944×1536 pixels 2 but the following procedure deals with two lower scales. Coarse graining into superpixels of size 2×2 (and 4×4) is carried out with the convolution of the images by a Gaussian kernel with a characteristic width of 2 (and 4) pixels. The resulting image is downsampled over a coarse 2×2 (resp. 4×4) regular square grid to create smaller images (called image at scale 2 and scale 4). The projections are obtained after flat field normalization and standard beam hardening correction [48] due to the high absorption of the cast iron with a third order polynomial.

Reconstructions and projections are performed with ASTRA toolbox [49]. It is obtained with Feldkamp-Davis-Kress (FDK) procedure suited for cone beams [6]. The geometric parameters for the initialization are given by the offline (i.e., standard) calibration of the tomograph that uses a pattern composed of vertical calibrated steel balls. This offline procedure is very fast (i.e., approximately one minute) and has been performed once at the beginning of the experiment. It can be noted that this initialization is not needed and could be replaced by the use of the proposed approach only.

The entire test of 16 loading steps has been performed over few days because each acquisition requires about 2 h. The acquired projections and the reconstructed volume are shown in Figure 7.1(a,b). Because of cone beam artifacts at the top and bottom of reconstructions, the projected ROI definition is 90 < z < 897 pixels for scale 2, and 45 < z < 449 pixels for scale 4.

Overnight, the grip displacements are held at a fixed value and nothing is supposed to move. In spite of this assumption, a DVC analysis of the test shows large motion and The proposed application procedure deals with the two scans in order to find the correct geometric parameters. The two previous identification methods are used. The first test case, in Section 7.3.2, is step 2 at scale 4. For this calibration the correction is large hence it may require many displacement modes. However because the total number of degrees of freedom is small (i.e., 6 RBMs × 4 time functions), a full measurement is performed. The second application, in Section 7.3.3, is the correction of step 1 at scale 2. Because the displacement level is not very high, the kinematics is captured with the modal decomposition constrained by time functions. It is shown that three modes are sufficient to capture almost all the kinematics. Instead of working with the 1,000 radiographs, for computation cost and because no high frequency displacement is expected, a subset of 100 regularly spaced radiographs is used.

As previously discussed, the amplitude vector a can be scaled by an arbitrary factor if p is scaled by its inverse. To make the results comparable and easily understandable, this scale factor was chosen so that α(θ k ) varies within the interval [-1, 1]. The associated space mode amplitude thus corresponds to the maximum displacement amplitude.

Full identification

The first test case is step 2. Overnight, the geometric parameters changed. Because large motions are anticipated, the analysis is performed at scale 4 and a full identification (i.e., determination of all space and time degrees of freedom, as earlier defined) is performed. The space regularization used in the modal measurement approach is composed of the six RBMs, and the temporal shape functions consist of sine and cosine functions with a period of 2π, a constant offset and a linear change. The sine and cosine functions are suited to model for example a steady rotation axis in the lab frame. The linear component captures the kinematics of non-periodic motions such as sample motions. Thus, the kinematic basis is composed of 24 degrees of freedom (i.e., 6 RBMs × 4 time functions).

The initial and final residuals are shown in Figure 7.4. The mean SNR increases from 20.3 to 28.5. This huge enhancement means that a very large part of the actual kinematics is captured. The residuals on the edges of the sample are not totally corrected but may come from the blurred volume itself. The axis of rotation is extracted from the displacement field. It is plotted in Figure 7.6 with the previously identified axis at step 1. It can be seen that the displacement is composed of a large translation (i.e., radius of 3.9 voxels) and a rotation that is similar to the previous identification. The displacement consists of an important motion along the X direction because the procedure is not very sensitive to motions in the direction of the X-ray beam (i.e., the X direction of the laboratory frame is that of the central axis of the cone beam).

The corrected volume shows huge improvements in terms of visual quality (Figure 7.7) with sharp nodules and edges. The Shannon entropy (initially at 6.12) decreases to 5.79.

Modal measurement

Initial residual fields and regularization

The second test case is step 1. Because small motions are expected, the analysis is performed with a modal measurement at scale 2. The initial residual fields are shown in Figure 7.8. In these fields, the top and bottom parts are composed of high values due to cone beam reconstruction and projection artifacts. These areas are not taken into account in the identification procedure and are masked. The residual field is made of positive and negative patterns, which is a signature of displacements, and more precisely, a rotation on the top and bottom parts. The mean SNR level is 24.5.

Displacement mode identification

The identification is composed of three space-time modes. The mean SNR after the correction of each mode is shown in Table 7.1. Three modes are sufficient to capture the kinematics because the SNR does not increase much after mode 3. The angular SNR is shown in Figure 7.9. Low SNRs at every π/2 increment (i.e., alignment of the sample edges with the beam direction) are interpreted as phase contrast.

It is worth noting that the first two modes are complementary and affect different angular sectors. In order to validate the choice of the time basis, the result of the mode identification with Dirac time functions (i.e., with no temporal regularization) is shown in the same plot. This temporal measurement is performed with the same spatial modes (it can be noted that in the case of a purely free identification, the results would have been different because of different spatial modes). It is observed that, for the first two modes, the free and constrained temporal histories are close to each other. Hence the chosen time functions provide an appropriate basis. For mode 3, the vertical component of the displacement field is close to zero, so that the sensitivity is low for some projection angles. The difference of time history between the free and constrained functions is on these non-sensitive projection angles (i.e., π/2 and 3π/2). For those angles, the time regularization is important to correctly capture the kinematics. Because these angles lack sensitivity, they bear very little weights in the procedure (see Figure 7.13) and the regularized time function is not a mere least squares fit of the free form but is weighted by sensitivities. These results also show the benefit of reducing the number of temporal modes so that the conditioning of the system is enhanced and physically realistic modes are extracted.

Correction and reconstruction

The 3D displacement of the volume gives access to the position of the real rotation axis. Because it is defined in a basis linked to the rotating sample, the motion has to be projected to the laboratory frame. The mean projected displacements of the top, center and bottom parts of the sample are shown in Figure 7.14. The red line is the mean position of the sample axis during rotation. The volume can now be corrected with the projection motions. Hence a new set of projections is obtained and the reconstructed volume is updated accordingly. Figure 7.15 shows the corrected reconstruction. The information entropy for this corrected volume is 5.74 (to be compared to the initial value of 5.93), meaning that the edges and microstructural details are sharper. After this step, a new correction iteration can be performed. However, because the displacement magnitude is small, most of the kinematics was captured in this first step.

Conclusion

An online calibration procedure for cone beam X-ray tomography based on projectionbased Digital Volume Correlation has been proposed. This technique, which does not require any change from the standard acquisition protocol, repositions an initially reconstructed volume in a given 3D frame so that each of its projection matches, for each angle, the initial sinogram. The position is given by a displacement field decomposed over space (i.e., rigid body motions) and time (i.e., angular) basis functions. The identified position of the volume over time reveals the instantaneous rotation axis while scanning the sample and the drift of the set-up.

In order to preserve a good conditioning of the system, a model reduction technique based on Proper Generalized Decomposition is introduced so that key modes are progressively added to reduce the projection residuals.

The two application tests, which are extracted from a 4D mechanical test, with two different calibration procedures show very significant improvements in the reconstruction quality and much lower residual fields, thereby proving that the corrections account for most of the initial inconsistencies. With much less reconstruction artifacts, the updated volumes can be further used for quantitative kinematic measurements. In addition to the norm of the residual fields, Shannon entropy is used to assess the calibration improvement on the updated volume.

The method works with small displacements compared to the microstructure length and yields sub-pixel uncertainties. In cases of larger displacements, the 3D volume can first be updated with a coarse initialization (e.g., obtained from previous computations, multiscale procedures or cross-correlations). An update of the procedure with the corrected reconstruction could reveal necessary if the displacement is not well captured at the first iteration.

An additional improvement of the present procedure would be to use algebraic reconstruction methods for which the metric is already an L2-norm between the re-projected reconstructed volume and the available projections.

In the presented test cases, a simple spatial regularization composed of rigid body motions is used (and is sufficient) to reposition the volume. However, the method is not limited to such cases. For instance when analyzing mechanical tests with continuous loading [5], the spatial regularization could be complemented with other elementary modes including deformations (e.g., tension or compression). This could also be performed to account for viscoelasticity or viscoplasticity during the scan.

Coupled with fast acquisition devices and applied to the measurement of displacement fields during continuously loaded mechanical tests, this method could give access to ultrafast mechanical identification that could not be performed with classical means such as 3D or 4D DVC.

Acknowledgement

This work has benefited from the support of the French "Agence Nationale de la Recherche", through the "Investissements d'avenir" Program MATMECA (ANR-10-EQPX-37 grant), and the COMINSIDE project (ANR-14-CE07-0034-02 grant). The authors would like to acknowledge Dr. Matthieu Vitse for precious advices on PGD methods.

Bibliography 8.1 Introduction

Tomography is a non-destructive imaging technique that enables the visualization of the bulk of the observed specimen. Tomography is now widely used in many fields (e.g., medical imaging for diagnostic [1], biology, material science [2], etc.), performed with various waves (e.g., X-ray, neutron, electron [3], terahertz, optics, ultrasound, etc.) depending on the experiment and material absorption and or scattering. Different instruments have been developed with different flux, space and time resolutions (e.g., for X-rays medical scanners, synchrotron, lab-CT, etc.) giving access to a wide range of imaging devices and performances.

In order to image the 3D structure, the specimen rotates over 180°or 360°with respect the source-detector pair and at a series of distributed angles radiographs are acquired. Radiographs are transformed with dark-fields and white-fields, to extract the relative beam absorption, transformed with a logarithm (Beer-Lambert law) or more sophisticated treatments for beam hardening, in order to obtain so-called projections of the local coefficient of absorption of the sample. The collection of projections at all angles constitutes a so-called sinogram. Then, from the sinogram, reconstruction algorithms [4] have been developed to reconstruct the 3D imaged volume. This technique relies on the strict satisfaction of conditions, in particular concerning the geometry of the set-up and the motion of the sample as a rigid rotation with the prescribed axis and angles.

The required time for a full 3D scan varies depending on the flux (and exposition time), type of camera and rotation speed of the device. Since the beginning of the development of these techniques, the time required to acquire a tomographic scan has constantly decreased [5]. Recent papers have reported on ultra-fast tomographies, at up to 20 Hz in synchrotron beamlines, that allow extremely fast processes to be captured [6,7].

Motion of the sample during the scan is one of the main issue of tomography that leads to poor quality, blurry volumes [8]. This is the case for medical imaging (as the patient or imaged organ may move), in vivo measurements [9], for electron tomography [3] (because of the extremely small scales of observation, one cannot guarantee a fixed rotation axis at nanometer accuracy) for usually minute to hour long acquisitions, fast mechanical behavior or continuous in-situ experiments [6]. Wrong or imprecise estimates of the calibration parameters (that may even vary along the scan) can also be seen as motions in the sinogram space and have the same deleterious consequences for the volume reconstruction. Sophisticated methods have been developed to avoid or limit motion perturbations [10] especially for periodic motion, for instance using a trigger for acquisition of radiographs based on a specific signal to captures always the same phase as can be done for cardiac or respiratory motion in medical imaging [11].

Many works have been devoted to correcting imperfect acquisitions as a post-processing treatment. For automatic (re)calibration, online methods have been applied as a postprocess after reconstruction to evaluate a corrected set of calibration parameters [12]. The identified motion of the specimen is often regularized as being rigid body motions [13,14]. Very early, corrections were also applied in the sinogram space ( [15]), with affine transforms [16].

Projection-based measurement methods (e.g., Projection-based Digital Volume Correlation (P-DVC) [17], 3D-2D registration [18,19]) have been developed to correct for rigid body motions (due to a rigid patient motion or variation of calibration parameters) from the radiograph data directly.

Yet, a deforming body with a significant strain and variation in time is a much more demanding case. Projection-based Digital Volume Correlation (P-DVC) has been shown to address part of the problem with complex 4D -3D space + time -kinematic identification [20,21,22]. First if the reference 3D geometry is well known, the displacement field can be evaluated on the fly as the sample is being deformed. This method requires a high quality reference volume and a series of deformed projections. A single projection per motion state is required to capture the full 4D (space-time) kinematics. Alternatively, imperfect acquisition conditions (but no sample strain) can also be corrected using a similar technique, without a pre-determined 3D reference geometry [17], considering that the deforming projection stack is the one used for the reconstruction.

Similar developments have been carried out very early in the context of medical imaging where periodic motion is frequent (heart beat, breathing). In particular, Refs. [23,24,25,26,27] have proposed to determine the motion of the sample from projection data. Small amplitude displacement fields with a periodic modulation in time were considered and identified using highly regularized kinematic models.

However, very often, a reference reconstructed volume is known, and is used as a prior for determining the motion [23,24,28,25,27]. This is often the case for radiotherapy treatment where the key issue is to irradiate the targeted region, in spite of a spurious motion, and hence the goal is to identify the displacement field in 3D, and a fast determination is more valuable than a very precise one.

In a similar spirit [25,29,30] do not consider a reference to be known but rather use a phase signal (say from an electrocardiogram) to extract from a long sinogram projections coming from a similar phase of the motion, and reconstruct a low quality volume for a series of phase. Registration of the reconstructed volumes [29,30] allows the displacement field to be estimated and interpolated for the entire range of accessed phase. Then, backcorrecting for this motion a deformed reconstruction grid is obtained [31] on which the projection data can be backprojected using a classical FBP/FDK algorithms [4]. In this way, each ray follows the deformed sample at each projection angle. The obtained volume has a better quality than the initial one (more details and sharper edges). Ruhlandf et al. [32] recently proposed an approach along the same vein without prior knowledge of a phase for each angle, nor of a reference volume, developed for phase contrast imaging at a synchrotron facility.

In most of these studies, the displacement and strain fields between scans was relatively small (strain of approximately 1% and uniformly distributed), and often the time (or phase) is believed to be known.

The present study proposes a strategy to reconstruct both the reference geometry and its large motion from a single sinogram. No periodic signal is used to constrain the kinematics. The recorded projections are the data that drive the measurement of the kinematic field, as is proposed in P-DVC. This however requires a "model," here a reconstructed 3D volume, to be known in order to measure the displacement field. It is proposed here to "learn" this model from the projection data itself using a multiscale approach.

The standard reconstruction methods are briefly presented in section 8.2, so that the introduction of motion can be cast in a similar framework. Section 8.3 details the joint determination of the reconstructed image and the motion experienced during the scan. The latter algorithm makes use of ideas comparable to those of P-DVC for the motion, and Algebraic Reconstruction Techniques (ART) for the microstructure and exploits a multiscale approach to disentangle microstructure and motion from the sinogram. Two virtual test cases of moving samples validate the procedure (section 8.4). The first example is performed with the Shepp-Logan phantom with large deformation up to 20%. The second example is a checkerboard with a more complex temporal pulsating motion.

Motionless X-ray tomography

Tomography reconstruction is based on the relative beam intensity attenuation for each discrete detector position r = [r, z] (where z is parallel to the specimen rotation axis, and r is perpendicular to it) and rotation angle. For simplicity, and because the present paper is a proof of concept, the displacement field is assumed to lie in a plane perpendicular to the rotation axis, so that each slice z remains independent from its neighbors, and the problem turns two dimensional. Hence, only one line of the detector is considered, for a unique value of z (omitted from now on).

Let us briefly recall the principle of tomography for a parallel beam: a projection p(r, θ) is defined as the line integral of f (x) along a direction e θ , or

p(r, θ) = D(r,θ) f (x) dx (8.1)
where D(r, θ) is the line parallel to e θ hitting the detector plane at position r. Different projection and interpolation algorithms exist. In the following procedure, the Matlab function radon.m is used. Tomography consists of recording a set of N θ projections p(r, θ) for a collection of angles θ(t) as the sample is rotated over a complete (or half) rotation about a fixed axis parallel to the detector plane. For a still sample, and a continuous rotation, p(r, θ(t)), written p(r, θ), is the Radon transform of f (x), p(r, θ) = R[f (x)] and hence the f (x) can be computed from an inverse Radon transform, f (x) = R -1 [p(r, θ)]. The transpose of the projection operator, divided by the length of the intersection of the ray with the domain where the reconstruction is sought, called backprojection for any angle, θ, is written B θ .

Tomography is now a very mature field and numerous powerful algorithms have been devised in order to deal with a discrete set of angles, with fan-bean or cone-beam projections [33], with laminography, etc. However, f (x) is always assumed to stand for a rigid and still object (independent of time or rotation angle). From the collection of acquired projections, different algorithms exist to reconstruct the 3D volume [4] and fall into two categories: Fourier-domain algorithms and algebraic algorithms.

Fourier space reconstructions With Filtered Back-Projection (FBP), each projection, p(r, θ) is first "filtered" with a ramp, or Ram-Lak filter, eventually windowed. Ignoring such windowing, in Fourier space, F[p(r, θ)](k, θ) is multiplied by |k|, inverse Fourier transformed, and then back-projected in real space, thereby producing a field g θ (x) that is invariant along the direction e θ . These fields g θ (x) are simply summed over all visited angles θ, producing the sought initial image, f (x)

f (x) = N θ θ=1 g θ (x) (8.2)
Iterative reconstructions Other reconstruction methods have received much attention, namely iterative algebraic approaches which tolerate deviations from the ideal conditions of the previous Fourier space reconstruction such as for instance having access to a continuous range of angles, covering the entire half (or full) rotation. Those methods exploit the linear structure of the problem to solve, but for computational efficiency, they avoid the writing of the linear system. They are based on the minimization over volumes, ψ(x), of the functional, Γ ART [ψ], equal to the quadratic norm of the difference between the projected reconstructed volume and the acquired projections

Γ ART [ψ] = r,θ R θ [ψ(x)] -p(r, θ) 2 (8.3) then f = Argmin ψ Γ ART [ψ] (8.4)
Additional prior information may easily be added to this functional through regularization, in order to compensate limited angle range for projections, or coarse sampling for example. This generally leads to better quality reconstructions than FBP algorithms at the expense of a higher computational cost.

To solve this huge linear inverse problem, ART algorithms essentially consist of iterative updates of the volume. Successively visiting each angle, the projection of the volume is compared with the acquired one. The difference is back projected and used to correct the volume (sometimes multiplied by a damping coefficient, not considered in our case). Faster convergence rate is observed when angles are not sampled in consecutive order but rather with a large difference between successive angles. This can be achieved for instance with a permutation of the angle order. A convergence criterion on the functional value can be used to stop the number of iterations (Γ ART [f ] < ), with , a threshold value with respect to noise and artifact acquisition. Generally few iterations (N ART ) are required for convergence. The algorithm for this method is detailed in Algorithm 6.

N θ do m ← π(k) ρ(r, θ m ) ← p(r, θ m ) -R θm [f (n) (x)] ∆f (n+1) (x) ← B θm [ρ(r, θ m )] f (n+1) (x) ← f (n) (x) + ∆f (n+1) (x)
Implement additional constraints on f (n+1) (e.g., positivity) n ← n + 1 end for end while During the reconstruction procedure, additional information, defined as constraints, can be added. Those regularizations allow the reconstruction of high quality volume with few or missing angles, noisy projections and artifacts etc.. This may come from prior knowledge on the different phases of the sample (as DART algorithms proposed by [34], reconstruction with binary images [35], Total Variation [36]), dictionary learning [37], etc.. However, because those regularizations are independent from the following proposed reconstruction with motion compensation, it is not considered hereafter apart from the positivity constraint ψ(x) ← max(ψ(x), 0).

Data driven reconstruction of non-rigid samples

It is proposed to study a specimen that moves during the acquisition with a space/time displacement field u(x, t) such that, at any time, the sample is expressed with respect to a reference state f (x + u(x, t)).

For a still object f and p are bijectively related to each other, hence the introduction of motion makes the problem ill-posed. The reconstruction of the volume from the previously introduced algorithms leads to a low quality, blurry, volume.

It is to be noted that the FBP reconstruction procedure has been extended to take motion into account in [29,30,32] The driving idea is to apply the back-projection step on the currently deformed geometry of the to-be-reconstructed sample, or equivalently to transport the back-projection onto the initial geometry, unwarping the motion, so that the X-ray beam would then follow non-straight paths. In Ref. [32], the motion is estimated from the registration of two reconstructions of the volume at different instant of time and linear interpolation.

f (x) = N θ t=1 g θ(t) (x -u(x, t)) (8.5)
Because this approach requires different volumes to estimate the displacement field, it is not suited when the motion is very fast and when a single scan can be acquired.

Moreover, it is difficult to estimate a quality criterion but visual on the reconstructed volumes thus on the measured kinematics.

A recently developed Digital Volume Correlation (DVC) procedure called Projectionbased DVC [20] allows one to identify the 4D (space-time) displacement field of sample from an initially reconstructed volume and its moving projections. An extension of this method has been applied to an online calibration of the tomograph [17]. An initial (blurry) reconstruction was performed from a set of initial parameters. The comparison between the projection of the blurry sample and the acquired projections is, in addition to the acquisition noise and artifacts, the signature the erroneous projection geometry parameters that can be identified and corrected. The sample could be re-positioned for each angle by a rigid body motion. Because the motion was simple and of low amplitude, the correction could be applied on the sinogram itself leading to very significant improvement on the quality of the reconstruction. However, more complex displacements, or larger amplitudes (involving larger displacement variations perpendicular to the ray) would render the corrections on the projection inaccessible.

It is proposed to introduce here a new two-step algorithm based on ART reconstruction on the one hand and P-DVC on the other hand to identify both a complex and large displacement field and volume texture with a single scan performed on a moving and deforming object. The ART functional is naturally extended to account for the motion as

Γ motion-ART [ψ, v] = r,t R θ(t) [ψ(x + v(x, t))] -p(r, θ(t)) 2 (8.6)
where the summation over time extends over the N t acquired projections (and not just a full rotation) then (f, u) = Argmin

ψ,v Γ motion-ART [ψ(x), v(x, t)] (8.7)
The updating procedure (indexed by l) is split into two parts that are repeated alternatively:

• a volume reconstruction from an iterative dynamic ART algorithm assuming a known motion (described in section 8.3.1)

• an update of the motion from P-DVC with a given reconstructed sample (described in section 8.3.2)

However, as such, this procedure does not tolerate large displacement amplitudes. To increase the robustness and fast convergence, a multi-scale approach is coupled to the previous two-step procedure, resolving first the large scale features of both microstructure and motion, and progressively enriching the description with finer details. The complete multi-scale procedure is described in section 8.3.3.

Dynamic reconstruction

The dynamic reconstruction used in this article is an extension of the standard ART algorithm, and will follow the same structure as algorithm 6. Considering the inner "for" loop, at time t (and angle θ(t), the volume is warped with the measured displacement field

f (n-1) (x, θ) = f (n-1) (x + u(x, t)) (8.8)
(initially u(x, t) = 0). The computed projection of f (n-1) along θ(t) is compared with the recorded projection and the residual (i.e., their difference)

ρ (n) (r, θ(t)) = p(r, θ(t)) -R θ(t) f (n-1) (x, θ(t)) (8.9) is normalized and back-projected ∆ f (n) = B θ(t) [ρ (n) (r, θ(t))].
Finally the correction term is unwarped to the frame of the undeformed state,

∆f (n) (x) = ∆ f (n) (x -u
) so that it matches the reference configuration and it is added to the volume, f (n) = f (n-1) + ∆f (n) . Let us emphasize that theoretically, ∆f (n) should have been defined implicitly as obeying ∆f

(n) (x + u) = ∆ f (n) (x)
. The two expressions are equivalent only for small strains and rotations, otherwise the unwarping should involve the Eulerian rather than the Lagrangian displacement, and one can be computed from the other. Let us also note that for not too large strains and rotations, ignoring the difference between Eulerian and Lagrangian displacements simply slows down the convergence, but the final solution is not affected. In the present case, the choice was made to use the Eulerian registration to achieve the convergence for engineering strains as large as 20%.

A convergence criterion has to be chosen as in the ART procedure. Nevertheless, the criterion based on the functional value can not be used in this case as the reconstruction is unperfect. A convergence criterion based on the variation of the functional or a maximum number of iteration N DynART can be set. The procedure is described in Algorithm 7.

Algorithm 7 Proposed motion-corrected algebraic reconstruction procedure, DynART(p, u) Initialization n ← 1 Initialization f (1) ← 0 Choose a permutation, π, over N t acquisition times

for i ← 1 to N DynART do for t ← 1 to N t do τ ← π(t) f (n) (x, τ ) = f (n) (x + u(x, τ )) warp ρ(r, τ ) ← p(r, θ(τ )) -R θ(τ ) [ f (n) (x, τ )] ∆ f (n+1) (x) = B θ(τ ) [ρ(r, τ )] ∆f (n+1) (x) = ∆ f (n+1) (x -u(x, τ )) unwarp f (n+1) (x) = f (n) (x) + ∆f (n+1) (x)
Implement additional constraints on f (n+1) (e.g., positivity) n ← n + 1 end for end for As earlier mentioned, additional priors can be added in this procedure at the end of the inner "for" loop. In the following, only a positivity constraint for f is added at each iteration.

Motion identification

The full procedure is a staggered two-step process where alternatively the volume is reconstructed from a frozen displacement, and the motion is identified from a frozen estimate of the microstructure. The second step is described now.

At step l, the reconstructed volume, f l (x), although imperfect, is now considered as reliable. The projected residual fields ρ l (x, t) (computed at the end of the previous procedure when the volume is no more updated) contains patterns that are the signature of an incomplete motion correction. For the identification of the displacement field, the functional for a given f can be linearized around the previously identified displacement field

u l = u l-1 + δu δu = Argmin δv r,t R θ(t) [∇ fl (x, t)δv(x, t)] -ρ l (r, t) 2 (8.10)
For a better conditioning, the space and time dependencies of motion may be regularized, either using "weak regularization", with a penalty on spatial or temporal rapid variation of the displacement field to be added to the above cost function, or reverting to "strong regularization" by choosing a parametrization space composed of smooth functions of space and time. At this regularization step, any additional information pertaining to the experiment (e.g., synchronous measurements from sensors of different modalities such as force, pressure or temperature measurements, cardiac phase etc.) can be incorporated in the kinematic model through functional dependencies on such parameters. Qualitative features may also be incorporated, for instance, the sudden occurrence of a crack, may be accounted for by allowing a temporal discontinuity in concerned degrees of freedom for the kinematics.

The chosen reduced basis is composed respectively of N τ time functions, ϕ i (t), and N s vector spatial shape functions Φ j (x) such as

u(x, t) = Nτ i=1 Ns j=1 α ij ϕ i (t)Φ j (x) (8.11)
with α ij the time and space amplitudes that weight the basis functions. Setting ϕ i (0) = 0, the reference state is at initial time or angle θ = 0, u(x, 0) = 0.

The minimization of the functional with respect to the displacement parameters δα is performed using Newton's descent method. This procedure requires the computation of the advected image gradient and Hessian of Γ. They are built from the projected sensitivities

S ij (r, t) = ∂R θ(t) f (x, t) ∂α ij = ϕ i (t)R θ(t) [Φ j (x)∇ f (x, t)] (8.12)
Numerically, the sensitivities are computed from finite differences. The Hessian matrix and second member built from those sensitivities is

H ijkl = r,t S ij (

Multi-scale approach

If displacement magnitude is bounded by a length scale λ, one expects that the reconstruction is fair at a scale larger than λ. Hence, if the original image is convoluted with a Gaussian of width λ, it should well match its sinogram. One convenient property of the projection is that the projection of the convoluted image is the convolution of the original projection with a Gaussian of the same width. However, because of motion, this matching is not perfect but just fair. It means that one may estimate a better match by treating the deformation as a slight perturbation.

More precisely, the recorded projections are convoluted by the Gaussian of width λ,

pλ (r, θ) = r G λ (r )p(r -r , θ) (8.16)
where, G λ (r) = 1/(2πλ 2 ) exp(-|r| 2 /(2λ 2 )). The reconstructed volume from convoluted projections is written f λ (x).

Using the progressively identified displacement field, a more accurate determination of f can be achieved using the above described reconstruction. Because a large part of the displacement is expected to be captured in u, the idea is to repeat the above procedure but with a smaller gaussian filter, namely cutting down λ by a factor of two. Thus at each iteration, the displacement correction being smaller and smaller, convergence to the actual displacement field is expected. A convergence criterion is chosen on the norm of the residual variation or on the norm of the displacement correction.

The summary of the complete procedure is described in Algorithme 8.

i ← 1 to N t do for j ← 1 to N s do f (x, t) ← f λ (x + u(x, t)) Volume advection χ j (x, t) ← R θ(t) [Φ j (x)∇ f (x, t)] S ij (r, t) ← ϕ i (t)χ j (x, t)
Projected sensitivities end for end for

H ← r,t S ⊗ S Hessian b ← r,t ρS Second member α ← H -1 b
Motion identification end while

Test case

Two numerical test case are proposed to validate the procedure. In order to build the input data, two geometries are chosen, and two kinematics (one per case) deformed and projected at all considered angle. The obtained projections are then corrupted by a white Gaussian noise (standard deviation of 1% of the gray level dynamic of the projections and are used as the virtual experimental inputs for our procedure.

Both examples are carried out on 512 × 512 pixel images. The beam is parallel, and only N θ = 300 projections are acquired over a single 360°rotation.

• The first application corresponds to a moving Shepp-Logan phantom with large displacement magnitude (up to 37 pixels) and large engineering strains (27%). Large strains are chosen here in order to highlight the robustness of the proposed procedure as compared with previously studied examples where strains were about 1% [32].

• The second test is performed on a checkerboard with smaller displacements but a more complex time evolution composed of two separated modes: a steady drift superimposed to a high frequency pulsating motion.

In both test cases, the displacement bases chosen for the inverse problem were similar to the ones used for performing the direct problem, so that no additional model error (apart from noise) is introduced. The space functions Φ(x) are composed of four C4 mesh elements (4-node square elements with bilinear interpolations). The space basis N s is hence composed of 18 degrees of freedom.

Shepp-Logan phantom case

In this test case, the Shepp-Logan phantom is used and deformed up to 27%. For this test case, a single time evolution (linear drift in time) is applied. The imposed displacement field can be written

v(x, θ) = θ/N θ • V (x) (8.17)
The reference and deformed phantoms are shown figure 8.1. The maximum displacement amplitude (on the central node) is 37 pixels. The first reconstruction of the image (standard ART procedure), presented figure 8.2(a) is very blurry. Some parts of the phantom are split in two. The initial projected residual fields are very high and stresses that the reconstruction is not properly performed.

Before using the proposed procedure, the multi-scale procedure presented section 8.3.3 is applied to the projections to willingly blur the reconstruction. After 60 iterations (i.e., volume updates), the displacement field has converged. The corrected reconstructed volume is presented figure8.3. The edges are sharp and the gray level amplitudes are correct. The projected residual fields (true metric of our procedure) is mostly composed of the white Gaussian noise meaning that the proposed procedure has been successful.

The displacement error computed on the nodal values displays a standard deviation of 3.10 pixel. This result validates the procedure.

As a last validation of the phantom reconstruction quality, the reconstruction is compared to the reference volume f . It is shown in figure 8.4 that the reconstructed shape and positioning is very good. The final difference displays a "ghost" of the phantom that points out a small intensity error that does not appear in the residual fields.

Pulsating checkerboard case

This second test case is here based on a checkerboard composed of 8 × 8 squares of 35 × 35 pixels each. This square shaped pattern is chosen to exhibit reconstruction errors very clearly since sharp and straight boundaries are very easily detected, and hence the visual perception is a very severe test.

In this example, the imposed (supposed unknown) displacement field is composed of the sum of two parts • a pulsating motion: Temporally, a shifted cosine function (1 -cos(...)) (obeying the constraint of being null at time 0) evolution with a non-integer number of periods to avoid symmetry (here 2.35 periods during the full-rotation scan). Spatially, the displacement field is a centered dilatation/contraction.

• a linear drift in time for all nodes with random directions and amplitudes

The applied displacement field can be written

u(x, θ) = (1 -cos(2.35 • 2π • θ/N θ )) • V 1 (x) + θ/N θ • V 2 (x) (8.18)
The nodal displacement vectors V 1 (x) and V 2 (x) are shown in figure 8.5(a). The maximum amplitude of u is 30 pixels. The reference image, the deformed one at the end of the scan and the chosen C4 mesh are shown in figure 8.5. The maximum strain is about 25%. Because of the large motion amplitude, the initial reconstruction (i.e., obtained from a standard ART procedure for which u(x, θ) = 0) is fuzzy and the quality is very poor as can be judged from Figure 8.6. The projection of this blurred volume is compared with the initial projection to generate the initial projected residual fields ρ(r, t).

The first reconstruction without displacement field correction and the corresponding projected residual field are shown in figure 8.6. The initial residuals are very high. After 60 iterations, (i.e., 60 updates of the reconstruction) -performed in approximately 2 In order to correctly appreciate the quality of the achieved volume, a difference with the initial perfect one is shown figure 8.9. This difference highlights a perfect positioning of the reconstruction, and only slight discrepancies of the gray level intensity on the bright squares are visible.

Discussion and Conclusion

An innovative algorithm is presented to perform simultaneously a dynamic reconstruction of a moving sample with the identification of the full 2D space and time displacement field. The method is derived from Algebraic Reconstruction Techniques coupled with Projection based Digital Volume Correlation. The iterative algorithm is based on two steps

• for a given displacement field, a dynamic algebraic reconstruction algorithm is proposed. Each iteration of the procedure consists in comparing the acquired projection with the projected warped volume (deformed with current displacement field). The projected residual is backprojected, unwarped to match the reference space and added to the volume. • for a given reference volume, a P-DVC analysis allows the displacement field to be identified. The projection of the (unperfect) warped volume is compared with the acquired projections. The residual can be read as motion using the computed sensitivity fields. An update of the displacement field is then performed.

A multiscale procedure has been proposed as an essential ingredient to properly correct large displacements. The acquired projections are first convoluted with a Gaussian kernel of large width (low pass filter) to increase its correlation length and capture large corrections from the linearized P-DVC functional. The Gaussian filter is then progressively reduced, following the residual norm evolution, to identify finer details.

The post treatment procedure, that exploits the same data as a standard acquisition (same number of projections and standard projection operator), has been tested with two challenging numerical examples (with large displacements and strains). The first is a Shepp Logan phantom with large displacement fields (up to 1/4 of the phantom length). The second is a checkerboard with a pulsating motion in time. Both examples are corrupted by a white Gaussian noise that probes the robustness with respect to the acquisition noise. The two applications show a nearly perfect identification of the displacement field and dynamic reconstruction. Performed with a parallel projection algorithm for simplicity, the exact same method can be applied with any projection model.

The proposed dynamic reconstruction algorithm has been devised as an extension of the ART algorithms. It is convenient with those approaches to include in the process an a priori knowledge of the scanned specimen (assumption on the gray levels, its variations, the number of phases, its sizes, etc.). Many different regularization have been proposed in the literature that enable to obtain high quality reconstructions, with less artifacts, from less projections or missing angles, etc.. Because those regularizations are independent of the current algorithms, it was chosen not to implement them and focus on the proposed method performances without any 'additional help'. Nevertheless, they are fully compatible with the proposed approach and can be implemented in a transparent fashion. When aiming to perform ultra-fast acquisitions with few angles, they would certainly be very precious to accelerate convergence, and improve reconstruction quality.

In the proposed examples, the optical flow was kept constant. Some applications may require to include a gray level variation model. A perspective of this work could be the scan of in situ mechanical test with high strains, the identified deformation could be used to correct for absorption evolution of the material considering a constant beam intensity.

The proposed procedure shows performances that can be beneficial to numerous fields. The clear reconstruction of the moving sample allows for qualitative and quantitative analyses:

• combined with Digital Volume Correlation [38] between well reconstructed volumes • combined with image segmentation for diagnosis from radiology

• combined with ultra-fast tomography acquisition as recently available from some synchrotron beam-lines [39,6]. This is key for data assimilation [40] and model identification and validation in material science [41] with CT-scan as already developed with MRI [42].

Chapter 9

4D dynamic CT

This chapter is an ongoing work deliberately concise that aims at developing a complete dynamic reconstruction framework applied to medical images. This paper proposes a 4D dynamic tomography framework that allows a moving sample to be imaged in a tomograph as well as the associated space-time kinematics to be measured with nothing more than a standard scan acquisition. The method exploits the consistency of the projection/reconstruction operations through a multi-scale procedure. The iterative procedure is split into two coupled parts: a motion compensated reconstruction algorithm and a projection-based measurement procedure that reads the displacement field directly on each individual projection. The procedure is tested on the moving tomography of a breathing chest, generated synthetically from static 10 real CT-scans of a breathing patient. At convergence, the initially blurry reconstructed volume has been cleaned from motion artifacts resulting in an improved reconstruction quality showing sharper edges and finer details.

Introduction

Being able to capture the inner structure of a sample in a non destructive way is of utmost importance for medical diagnoses or for quantitative measurements. Among many different 3D imaging techniques (Magnetic Resonance Imaging (MRI), Ultrasound imaging, ...), tomography is the most frequently used. Initially developed for medical imaging [1,2] as it enables distinguishing the different soft and hard biological tissues , tomography is now widely used in many other fields (e.g., biology [3], material science [4,5,6], etc.). Performed with different instruments (e.g., X-rays medical scanners, synchrotron beams, lab-CT, etc.) and waves (mostly X-rays but also neutron, optical, positron, etc.), a wide range of imaging (space and time) resolutions and contrast suited to almost any material is accessible.

A tomography scan is performed as follows, the specimen rotates over 180 • or 360 • with respect to the source-detector pair (often a rotating sample in material science and rotating source/detector in medical imaging) and a series of radiographs at distributed angles are acquired and normalized into projections. The collection of projections at all angles constitutes a so-called sinogram. From the latter, inverse reconstruction algorithms [7] have been developed to reconstruct the 3D imaged volume. The technique relies on the strict satisfaction of conditions, in particular concerning the geometry of the set-up and the motion of the sample as a perfectly controlled rigid rotation with prescribed axis and angles.

Motion in tomography is a complex topic that offers the two sides of the same coin:

• If the motion happens during the imaging process, it leads to a poor quality reconstruction with blurry edges [8]. The reconstruction can hardly be exploited, even visually. This displacement may come from various reasons: motion of the sample itself (e.g., due to the experiment itself, creep, heartbeat and breathing [9], etc.), unexpected motion of the setup, uncalibrated tomograph (although this is not exactly motion, it can be translated to such an artifact), or related to quick phase and texture changes [10,11]. If the motion is not too fast with respect to the acquisition time (sub-second scale in a synchrotron, minutes/hours in a lab µCT), the experiment can be performed continuously. Otherwise different solutions have been found to reduce the scanning motion: using smart sampling in time for periodic motions (e.g., based on phase measurements [12,13]), accelerating the scanning procedure (e.g., faster setup, brighter sources [14], multi-sources [15,16], developing ultrafast scanning procedure in the synchrotron [17,18]), using, in medical imaging, sedative, beta-blocker to reduce the heartbeat rate (especially for babies or small children [19]), requesting the patient not to move, apnea etc.

• On the contrary, if this motion is measured and controlled, it gives access to extremely rich quantitative data that are precious pieces of information for diagnoses, model identification and data assimilation (e.g., in material science [5,20,21], biomechanics [22,23,24]). Finding a way to image the microstructure and not be polluted bu the displacements during the scan is a challenge. The solution is often to image the sample at different static states (i.e., performing a 4D (space/time) measurement procedure [21,25,9]). Spatio-temporal deformable registration has received considerable attention in literature, in material science for the identification and validation of mechanical models (e.g., [26,21]) and in medical imaging, mostly in cardiac image analysis [27,28], but more recently also for respiratory-correlated imaging of the thorax [25,9]. Those space-time procedures take time and the scanned samples are subjected to huge quantities of radiation doses. Those techniques lead to highly redundant acquisitions as the same microstructure is imaged multiple times (only its deformation evolves but its complexity is generally low compared to the microstructure).

Although the sample is imaged with motion or an uncalibrated tomograph, different methods have been developed in the literature to get a satisfactory reconstruction quality. The methods are often composed in two successive procedures: (P1): evaluate the motion (or calibration parameters), based on either additional sensors, or full field measurement methods on the acquired images (e.g., image registration techniques such as digital image/volume correlation [29]).

(P2): correct the reconstruction procedure from the results.

For the calibration of tomograph parameters (i.e., positions of the detector, source and sample at any time), online procedures [30,31,32,33,34] use the scanned moving sample itself to retrieve all the calibration parameters. Those measurements are often based on affine transforms. The obtained corrected parameters improve the quality of the specimen reconstruction [35].

When the sample is subjected to a more complex motion, the reconstruction can not be corrected by only hanging the calibration parameters. After having measured the displacement field, the inverse reconstruction algorithm has to be enriched with motion. Techniques using a motion-compensated Filtered Backprojection (FPB) algorithms have been developed [36,37,38]. Each filtered back-projection is warped with the timeassociated measured space-time displacement field so that the sum of all those warped back-projections gives access to the volume in a static configuration.

The topic of dynamic tomography is currently attracting a lot of interest. Three very recent papers tackle the CT imaging of moving samples. In [39], the authors studied a burning match stick imaged at the PSI synchrotron. With 18800 acquired fast projections (in 18.8 s), many 3D motion-blurred volumes could be reconstructed. When applying optical flow analysis, the measured displacement between the volume at time t -1 and t + 1 is linearly interpolated at time t and thereby allowing the central volume at time tt to be corrected. A full 4D space-time 'movie' of the phenomenon could be obtained. This method however requires the use of many acquired 3D volumes (at least 3) thus a high dose. The measured displacement has a relatively small amplitude compared to the volume texture characteristic scale. A criterion based on the image reconstruction quality is not easy to set and the quality has to be appreciated visually. A similar recent technique [40] deals with the correction of a volume using Digital Volume Correlation and an extended Simultaneous Algebraic Reconstruction Technique (SART) algorithm. To be able to correct a single rotation volume, the authors sub-sampled the acquired projections in 2 sub-acquisitions from which the motion is evaluated and further involved in the reconstruction strategy. This method is however not suited to large and irregular displacements. In [41], the volume sub-sampling is performed more easily because of a specially designed sampling acquisition strategy with many rotations (that cannot be adapted to any tomography). One displacement field, constant in time, is estimated each successive pair of reconstructed volumes and is used to correct the reconstruction procedure. In this latter reference, although the tackled displacements and deformations are important during the entire test, an important quantity of reconstructions could be performed thanks to the acquisition procedure leading to small incremental displacements.

In all those presented very recent developments, the displacement field measurement (carried out from 3D registration methods) and corrected reconstruction can hardly be performed with a single 360 • acquisition if the motion is important and irregular. Moreover, the final quality has to be appreciated qualitatively and no objective criteria on the reconstruction quality or on the motion estimation are proposed.

A recently developed full field measurement procedure consists in reading the displacement field of the sample directly from the acquired projections. The projection-based measurement is very well fitted to the above mentioned first point (P1) dynamic tomography as it provides a space-time displacement for each frame [42]. Developed in the medical imaging field and called 2D/3D image registration [43,44,45,46,47], a similar approach has also been developed in material science and called Projection-based Digital Volume Correlation, P-DVC. In those approaches [48,49,50, ?], an initially well reconstructed volume is warped in space and time such that its projections match the recorded deformed projections. This procedure offers huge savings in acquisition time because the experiment is followed by few projections (only two projections per state in [48,49,50] and a single one in [33]) instead of complete volumes.

It is to be noted that the projection-based measurement functional is very similar to the one used in standard algebraic reconstruction techniques (ART). A coupling between ART and P-DVC appear thus natural. It is proposed to combine both of them to achieve the desired dynamic tomography. An initially motion-blurred reconstructed volume is used as the reference volume and leads to an evaluation of the displacement field. The measured motion is then included in a motion-compensated ART algorithm and these two steps are performed iteratively. In [33,34], the calibration of parameters of the tomograph was estimated using affine transforms on the sample. In [51], the coupled measurement/compensated-reconstruction framework was developed with 2D synthetic examples of moving phantom with known simple motions. Those synthetic test cases resulted in high quality reconstructions as well as the full space-time kinematics to be identified.

We propose in the present paper a fully coupled projection-based measurement method with a motion-compensated ART algorithm embedded in a multi-scale framework. The method is applied on the CT-scan of the chest of a breathing patient. The proposed method provides the 4D regularized motion of the lungs along the 4 respiratory cycles and the un-blurry reconstruction to be carried out.

Coupled reconstruction-measurement problem

The collected intensity for each detector position r = [r, z] (where z is parallel to the specimen rotation axis, and r is perpendicular to it) and rotation angle is to be related to the relative beam intensity attenuation. The Beer-Lambert law relates the line integral of the material absorption along the X-ray path L(r) from source to detector at position r, to the recorded intensity I(r, t) at time t (and rotation angle θ(t)):

I(r, t) = I 0 (r, t) exp x∈L(r) µ(x, t)dx (9.1)
with µ(x, t) the linear attenuation coefficient and I 0 (r, t), the intensity recorded without the sample on the beamline called often flatfield or white field. To obtain the projections, one has to normalize the recorded intensity with the flatfields (after possibly subtracting off darkfields acquired in the absence of beam), and compute the co-logarithm, as written in (9.2)

p(r, t) = -log[I(r, t)/I 0 (r, t)] (9.2)
Possibly standard tomography artifacts [52,53] (e.g., beam hardening, intensity variations, etc) can be corrected at this stage. The collection of N θ projections p(r, θ(t)) for each angle is called the sinogram. The problem can be written in its discrete form as a simple linear system. In the present work, the projection operator as defined in [54] is performed using the ASTRA toolbox and in particular, its GPU implementation. When motion u(x, t) is considered during the experiment, the projection of the deformed volume, defined for every voxel

x = [x, y, z] is written p(r, t) = Π θ [f (x + u(x, t))ψ(u(x, t))] (9.3)
with ψ(u(x, t)) = 1 -Tr(E) the local volume variation, computed from the Green-Lagrange deformation tensor E. This correction term is used to modify the gray level with respect to the absorption (thus µ(x, t)) changes of the material. The functional that can be defined is then, based on a least square minimization

Γ(f, u) = r,t (Π θ [f (x + u(x, t))ψ(u(x, t))] -p(r, t)) 2 (9.4)
This functional exploits the so called consistency equation [34] and provides a validation metric of the procedure (e.g., quality of the reconstruction and measured displacement field): the residual field ρ(r, t) = Π θ [f (x + u(x, t))ψ(u(x, t))] -p(r, t). The general coupled problem is split to two minimizations and solved in a fixed point framework, iteratively such that

f = Argmin φ Γ(φ, u) (9.5) u = Argmin v Γ(f, v) (9.6)
Figure 9.1 is a scheme of the staggered two step procedure. The framework is separated in two iterative procedures: the measurement of the motion from a reconstructed volume and the corrected reconstruction. At the beginning, the motion is initialized at 0.

Reconstruction with motion

The first part of the two step procedure consists in reconstructing the 3D microstructure of the sample from the acquired sinogram and the measured displacement field (initialized at 0 unless additional sensors information or initial guess can be exploited). The use of an ART reconstruction algorithm is natural as it results from equation (9.5).

It was chosen here to use a SART algorithm. At each iteration, indexed by (n) the volume is warped with the measured displacement field and compared with the recorded projection, the difference of which defines the residual field that generates the projected residual fields

ρ (n) (r, t) = p(r, t) -Π θ [f (n-1) (x + u(x, t), t)ψ(u(x, t))] (9.7)
which is normalized and back-projected ∆f

(n) = B θ [ρ (n) (r, t)
] with the backprojection operator B θ . Finally the correction term is unwarped to the frame of the undeformed state, so that it matches the reference configuration and it is added to the volume,

f (n) = f (n-1) + ∆f (n) (x -u).
A convergence criterion has to be chosen as in the ART procedure. However, functional value can not be used directly in this case as the reconstruction is imperfect. A convergence criterion based on the incremental variation of the functional (or an escape condition limiting the maximum number of iterations) can be set.

Projection based measurements

The second part of the two step procedure is the projection-based displacement field measurement based on Projection-based-DVC. This procedure uses as reference volume the previously reconstructed volume and minimizes, with respect to the displacement field, the difference between the projection of this warped volume with the recorded projections. DVC algorithms are known to be ill-posed [20], consequently the displacement field has to be regularized. Among many possible choices to regularize the displacement field, our choice is to write the displacement field on a reduced kinematic basis composed of few degrees of freedom u

(x, t) = ij a ij φ i (x)σ j (t) (9.8)
with a ij unknown amplitudes associated with φ i (x) and σ j (t) respectively space and time functions that characterizes the space-time evolution. Those fields may come from a prior knowledge on the kinematics. The space can be expressed on subsets (called local DIC [55]), on a mesh kinematics (called global DIC [57]), on b-spline interpolations, etc. The time can be, for example, imposed by additional sensors measuring the force measurements (used in most of mechanical experiments) or a phase signal (accessible in cardiac imaging). Finally, both space and time evolutions can be coupled using a model-based kinematics (called integrated DIC [58]) or using knowledge from prior similar experiments then reduced with singular value decomposition [59].

The choice of the two space and time functions φ i (x) and σ j (t) and their numbers are very important. Too many degrees of freedom and the system becomes ill conditioned and too few introduces model errors. The residual fields of this procedure ρ, as defined previously, inform on the quality of the procedure and on the chosen model. From those residual fields, the model and procedure can be (un)validated and if required, the model can be enriched with additional degrees of freedom.

The non-linear problem is solved using a Newton-descent algorithm, i.e., successive linearizations around the found solution, until a convergence criterion is reached. All components of the gradient of the cost function with respect to all degrees of freedom are computed using finite differences. Because of a small number of degrees of freedom, each linearized problem can be easily inverted (contrary with other convergence procedures based on multiple research directions on separated spaces [60])

Multi-scale procedure

The key to be able to solve the previous two step framework, even with large displacements, is to apply a multi-scale coarse-to-fine procedure. Indeed with large displacements, the above linearized functional will not be able to converge. The measurement procedure has first to rely on the global shape of the scanned sample (e.g., the low spatial frequency patterns) and progressively refine on the fine-scale detail texture. For that, a pyramidal approach (see figure 9.2), starting with the coarse scales and ending with the finest is used with two parallel enrichments

• the projections are smoothed from the convolution with a Gaussian kernel whose characteristic length can be adapted to the displacement amplitude. The high frequency patterns of the projection and reconstruction are erased and the functional is hence driven by the global shape motion. The projections can also be downsized, in space and time to accelerate the procedure as downsampled images carry the same information as the filtered one.

• the kinematic model can also start from a very coarse kinematics (e.g., simple rigid body motions, affine transforms) and progressively be enriched while the microstructure appears.

For each scale, the procedure (with a chosen mesh and volume resolution) is computed until convergence of the residual fields and displacement variation. In the following application, 3 scales in images and 2 in models are used (the two last volume scales are computed with the fine mesh).

Results

The application is performed on the CT scan of a breathing patient. The data are generated synthetically from 10 real reconstructed CT-scans [9,12] of the respiratory phases from the CREATIS laboratory is projected 8 times successively at 8 angles (and then this projection operation is iterated) in order to be correctly sampled in time without too many respiratory cycles (that would not be realistic). The obtained sinogram is hence the input of our procedure. The projection operator used to generate the data is a parallel beam for simplicity issues and performed with the ASTRA toolbox [54]. Although common medical CT beams can be modeled by a fan beam, it only affects the projection Π θ and does not modify the principle of the current procedure.

4D kinematic measurements

The time evolution can be based on external sensors (force measurements, breathing signal, etc). It was proposed in this application not to include any additional information.

Because the bottom part of the scan moves with the breathing, the mean intensity of the image, computed over a region that contains mostly the bottom part of the image, varies and is assumed, at first order, to be proportional to the motion. This measured signal is hence used to regularize the time evolution. In addition, the absolute value and square of this signal and a 1 period sinusoidal function (that could be related in real cases to the scan rotation issues) are added to the time library. The space kinematic model is embedded in a multi-scale framework. The kinematics starts for the first kinematic scale with a simple C8 element (i.e., 8 noded-cube elements with tri-linear interpolations) composed of 24 degrees of freedom (i.e., a total of 96 dofs for a ij for the 4 time steps). Then, when the kinematics allows a better reconstruction of the volume, a finer mesh composed of 75 T4 (i.e., 4 noded-thetraedra elements with tri-linear interpolation) is used creating 900 degrees of freedom. In the application, the volume has 3 scales (downsizing of 4 -2 -1) and are associated with mesh (coarse C8fine T4 -fine T4). Each scale is computed until convergence (4-5 iterations)

The measured displacements are shown figure 9.4. The vertical motion is important around the diaphragm and almost 0 at the top part of the chest. A maximum vertical displacement amplitude is measured at approximately 7 voxels and corresponds to 14 mm.

Reconstruction

Three image scales are used where the projections are filtered and downsized by a factor 4, 2 and 1 (initial projections). The time sampling is also downsampled by the same factor. The fine reconstruction is performed with all the 360 angles. No additional regularizations are added in the procedure. Slices of the reconstructed volume before and after are shown in figure 9.5. The correction naturally focuses on the large displacements area, i.e., the bottom of the chest around the diaphragm. It can be noted that the diaphragm is much better reconstructed after the displacement correction. Initially smooth and blurry, the improved edges become much sharper. The bottom plot shows two intensity curves corresponding to the two yellow lines in the zoomed areas. After the displacement field correction, the diaphragm is clearly visible while the initial curve is much more fuzzy. Segmenting the volume for example would be easier with the corrected volume.

Objective quality criterion

The true metric of our procedure is the norm of the projected residual field. This residual field is expected to be null at convergence (see figure 9.6 before and after the procedure correction). The residual field contains what was not captured by the kinematic correction procedure thus can be composed of convergence issues, noise and artifacts but includes also all model errors (approximate kinematic model, inaccurate interpolations, etc). It can Figure 9.5: Sagittal and coronal slices of the reconstructed chest without (left) and with (right) motion correction (the two red lines indicates the slice positions). It can be observed (especially in the zoomed area) that the diaphragm is better reconstructed and sharper (as pointed out by the two white arrows). On the sagittal view, a horizontal plane that polluted the reconstruction disappeared (white arrow). The bottom plot shows 2 intensity curves extracted from the left and right volume on the zoomed area (yellow lines). The dashed curve, corresponding to the non-corrected volume, is smooth. The solid curve is sharper because of the displacement correction. be seen that the residual values are very high in the beginning. On the bottom part, large positive and negative values are visible resulting from the oscillating respiratory motion. At the end of the procedure, most of the oscillating motion has been erased meaning that the estimated displacement field was trustworthy. The norm of these residual fields decreased from 3.2 down to 1.1 at the end of the procedure. Because the projection of the warped volume corresponds to the recorded deformed projections, the procedure is deemed validated.

Discussion

A 4D displacement field measurement coupled with a motion compensated reconstruction procedure framework has been developed. The method is based on nothing more than a standard sinogram of a moving sample. The procedure consists in two parts: (i) a motion compensated procedure extended from a SART algorithm and (ii) a displacement field measurement from the reconstructed volume and the initial sinogram.

The procedure is tested from an artificial time series composed of projections issued from 10 real volumes acquired at different phases og the respiratory motion. Mixing all states together, a moving sinogram is synthetically generated whose reconstruction leads to motion blur artifacts. The kinematics is regularized with 4 time functions and a 75 noded T4 unstructured mesh creating a 900 degrees of freedom problem. The multi-scale procedure, performed for the volume and model, results in a fast convergence. Finally the quality-improved reconstruction was cleaned from a large part of the motion blur artifacts leading to projected residual fields of very small amplitude.

This procedure opens up new avenues in dynamic tomography enabling to image moving samples, even if the displacement is large and without additional data as compared to a standard measurement. Moreover, any additional sensors and detectors could be used transparently in the above method to regularize, either the volume reconstruction (e.g., assumptions on the phases [61,62], TV regularizations, dictionnary of commonly encountered tissue morphology, etc.) or on the kinematic model (e.g., force and displacement sensors, etc.). Using dual sources scanners is also a promising perspective, especially when the sources are orthogonal [15,16] so the displacement field could be completely estimated for each acquisition time (i.e., a better conditioned problem).

Bibliography

Conclusions and perspectives

The common thread of this manuscript is the data reduction in the initial Radiographies -Volume reconstruction -Measurement -Parameter identification sequence presented in the Introduction. Being able to short-cut the standard reading sequence with projectionbased measurements opens up new opportunities when designing a mechanical (or more general) test. The quantity of required data has been drastically decreased with the successive developments (leading to less acquisition time and reduced doses). This PhD dissertation is composed of two parts: different developments of the projection-based identification when a reference volume is either known (Part I) or unknown (Part II).

The first part, comparable with stereo-setups, offers gains of more than two orders of magnitudes in acquisition times and doses. From multi-view identification in Chapter 2, the procedure has been developed with a single acquisition per state, carried out continuously with simultaneous loading and rotation of the sample. Fast behaviors (time constant longer than 1-2 s) can be measured in a lab-CT (Chapter 3-5).

The second part aims at identifying both the kinematics and microstructure. The method called dynamic tomography is a coupling between P-DVC methods using a single projection and a motion-compensated reconstruction procedure. These extremely promising developments have a huge potential in materials and medical imaging.

In figure 10.1, the diagram proposed by E. Maire and P.J. Withers [1] has been enriched 209 with 4 experiments reported in this dissertation (blue crosses) and one that is only based on a virtual (synthetic) test (yellow cross). Few remarks are to be made: (1) The points do not take into account the reference acquisition (required in the development of the first part). ( 2) The target of the previous work was focused on the methodology and development and not on the acquisition time (e.g., acquisition of 8 projections instead of 2 in [2], average of 5 frames in [3,4] even if a single one could have been used, etc.). ( 3) The dynamic tomography of the medical CT has not been added as the data were re-generated synthetically and because the helicoidal acquisition is a different experimental technique.

(4) Space resolution has never been a goal with P-DVC procedures. ( 5) Even if the 3D volume is not imaged in the deformed state, it could be computed using the reference one, the displacement field and possible gray level variations. The values are summarized in the Table 10.1. The number ( 5) is not related to a real experiment. However it has a huge potential when dealing with periodic motions and has been tentatively included in the diagram. The first part of the manuscript is dedicated to the development of Projection-based approaches with the use of a reference volume. The target of the successive chapters is to reduce the quantity of data required in the identification procedure (kinematics, materials and fluid shape). Initially performed with 2 projections per state [7,8,9,10] (multi-view framework), a single projection has been shown to be sufficient thanks to time regularization, offering the opportunity of a continuous evolution (i.e., continuous acquisition, loading and rotation). Applied to different experiments (compression with cracks, tension with non-linearities (plasticity), imbibition, vibration) and materials (plaster, cast iron, sandstone, synthetic materials), the developments lead to huge gains in required data and acquisition time.

• Chapter 2: extracted from [2], P-DVC was extended with cone beam in a lab-CT with 2 projections per time state. The application of the method was performed on a DCDC plaster sample, filled with copper powder to enhance the contrast. With a 14 degrees of freedom kinematics: 12 boundary conditions of a Finite Element elastic model and 2 crack positions, the measurement was performed on 13 loading steps with and without cracks. The measured displacements were used in a weighted FEMU-UF identification procedure to extract 2 material parameters: E and K IC that turned out to be in accordance with ex-situ measurement values. The very low level of residuals at the end of the procedure validates the method (some detector artifacts could even be seen in the projections meaning that most of the kinematics residuals were corrected). In this experiment, a gain of 350 in saved data was reached compared with standard DVC requirements. Nevertheless, few points were limiting and were improved in latter works: (i) the acquisition of 8 radiographs when only 2 are required, (ii) the experiment was stopped at each acquisition, (iii) an average of 50 frames was performed to increase the signal to noise ratio (unnecessary because the procedure is not sensitive to uncorrelated noise).

• Chapter 3: extracted from [3], P-DVC was extended to a regularized framework where the identification could be performed using a single projection per time step. The development of time regularization allows imposing both space and time evolutions. As in space where it is very natural to couple all individual pixel/voxel information (e.g., in a Global framework), coupling time steps enables a global resolution of the problem. This means that all (or at least more than 1) time steps / projection angles contribute in the identification of the reduced time basis. For example with the use of polynomial time evolutions, even if some degrees of freedom have no sensitivities in a direction of the space, the previous and the next direction will allow solving the problem. The method was applied to the measurement of the kinematics for a slender sample. The kinematics basis was chosen to be vertical sections animated by rigid body motions, in the spirit of the beam theory (although a global FE mesh could have been used as well). The identification was performed using a PGD algorithm: the displacement field was written as a sum of space/time modes which are identified successively. This modal identification naturally starts with the most salient to the least sensitive modes and thus can be stopped when a residual convergence criterion is reached. The application is performed with two steps measurements: a rigid body motion measurement that cleaned up a large part of the residuals and the deformation field. The results show the large top and bottom displacements of the sample and the central non-linearities (plastic deformations).

Although only the vertical displacement was plotted, the horizontal displacements (bending, torsion, etc.) was also captured. The signal to noise ratio (thus residual fields) evolution with the successive displacement field corrections validates the procedure. Because a single projection per state was required, the experiment was performed with a continuous load rotation and acquisition. This continuous method is thus not affected by relaxation or creep that happened in similar experiments performed for sequence volumes measurements. This 5-minute P-DVC experiment which consisted of 127 loading steps allowed gains in the experimental time of more than two orders of magnitude.

• Chapter 4: extracted from [4], P-DVC is extended to integrated approaches. Based on the same experiment and full field measurement results of Chapter 2, two identification procedures of a model have been developed: FEMU and Integrated approaches. The chosen model was a simple 1D free-shape plastic beam model driven by 8 control points in the (ε zz ,σ zz ) space. The FEMU procedure uses the measured displacements to identify the behavior. The small displacement residual at the end of the identification allowed the model to be validated. The integrated approach identified the 8 unknowns directly from the projections. The parameter variations that have a signature in motion are here 'decorated' with the microstructure and projected in order to generate the gray level sensitivities of the constitutive parameters. With this approach, the time evolution is driven by the force measurement and the variation of the control points modify the kinematics. The short-circuit of the entire reading sequence is here completed. The residual fields at the end of the procedure are in the same range as obtained with the kinematic measurements showing that the model was both relevant and well identified.

• Chapter 5: extracted from the work performed in Lund University, Sweden [6], in collaboration with Pr. S.Hall and his team, the projection based measurement was here extended to gray level variation tracking. The fluid invasion on a porous sandstone sample was imaged using neutron tomography at 5 fps. With neutron, the high fluid absorption is visible on the projections and evolves from the bottom of the sample to the top. The residual fields could thus not be corrected with the gradient of the microstructure. In order to model the phenomenon, the design of a fluid column was performed based on a time evolution (low order polynomials), the 3D height of the front (low order bivariate polynomials) and a front saturation profile.

The three coupled problems were solved using a PGD framework. Two examples were treated: (i) measurement of the front during a single 3D scan (i.e., radiographs acquired at 5 Hz) and (ii) during the entire experiment with the fluid starting from the top to the bottom of the sample. In both examples, the measurements allowed the projected residual field to be cleared. The identification performed 300 times faster than the 3D methods allowed small variations of the front (less than 1 voxel) to be captured. A comparison with 3D analysis showed very good agreement.

This work is very different from the previous developments as the correction model is not restricted to the advection of an initial microstructure but could be designed with a gray level variation correction model. Being able to modify the 3D gray level is interesting for example when dealing with phase transition, cracks voids or cavity or bubble nucleation, etc. It is thus to be noted that the reference volume is performed at the end and only used to give an initial fluid column intensity. This 3D scan was not necessary and everything could have been performed without reference volume. (and this chapter transfered to Part2).

• Chapter 6: is a proof of concept to perform an in-situ vibration modal measurement in a tomograph. Because an appropriate sampling in time is not possible (neither with volumes nor with radiographs), this goal is challenging. However such sampling is not required. Because the motion is periodic, every projections acquired at random times are linked with the modal basis that could be retrieved with a large number of projections. The projections should also be acquired at different angles to have sensitivity in the different space directions. The key lies in being able to give a link to bind the different angles together. A statistical assumption on a representative sampling per projection angle of a steady state distribution is made and allows the procedure to work. The measurement of the modal basis and the random amplitudes is performed on a synthetic test case with 3 projection angles and 150 frames per angle. The results validate the developed measurement procedure. Because all projection angles may not be equally sensitive to the displacement (e.g., the displacement of one mode may be predominantely in the direction of the projection), the measured amplitude of each mode may not be equally accurate but this leads to no detrimental consequence on the identification of spacial vibration modes.

Dynamic tomography

The second part of the manuscript is the development of dynamic tomography. For phenomena that can not be stopped (e.g., heart beating, respiratory motion, material melting etc.), the P-DVC procedure described in the previous chapters can not work as the reference image can not be acquired. However, because it is possible to measure the displacement field from projections, this motion can be used to correct the reconstruction procedure. The motion blur generated by a static reconstruction procedure is cancelled and the reconstruction can be made much sharper.

• Chapter 7: extracted from [5], the online calibration of the tomograph using P-DVC techniques is performed. The application is extracted from a 4D in-situ test performed by Amine Bouterf and Ante Buljac. During the night the sample is maintained loaded in the tomograph. Because of uncontrolled motion, thermal variations etc., the calibration of the tomograph changed leading to a blurry reconstruction. A P-DVC procedure is performed between the blurry reconstructed volume and its own radiographs highlighting the inconsistencies in the reconstruction procedure (i.e., the motion or wrong calibration parameters). With the space regularized using rigid body motions and a time regularized using sinusoidal functions, a PGD procedure is implemented to capture the first motion modes. Then after few modes identification, the reconstruction is updated with the obtained calibration. The updated volume is of good quality with sharp edges. The procedure, that does not require more than a standard acquisition, is thus validated.

• Chapter 8: extracted from [11]. A dynamic tomography method is developed in 2D with synthetic phantoms. The previously developed method for motion identification from blurry volumes is extended in a framework where the identified motion is used simultaneously in the reconstruction procedure. A motion compensated SART procedure is hence coupled with P-DVC such that every identified displacement field allows the reconstruction to be sharper. The method is applied on 2D examples using a parallel projection. The first test is a Shepp-Logan phantom with large displacements and strains but a simple linear time evolution. The second test case is a pulsating checkerboard with two temporal modes. In both examples, in which noise has been added, from an initially blurry reconstruction with large projected errors, the projected residuals are completely erased. Moreover the comparison of the reconstructed images with the ground truth shows that the reconstruction is almost perfect and well positioned and thus the approach is deemed validated. At the end, the full 3D (2D space + time) kinematics and high quality microstructure are identified altogether.

• Chapter 9: The developed dynamic tomography is extended to 4D for medical imaging. A collection of 10 medical CT-scans of the different respiratory phases (downloaded from the CREATIS website 1 ) are used to generate a dynamic test case with a patient breathing during the acquisition process. The method does not require more than a standard sinogram acquisition. The dynamic tomography framework is developed with gray level corrections related to the volume (thus X-ray absorption) variation and carried out in a multi-scale procedure starting from coarse scales (i.e., downscaled volume and reduced model) to the finest (i.e., full resolution volume with a complex T4 mesh). The clear improvement of the reconstruction quality is visible especially in the large displacement areas around the diaphragm. The edge of the diaphragm/lung fronter is sharper. The evolution of projected residual fields (corresponding to a residual sinogram) highlights the accuracy of the method because a large part of the oscillating motion error are corrected.

Regularization vs Freedom: One difficulty in the measurement procedure is to select an appropriate kinematic basis. Too many degrees of freedom and the system will be ill conditioned thus very difficult to solve and very long to compute using P-DVC. Not enough freedom and it will introduce model errors that obviously pollute the entire identification. The correct balance is not easy to find. Different strategies can be developed:

• Penalizations: different penalizations can be introduced in the space and temporal fields. The easiest to implement is the Tikhonov regularization [12] although the amplitude factor may not be easy to tune. Local elastic constraints [13] can be used instead and are more physical. Penalizing large velocities when dealing with a space/time problem can also be implemented.

• Models (i.e., integrated approaches [14,15]): with integrated methods, the regularization is performed based on a physical model driven by very few parameters. With a relevant model, the displacement field should match the observed kinematics. The difficulty is to select a relevant model.

• Trial-error: based on the residual field with localized error, the user can update the model and add progressively degrees of freedom [16]. This method is usually performed with crack initiation and propagation. In such cases the mesh can be opened, nodes can be added or removed, etc.

Sample texture: The sensitivity to the parameters to be identified is the key parameter that has to be carefully considered. Although we aim at reducing the number of experimental data, we do not want to loose too much sensitivity. Sensitivity can thus be tuned prior to the experiment. An optimization of the test would allow selecting specific angles, loading rate, geometry, etc. Concerning the sample itself, it is possible to enhance the sensitivity using tracers. These tracers need to have a different X-ray absorption than the specimen to be visible. They can be included in the sample whenever possible, as speckles (e.g., absorbing powder in melted specimen as in the plaster test case, medical tracers, etc.), by modifying the geometry and shape of the sample (adding holes, notches, using topology optimization, etc.) or applied on the surface (e.g., a simple white paint composed of absorbing titanium dioxide, particles or aggregates). The two latter cases would consist in a P-DVC based on the 3D surfaces of the specimen. Working with X-ray untextured volume with (sensitive) edges could hence be possible. It is noteworthy that the important term is, in P-DVC procedures, the projected sensitivity Π θ [∇f (x)ψ(x)] and not the 3D sensitivity. A well textured sample with very small and numerous particles (say 0.5% of the cross section size) could be (depending on the kinematics ψ(x) of course) considered as a uniform projected gradient thus without projected sensitivity.

Imaging procedure: The P-DVC procedure is based on selected projections in order to identify parameters of a model. Many parameters could be optimized to enhance sensitivity: the choice and number of angles, the rotation speed, the loading speed.

All projections are not useful for the model identification. In the DCDC test case for example, the projection facing the hole was the most important as this view was highly sensitive to the crack opening. Although, the second projection angle had no sensitivity to the crack, it enabled the whole kinematics to be captured. Being able to select the most meaningful projections is a challenge when considering many parameters. A scalar criterion has to be selected for the optimization procedure. Another strategy could consist in selecting successively projections focused on the minimum parameter sensitivity until all values are above the noise uncertainty.

For the optimization of the sample shape for the sensitivity to a multi-parameters model (that is another topic), a criterion has been proposed in [17] where the minimum eigenvalue of the Hessian matrix was maximized.

Selecting the rotation and loading speed is also a perspective. Coarsely, an appropriate sampling (to limit redundancies) in time is successive complementary projections (e.g., 0 • -90 • -180 • -270 • ...). For the measurement of smooth time evolutions, the rotation does not need to be fast and the sampling does not need to be important (in Chapter 2, 1.5 turns during the loading evolution was enough).

New Experimental procedures

The use of P-DVC allows the test to be performed very quickly. Different experimental procedures optimized up to now for standard full field measurements can be revisited. Various experimental procedure: The development of P-DVC opens up new way of designing an experiment. Three possibilities:

• Multi-view: performed with a single beam and at least two acquisitions per step, using a multi-source scan (as developed in medical scanners) or with a beam splitter. The different fields of view do not necessarily need to be at 90 • . The rule is to have enough sensitivity to be able to identify the motion (e.g., 10 degrees with a well textured sample and a well regularized kinematics could be enough, at the expense of an anisotropic resolution, see also TEM tomography where the tilting angle is limited [18]) if a rotation is superimposed. Multiple beams also offer additional benefits: two X-ray energy may provide an additional chemical contrast. X-ray and neutrons tomography can be combined as in the NeXT project. This last perspective is discussed further down.

• Single view acquisition: the experimental setup is simpler but a strong knowledge on the expected kinematics is needed to correctly choose the regularization. In this kind of experiment, the rotation speed, the loading rate and acquisition rate have to be correctly chosen. For the measurement of a smooth phenomenon, the rotation can be slow but for a fast experiment, with time localized phenomenon, the higher the rotation rate, the easier the identification.

• Dynamic tomography: still an ongoing project, dynamic tomography is a very promising topic when dealing with very fast phenomenon and medical imaging.

In all of those applications, every additional data can be introduced to help the kinematics to converge. An additional visible camera for example or a laser measurement are very accessible setups that could easily be used.

Obviously the P-DVC procedure can also be included in a standard 4D DVC experiment and would allow following the motion between acquisition steps (and follow a crack evolution for example).

Unconventional in-situ machine testing: One limitation when designing an in-situ testing machine is the frame of the machine that has not to move in front of the beam. Three solutions are developed: using a transparent frame (for tensile test [19,20] or 3 points bending [21]), enabling the sample to rotate in the machine (with at least one motor) or using a laminography setup for thin samples. It is to be noted that almost all mechanical in-situ testing machine are uniaxial testing machine. However, with the developed P-DVC procedures, the sample has not to be imaged and rotated during the experiment. This opens up new experimental in-situ procedures with standard machines. The sample could be imaged before the experiment, without the machine. Then, positioned in-situ , inside of the testing machine, the sample could only rotate in the possible range such that the acquisition is not disturbed by the frame shading. For example, one could even imagine using a bi-axial testing machine such as the one used in laminography experiments and in [22,23] as soon as few windows allow a projection acquisition at different angles. A test performed with a small hexapod testing machine that enables very complex multi-axial controls (see figure 10.4) could be a very interesting experiment. Real time tests: such as hybrid test that consists in driving a test, on the fly, from a numerical simulation with specific quantity of interests measured or computed. A real dialogue between experiment, measurements and numerical simulation is performed. Among the different categories,

• sub-structured hybrid tests [24,25] where a large structure is simulated and a small element of this structure that is critical or complex is tested experimentally. The computed boundary conditions are injected physically in the experiment through more or less complex actuators. The loop is closed when the measured quantities that are injected in the simulation.

• clone hybrid tests, where the simulated sample is the same as the measured one [26,27]. This method allows the boundary of a sample to be controlled from computed quantities of interest (e.g., following a predefined crack path) and generates very discriminant tests.

Both hybrid techniques are not suited for full 4D measurements because the acquisition is too long. However, P-DVC could be an interesting way to control a standard 4D experiment. The acquisition of the reference volume and sensitivity computation could be performed in an offline phase. Then, at the beginning of the experiment the linearized problem could be solved very quickly. During the acquisition of the full 3D scans, the volume and sensitivities could be updated.

Developing hybrid tests may require the used of complex multi-axial testing machine when modifying the boundary conditions in accordance with the numerical model (the use of multi-axial force sensors may be useful too). This perspective could be coupled with the unconventional in-situ machine testing one. Being able to perform such real time controlled tests is very interesting for parameters identification and to design very discriminating tests. Instead of acquiring volumes at predefined [σ,ε] positions (and discovering the results maybe too late), the acquisition states could be decided in accordance with the measured quantity evolution (e.g., damage, plasticity, crack position, etc).

Moreover, depending on the results, with a conic beam, a displacement of the sample and testing machine carried out with the turntable would allow zooming in interesting regions such as the crack tip area or on the contrary move back to have a more global determination of the boundary conditions. Then the projection operator should be updated with respect to the time evolution and could be given (at first order) by the tomograph sensors (as performed for stereo-vision in [28]).

Perspectives and future projects

NeXT-Grenoble: An extremely promising imaging instrument is the NeXT-Grenoble project2 developed at the 3SR laboratory and Institut Laue-Langevin (ILL) where a sample can be simultaneously imaged with neutrons and X-rays. Designed on the ILL neutron tomography beamline, an additional lab X-ray source is positioned perpendicularly to the neutron beam. Because of a large space in the casemate, sample can be imaged in-situ (under loadings and/or fluid injection for example, as developed in the PhD of Maddi Etxegarai at the 3SR -Grenoble). A multi-modality registration framework has been developed with X-ray and Neutron tomography data [29] allowing to exploit the complementarity between all modalities and find a 'common language' between X-Ray and Neutron images. The comparison of absorptions is presented figure 10.6 and shows the complementarity of the two waves. With this instantaneous double acquisition coupled with P-DVC, high speed phenomenon could be imaged without having to rotate the sample and exploiting all modalities advantages. A dried pine cone opening in a wet environment would be a funny (as funny as useless) and challenging kinematics-fluid measurement application of moving samples with gray level changes. In a similar spirit, medical dual-source or dual-energy CT [30,31] have a simultaneous acquisition with different energies (for example, 80 and 140 kVp), with orthogonal of common X-ray beam directions. This technology could also highly reduce the time resolution of the measured phenomenon. Images of the NeXT project and a dual-energy CT are shown figure 10.7.

Elastography in medical imaging: Biomechanical properties of living tissues and organs are of great importance as they contribute to or are responsible for tissue health and disease. The fine characterization of its properties (always elastic, often the Young modulus, sometimes the Poisson ratio) allows quantitative diagnosis (i.e., Oncology: mechanical properties of tumor [34]; Opthalmology: eye pathologies [35,36]; Applied to soft tissues: kidneys, breasts, liver, prostate, cells, fascia tissues, etc).

Over the past 25 years, a range of elastography techniques that replace palpation have been developed to image the mechanical properties of tissues [37,38]. Identification of the mechanical properties can be performed [39,40,41]. Two main categories [37,38]: static and dynamic approaches:

• Static / quasi static: Compression optical coherence elastography OCE [40,42,41]. . Image extracted from [32]. See also [33] for a material science example.

• Dynamic: Harmonic OCE [38], sono-elastography [43], Shear elastography [44,39].

However, most of elastography methods are limited in depth (3-4cm for ultrasound elastography, 1 up to 3mm for OCE depending on the material opacity, diffusion and scattering).

CT images have not, to the best of my knowledge, been used for elastography purposes and coupled with registration methods. A coupling between elastography and CT-P-DVC would allow, based on very low dose X-ray CT scans for both static and dynamic, a tissue characterization. This approach could be applied on standard CT scanners without additional setup (possibly with a specimen enriched with standard medical tracers in order to increase contrast and hence the sensitivity).

Dynamic tomography: Dynamic tomography has a huge potential in material and medical imaging. It gives access to the correct reconstruction of moving volumes. Applied in electron microscopy this would allow uncontrolled moving samples to be imaged. This mechanical tests, the procedure could be performed completely continuously. In medical imaging, the acquisition of fast moving phenomenon could be performed at low dose. This low dose approach could also be used for in-vivo motion measurements (high speed wing beats in a bee at approx. 120 Hz [45], modal in-vivo characterization of abdominal aortic aneurysm sealing at approx 1-20 Hz [46]) as well as ear-drum frequency characterization (much higher frequencies), etc.

Final words

The developed projection-based measurement allows, released from the volume reconstructions, the extremely fast 4D displacement field measurement of an experiment. The kinematics is read directly on the projections instead of volumes enabling a precious gain in acquisition time. Moreover, when exploiting the consistency of the projection/reconstruction operator, the reconstruction can be corrected from the motion-artifacts and the method can be extended to moving samples. Finally, when a gray level correction is included in the model of the procedure, the material and phase changes can be identified.

The method applied on plaster, cast iron, sandstone and fluid, medical images for various experiments: elastic brittle, plastic, capillary and pressure driven flow, vibration, CT calibration, CT motion compensation, has huge potentials.

This projection-based approach opens up new avenues when performing a test, faster, with reduced doses. The experiment are now focused on the key data (measurement, model parameters, etc), allowing a large part of the redundancies to be cut down. Each acquired data became thus more relevant for the experimental purpose. Coupled with other orthogonal CT developments such as synchrotron setups (whose acquisition could be performed at 270,000 frame rate!), multi-sources scanners, etc., or image regularizations methods, the procedure performance could be very much increased.

The 4D imaging of first in all times phenomena can be dreamt of.

Figure 1 . 2 :

 12 Figure 1.2: Standard data flow when reading tomography projections, composed of three successive inverse problems.

Figure 1 . 3 :

 13 Figure 1.3: (a) Flat-field projection showing a vertical gradient and white patterns (b) Radiograph with the same background (c) slice of the raw sinogram at z = 256 pixels without flat-field correction with bright vertical lines. Image extracted from[24]

Figure 1 . 4 :

 14 Figure 1.4: Reconstruction with different methods of a rat bone with (a) the ground truth, (b) the FBP with 20 projections, (c) SART, (d) TV regularization and (e) DART from[15]

Figure 1 . 5 :

 15 Figure 1.5: DVC application with on the top line the reference and deformed plasterboard volumes and on the bottom line the residual field and the measured displacement field. Figure extracted from[START_REF] Buljac | Digital volume correlation: Review of progress and challenges[END_REF]

Figure 1 . 6 :

 16 Figure 1.6: Comparison between the measured (first line) and the identified displacement fields (second line). The residual fields are presented in the third line. This figure is extracted from[START_REF] Pottier | Contribution of heterogeneous strain field measurements and boundary conditions modelling in inverse identification of material parameters[END_REF]

Figure 1 . 7 :

 17 Figure 1.7: Projection-based measurements: short-circuiting the data flow, from projections to measurement.

Figure 1 . 8 :

 18 Figure 1.8: Evolution of the vertical displacement field, measured from (a) DVC procedure, (b) P-DVC with all the 600 projections, (c) 48 projections and (d) only two orthogonal projections, extracted from[START_REF] Leclerc | Projection savings in ct-based digital volume correlation[END_REF]

Figure 1 . 9 :

 19 Figure 1.9: First line: two projections of the sample, at 0 and 90• . Second line: projected residual field at the end of the procedure. The top part of those residuals was masked and has not been corrected in the procedure and shows the initial error. This figure is extracted from[START_REF] Leclerc | Projection savings in ct-based digital volume correlation[END_REF]

Figure 1 .

 1 Figure 1.10: Projection-based identification of model parameters

  Figure 1.11: Projection-based measurement and identification procedures using a reference volume

  )

Figure 2 . 1 :

 21 Figure 2.1: (a) Drawing and 3D reconstruction of the sample and (b) the in situ testing machine [35] inside the tomograph with 1○, the testing machine with a carbon tube being supported by the rotating stage, 2 ○, the X-ray source and 3 ○, the detector. Note that the distance between the source and detector is not the one used in this article

Figure 2 . 2 :

 22 Figure 2.2: (a,b) Projection of the sample loaded with respectively Pb and Cu powder and (c) Different sizes of Fe and Pb powder (top and bottom)

Figure 2 . 3 :

 23 Figure 2.3: (a) Projection of crack propagation sensitivity and (b) norm of this sensitivity fields for different angles;

Figure 2 . 4 :

 24 Figure 2.4: Sinogram of the reference state at (a) s(r, θ = 0 • , t = 1) and (b) s(r, θ = 90 • , t = 1). z is the rotation axis, centered with respect to the image. 1 pixel → 30µm

Figure 2 . 5 :

 25 Figure 2.5: Sinogram of deformed specimen at ()a s(r, θ = 0 • , t = 8) and (b) s(r, θ = 90 • , t = 8). The vertical crack may be visible in (a) from the notches. The white dashed rectangle is the bounding box boundaries and the full white lines are the central mask boundaries.

Figure 2 . 6 :

 26 Figure 2.6: Artifact residual (a) ρ(r, 0 • , 1) and (b) ρ(r, 90 • , 1). The four quadrants of the detector are visible.

Figure 2 . 7 :

 27 Figure 2.7: Residual fields (a) ∆ρ(r, 0 • , t = 8) and (b) ∆ρ(r, 90 • , t = 8). The negative and positive patterns are the signature of a displacement

Figure 2 . 9 :

 29 Figure 2.9: Convergence of the variation of the mean SNR (left) and of the norm of the displacement variation (right)

Figure 2 .

 2 Figure 2.10: 2D residual field at step 8 after the displacement field correction for (a) θ = 0 • and (b) 90 • . The black disk in the center of the right image is a reconstruction artifact.

Figure 2 . 11 :

 211 Figure 2.11: Displacement field at convergence in voxels between the reference and step 8 in (a) X, (b) Y and (c) Z direction. The voxel size is 30 µm.

Figure 2 . 12 :

 212 Figure 2.12: 2D residual field at step 7, before the crack initiation, for (a-c) θ = 0 • and (b-d) 90 • . The first line is before the displacement field correction and the second after. The vertical lines that are apparent in (c-d) are a detector artefact (also noticeable at step 8).

Figure 2 . 13 :

 213 Figure 2.13: Mean SNR before and after the displacement field correction for all loading step.

Figure 2 . 14 :

 214 Figure 2.14: Comparison between the measured and computed force considering the identified Young modulus for the 13 steps. From step 11 and on, (and possibly earlier) the occurrence of damage is believed to make the comparison irrelevant

Figure 2 . 15 :

 215 Figure 2.15: Norm of the displacement field with the kinematic minimization only, u 1 , and with the force and kinematic minimization (including the identified Young modulus), u 2 , expressed in voxels

Algorithm 2

 2 4D-P-DVC fixed-point procedure while High residual norm do Initialize {a l } and {b l } Correction F (x) ← F (xu l-1 ) Compute updated projected sensitivities S i (r, t), Equation (3.10) while Ns j=1 ∆b l j Φ j (x) < p and ∆{α l } < α do Compute spatial mode {b l }, Equation (3.12) Compute temporal amplitude {a l }, Equation (3.15) end while Update displacement field u l , Equation (3.7) l = l + 1 end while

Figure 3 . 1 :

 31 Figure 3.1: (a) Radiograph of the sample (the rotation axis is at the center of the image). Gray levels have been normalized so that the maximum log-attenuation is set to 1, and that of air to 0. (b) Reconstructed volume, and (c) the proposed space regularization with rigid body motions on few sections of the beam

Figure 3 . 2 :

 32 Figure 3.2: In-situ tensile test with (a) the dog-bone sample used in the present procedure (measured size of the rectangular cross section of 1.31×0.91 mm 2 ). (b) Testing setup with 1 ○ the testing machine with a carbon fiber composite loading tube, 2○ the X-ray source, and 3○ the X-ray detector

Figure 3 .

 3 Figure 3.3 shows 3 selected projections at different times, angles and load.The particular choice of these angles will be discussed after the analysis of the results. The 60-pixel wide right and left edges of the radiographs and 100-pixel long top and bottom

Figure 3 . 3 :

 33 Figure 3.3: Projection at scale 1 of the sample at different angles and load during the tensile experiment. Projections are normalized between 0 (air) and 1 (maximum attenuation over a full rotation). (a) Time step 90, θ = 80 • , T = 630 N. (b) Time step 110, θ = -64.8 • , T = 715 N. (c) Time step 123, θ = -158.4 • , T = 736 N.

Figure 3 . 4 :

 34 Figure 3.4: One hundred twenty seven force measurements T (t) of the tensile test starting from 250 N. A first rotation is performed at constant load (corresponding to the first plateau until the 50th time step, or 100 s). Then the load is increased until failure. The three red dots are the projections shown in Figure 3.3 for time steps 90, 110 and 123.

Figure 3 . 5 :

 35 Figure 3.5: Initial residual field (keeping the projection normalization) composed of a large motion for (a) time step 90, θ = 80°, T = 630 N, (b) time step 110, θ = -64.8°, T = 715 N, (c) time step 123, θ = -158.4°, T = 736 N.

Figure 3 . 6 :

 36 Figure 3.6: Residual field after the rigid body motion correction. Tensile patterns not yet corrected are visible at (a) time step 90, θ = 80°, T = 630 N, (b) time step 110, θ = -64.8°, T = 715 N, (c) time step 123, θ = -158.4°, T = 736 N.

Figure 3 . 7 :

 37 Figure 3.7: Mean corrected translation measured with 8 rigid body motion modes for all loading steps with (1), (2) and (3) respectively corresponding to u x , u y and u z .

Figure 3 . 8 :

 38 Figure 3.8: Residual field at the end of the procedure. (a) Time step 90, θ = 80 • , T = 630 N, (b) time step 110, θ = -64.8 • , T = 715 N, (c) time step 123, θ = -158.4 • , T = 736 N.

Figure 3 . 9 :

 39 Figure 3.9: SNR of each projected residual field. The mean value, 9.9 dB for the initial SNR (a), 25.4 dB for the RBM corrected (b) and 26.6 dB for the corrected residual (c) shows an excellent description of the kinematics.

Figure 3 .

 3 Figure 3.10: (a) Measured vertical displacement field for each sections without rigid body motions, expressed in voxels, and (b) corresponding vertical strain. A large strained zone of approximately 8 -9% appears in the central part.

  Extracted from: A. Buljac, C. Jailin, A. Mendoza, J. Neggers, T. Taillandier-Thomas, A. Bouterf, B. Smaniotto, F. Hild, S. Roux, Digital Volume Correlation: Review of Progress and Challenges, Experimental Mechanics, 58(5):661-708, 2018

Figure 3 . 11 :

 311 Figure 3.11: Fig. 23 Projected residual field at angle θ = -150 • for step 110 (a) with no displacement field correction and (b) after the correction of the measured displacement. The black lines are the projections of the central cube element. Note that the gray level color bar (where gray levels are scaled to the full projection dynamic range) differs by a factor 2

Figure 3 .

 3 Figure 3.12: Fig. 24 (a) Measured displacement of the 8 nodes. During the first revolution (t < 50) a slight compression is visible, which is attributed to a relaxation phenomenon. In the increasing load part, the central region is stretched by approximately 11 voxels (or 110 µm). (b) SNR history for the initial residual fields (1), the residuals corrected by (x,y) rigid body translations (2), and the final residual fields with all corrections (3)

Chapter 4

 4 Identification from a single projection C. Jailin, A. Buljac, A. Bouterf, F. Hild, S. Roux, Fast 4D tensile test monitored via X-CT: Single projection based Digital Volume Correlation dedicated to slender samples, Journal of Strain Analysis for Engineering Design, 2018 Reproduced from

Figure 4 .

 4 1 shows 3 selected projections at different times (steps 90, 110 and 123). The first one (

Figure 4 . 1 :Figure 4 . 2 :

 4142 Figure 4.1: Projection at scale 1 of the sample at different angle and load during the tensile experiment. (a) Time step 90, θ = 80°, κ = 630 N; (b) time step 110, θ = -64.8°, κ = 715 N; (c) time step 123, θ = -158.4°, κ = 736 N

Figure 4 . 3 :

 43 Figure 4.3: Residual fields from volumes corrected by the rigid body motions. (a) Time step 90, θ = 80°, κ = 630 N; (b) time step 110, θ = -64.8°, κ = 715 N; (c) time step 123, θ = -158.4°, κ = 736 N

Figure 4 . 4 :

 44 Figure 4.4: SNR for all radiographs of the test with the different corrections. (a) Initial RBM-corrected residual fields, (b) residual fields with the kinematic measurement proposed in Ref. [36], (c) at the end of the proposed integrated approach

Figure 4 .

 4 Figure 4.5(a) shows the identified mechanical behavior of the studied sample using 3 different procedures (i.e., FEMU and Integrated P-DVC for the 8 parameter hardening law

  Figure 4.5(b). A displacement amplitude of 16 voxels (≈ 170 µm) is reached at the end of the experiment. High gradients appear in the central part, as expected from the geometry.

Figure 4 . 5 :

 45 Figure 4.5: (a) Stress-strain curve reconstructed from the converged parameters with different methods. (+, orange): initial guess, ( * , red): FEMU approach from displacement results of Ref. [1], ( , blue): integrated approach with 7 unknowns, (•, black): integrated approach with 8 unknowns. (b) Longitudinal displacement expressed in voxels (1 voxel↔ 10.7 µm), obtained from the identification of the mechanical law with the integrated approach with 8 parameters

Figure 4 . 6 :

 46 Figure 4.6: Residual fields from volumes corrected with the identified plastic law. (a) Time step 90, θ = 80 o , κ = 630 N; (b) time step 110, θ = -64.8 o , κ = 715 N; (c) time step 123, θ = -158.4 o , κ = 736 N

Figure 4 . 7 :

 47 Figure 4.7: Residual displacement field, expressed in voxels (1 voxel↔10.7 µm), obtained from the difference between the displacement field obtained from the integrated approach and that measured with the kinematic approach of Ref.[1]

  This work has been performed in Lund University under the supervision of Pr. S.Hall C. Jailin, M. Etxegarai, E. Tudisco, S. Hall, S. Roux, Fast tracking of fluid invasion using time-resolved neutron tomography, Transport in Porous Media, 1-19, 2018 Reproduced from

Figure 5 . 1 :

 51 Figure 5.1: Projections for (a) the reference state p 0 (r, t), (b) a state during the invasion process p(r, t) and (c) the difference that shows the projection of the 3D column of fluid.

Figure 5 . 2 :

 52 Figure 5.2: (a) Cylindrical sandstone sample with two, opposing flattened surfaces each containing a notch, which are offset to each other. In the pressure-driven flow experiment, the fluid entered the sample from the bottom. (b) Scheme of the experimental setup. The water tank can be moved up and down to keep a constant pressure during the experiment

Figure 5 . 3 :

 53 Figure 5.3: Comparison of the mean vertical intensity lines, from the model ζ ((a) black dashed curve), the blurred model with a Gaussian of two vertical pixels ((b) black bold plain curve) and experimental data ((c) red thin curve) extracted from a 3D reconstruction.

  The projections of the fluid column for the first and last radiographs are shown in Figure5.4. The difference between the first radiograph and the mirror symmetry of the last one indicate that the fluid front clearly moved during the 1 min acquisition. Such an evolution would blur the reconstructed 3D fluid field and, therefore, influence any subsequent 3D analysis.

Figure 5 . 4 :

 54 Figure 5.4: Projection of the 3D column of fluid at the first (a) and last (b) angle. The difference between the first and the (mirrored) last image in (c) shows the fluid advance during the scan.

Figure 5 . 5 :

 55 Figure 5.5: Results of the first mode identification with (a) the space field of the surface and (b) the time evolution of the front g 1 (t) in black and the free identification in red.

Figure 5 . 6 :Figure 5 . 7 :

 5657 Figure 5.6: Selected residual fields at (a) β = 0°or t = 0 s and (b) β = 180°or t = 60 s at convergence, with a flat and horizontal (constant height) surface.

Table 5 . 1 :Figure 5 . 8 :

 5158 Figure 5.8: Surface modes m = 2 and 3. (a-c) is the space surface F m (x, y) and (b-d) is the time evolution g m (t)

Figure 5 . 9 :

 59 Figure 5.9: Evolution of the front at different times of the scan acquisition, in voxels, reconstructed from the three modes for (a) t = 0 and (b) 60 seconds

( 7 -

 7 21-40) and 3 different angles β =[240°, 0°, 292°] are shown figure 5.11.

Figure 5 . 11 :

 511 Figure 5.11: Projections of the sample during the pressure-driven flow test (a) step 7, β = 240°, (b) step 21, β = 0°and (c) step 40, β = 292°T

Figure 5 . 12 :

 512 Figure 5.12: Residual field at convergence of the procedure, with the correction of 5 identified modes during the test (a) step 7, β = 240°, (b) step 21, β = 0°and (b) step 40, β = 292°

Figure 5 . 13 :

 513 Figure 5.13: Three first identified modes respectively (a-b-c) for the surface field and (d-e-f) for the time evolution

Figure 6 . 2 :

 62 Figure 6.2: Projection of the synthetic volume.

Figure 6 . 3 :

 63 Figure 6.3: (a-b) two noisy residual fields at the same angle and different acquisition times.

Figure 6 . 4 :

 64 Figure 6.4: Measured displacement field Φ(x) on the T4 mesh for the first test case.

Figure 6 . 5 :

 65 Figure 6.5: (a): comparison of the ground truth and measured amplitude for each projection α(t) of the three angles (1-2-3) respectively for the intervals [0-150-300-450] and (b): the ground truth in (1) and the error in (2).

Figure 6 . 6 :

 66 Figure 6.6: Two noisy residuals fields of the second test case acquired at the same angle and different time acquisition.

Figure 6 . 7 :

 67 Figure 6.7: Measured displacement field Φ 1 (x) for the first mode of the second test case.

Figure 6 . 8 :

 68 Figure 6.8: (a): comparison of the ground truth and measured amplitude for each projection α 1 (t) of the three angles (1-2-3) respectively for the intervals [0-150-300-450] and (b): the ground truth in (1-2-3) for the three angles and the error in (4)

Figure 6 .

 6 Figure 6.10: (a): comparison of the ground truth and measured amplitude for each projection α 1 (t) of the three angles (1-2-3) respectively for the intervals [0-150-300-450] and (b): the ground truth in (1-2-3) for the three angles and the error in (4)

Chapter 7

 7 Online tomograph self-calibration C. Jailin, A. Buljac, A. Bouterf, M. Poncelet, F. Hild, S. Roux, Self-calibration for lab-µCT using space-time regularized projection-based DVC and model reduction, Measurement Science and Technology, 29(2):024003, 2018 Reproduced from

Figure 7 . 1 :

 71 Figure 7.1: (a) Acquired projections for θ = 0°. (b) Corresponding reconstruction

Figure 7 . 2 :Figure 7 . 3 :

 7273 Figure 7.2: Slices of the reconstruction of step 1 for (left) y = 500 voxels and (right) z = 150 voxels. The top and bottom parts are slightly blurred due to motion or tilted rotation axis

Figure 7 . 4 :

 74 Figure 7.4: Initial and final residual fields after the correction of the three modes when (a) θ k = 90°and (b) θ k = 330°T

Figure 7 . 5 : 4 Figure 7 . 6 :

 75476 Figure 7.5: Initial and final angular SNRs for step 2 with the full measurement approach. The two arrows show the position of the extracted residuals of Figure 7.4

Figure 7 . 7 :Figure 7 . 8 :

 7778 Figure 7.7: Zoom on the reconstruction of step 2 at slice z = 100 voxels (a) for the initial state and (b) after correction of the rotation axis

Figure 7 . 9 :

 79 Figure 7.9: SNR(θ) for all treated angles after successive additions of modes. The two arrows show the position of the extracted residuals of Figure 7.10

Figure 7 . 10 :

 710 Figure 7.10: Final residual fields after the correction with three modes for (a) θ k = 90°and (b) θ k = 330°T

Figure 7 . 11 :Figure 7 . 12 :Figure 7 . 13 :

 711712713 Figure 7.11: First mode. (a) Temporal mode with (black circles) and without (blue points) time constraints. (b) Spatial mode in the x, y, z directions, expressed in voxels on the mesh (1 voxel ≡ 5.6 µm)

Figure 7 . 14 :

 714 Figure 7.14: Mean displacement of the top, center and bottom parts of the sample in the laboratory frame and its projection on the detector frame. The axis of the X-ray cone beam is X

Figure 7 . 15 :

 715 Figure 7.15: Zoom on the reconstruction of step 1 for slice z = 100 voxels (a) for the initial state, (b) after correction of the rotation axis

Algorithm 6

 6 Standard algebraic reconstruction procedure, ART (p) Initialization n ← 1 Initialization f (n) ← 0 Choose a permutation, π, over N θ indices while ρ(r, θ) > do for k ← 1 to

Figure 8

 8 Figure 8.1: a) Reference image and the 9-node mesh the node of which are subjected to a random displacement, assumed to be linear in time; b) deformed phantom at final time N θ .

Figure 8 . 2 :Figure 8 . 3 :Figure 8 . 4 :

 828384 Figure 8.2: a) Initial reconstruction with u(x, t) = 0; b) initial projected residual fields ρ(r, t). Note that the color amplitude that is saturated in this image has been selected to be the same with the corrected residuals shown further down (figure 8.3)

Figure 8 . 5 :

 85 Figure 8.5: a) Reference image (unknown) and applied nodal displacement field V 1 (x) in red and V 2 (x) in light blue; b) deformed checkerboard at time N t .

Figure 8 . 6 :Figure 8

 868 Figure 8.6: a) Initial reconstruction with u(x, t) = 0; b) initial projected residual fields ρ(r, t)

Figure 8 . 8 :Figure 8 . 9 :

 8889 Figure 8.8: Zoom in the reconstructed volume (a) with a standard non-corrected volume and (b) with the proposed procedure

Figure 9 . 1 :

 91 Figure 9.1: Scheme of the procedure that required only the commonly acquired sinogram.The framework is separated in two iterative procedures: the measurement of the motion from a reconstructed volume and the corrected reconstruction. At the beginning, the motion is initialized at 0.

Figure 9 . 2 :

 92 Figure 9.2: Multi-scale framework starting from the coarse scales to the finest. The volumes are blurred with a Gaussian kernel (and downsized) to enhance the sensitivity to the global kinematics. The model also starts with a coarse resolution and is progressively enriched with respect to the residual fields. In the following application, the two last meshes at scale 'medium' and 'fine' are the same.

Figure 9 . 3 :

 93 Figure 9.3: Synthetic test case generated from 10 real CT-scans. (a) comparison of the projection at 90 • of two 3D scans (the two red dots in the right figure) and (b) sequence of the simulated projections. Each one of the 8 angles is the projection of the same state so that the entire test is composed of 4.5 respiratory cycles that is realistic.

Figure 9 .

 9 3(a) shows the difference of two respiratory state at one projection (90 • ) highlighting the vertical motion of the diaphragm. This displacement is essentially vertical and localized at the diaphragm. The right plot (b) is the projection sequence showing the way the 10 CT scans are projected on the 360 projections.

Figure 9 . 4 :

 94 Figure 9.4: Measured displacement field for the time acquisition at angle 48 • . In the left image, the vector field is plotted on a coronal slice and in the left image, the vertical displacement field is shown. The displacement is expressed in voxels with 1 voxel in [x, y, z] ↔ [1.95,1.95,2.0] mm. In (c), the final displacement field is plotted with arrows on the coronal view.

Figure 9 . 6 :

 96 Figure 9.6: Projected residual fields (sinogram of 2D projections + angle) before (left) and after (right) correction by the measured displacement field. The oscillating positive and negative patterns are the signature of the periodic respiratory motion. After correction by the measured displacement field, a large part of the error has been corrected.

[ 1 ]

 1 W. A. Kalender. X-ray computed tomography. Physics in Medicine & Biology, 51(13):R29, 2006.

Figure 10 . 1 :

 101 Figure 10.1: Evolution of the spatial and temporal resolution for X-ray imaging, extracted from [1] and enriched with 5 crosses corresponding to the developments of this thesis. Blue crosses stand for actually performed experiments, yellow cross corresponds to a synthetic trial. The different cases are defined in Table 10.1.

Figure 10 .

 10 3 presents 3 new fast experimental measurement protocols described above.

Figure 10 . 3 :

 103 Figure 10.3: New experimental 4D procedure using P-DVC framework.

Figure 10 . 4 :

 104 Figure 10.4: Hexapod testing machine that could be positioned in-situ. The three red arrows indicate the possible acquisition windows where the beam is not masked by the 6 actuators.

Figure 10 . 5 :

 105 Figure 10.5: Scheme of a sub-structured hybrid test where a small part of the structure is simulated.

Figure 10 . 6 :

 106 Figure 10.6: Example of absorption with Neutron (top) and X-Ray (bottom) (images from PSI)

Figure 10 . 7 :

 107 Figure 10.7: (a) NeXT-Grenoble project combining X-ray and Neutron tomography and (b) Dual-energy medical X-ray CT (Siemens Healthcare technology). Image extracted from[32]. See also[33] for a material science example.

Table 1 .

 1

	1: Comparison Stereo-Correlation / P-DVC
	Stereo-Correlation	P-DVC methods
	Cameras	X-ray / detector
	Calibration matrices

  4 ,2H 2 O)). Two notches (3.80 mm long from the center and 200 µm thick) are introduced in the mold. Finally a 2.30 mm diameter hole is drilled in the center.

Table 2 .

 2 .1 summarizes the main characteristics of each loading step. 1: Loading steps information for the two parts of the experiment

	Number	1	2	3	4	5	6	7
	full 3D scan Yes No No No No No No
	load [N]	0	91.8 188 305 393 488 609
	visual crack No No No No No No No
	Number	8	9	10	11	12	13
	full 3D scan No No No No No Yes
	load [N]	628 638 651 669 691 625
	visual crack Yes Yes Yes Yes Yes Yes

Table 4 .

 4

	1: Eight initial control points of the free form hardening law (8 strain unknowns
	are calibrated)
	Parameters -γ 1	γ 2	γ 3	γ 4	γ 5	γ 6	γ 7	γ 8
	σ [MPa]	0 235 335 435 525 555 595 605 620
	ε [%]	0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6

Table 4 .

 4 2: Calibrated values of the 8 strain control points of the free form hardening law

	Parameters -	γ 1	γ 2	γ 3	γ 4	γ 5	γ 6	γ 7	γ 8
	σ [MPa]	0 235 335 435 525 555 595 605 620
	ε [%]	0 0.20 0.32 0.50 1.42 1.62 4.26 5.96 6.16

Table 7 .

 7

	1: Mean SNR as a function of the number of correction modes. A fast convergence
	with three modes is noted
	Mode l	0	1	2	3	4
	<SNR> l 24.5 25.8 26.7 26.8 26.8

  1 . The initial volumes where obtained from a time sampling based on the acquisition of a respiratory surrogate signal (Pneumo Chest pressure belt).

We propose in this paper to re-generate a moving sinogram by projecting the 10 volumes of 256 × 256 × 141 voxels successively at 360 different angles. Each displacement state

Table 10 .

 10

		1: Resolution enrichments
		Name	Ref. Space resolution Time resolution
	(1)	Multi-view P-DVC	[2]	15 µm	35 s
	(2)	Online calibration	[5]	2.8 µm	7.2 s
	(3)	Single view P-DVC	[3, 4]	10.7 µm	2.5 s
	(4) Single view neutron tracking	[6]	110 µm	200 ms
	(5)	Lab-CT modal P-DVC	-	>10 µm	up to 10 ms

A dynamic flatfield correction procedure has been proposed during this thesis[24] and enables, from the edges of the projections that are not shadowed by the sample, to estimate an optimum combination of the flatfields library.

https://www.creatis.insa-lyon.fr/rio

https://next-grenoble.fr/

Acknowledgement

This work has benefited from the support of the French "Agence Nationale de la Recherche" through the "Investissements d'avenir" Program under the reference ANR-10-EQPX-37 MATMECA and ANR-14-CE07-0034-02 COMINSIDE.

Acknowledgments

This work has benefited from the support of the French "Agence Nationale de la Recherche" through the "Investissements d'avenir" Program under the reference ANR-10-EQPX-37 MATMECA.

Bibliography

Algorithm 4 Modal P-DVC fixed-point procedure

Compute all sensitivities S i (r, t) ← Π θ(t) [f (x)ψ i (x)] while High residual norm do Initialize α k , such that α 2 l = 1 and possibly t l α(t) = 0 while ∆Φ k < u and ∆α k < α do Compute spatial mode u k (eq. (6.2)) Compute temporal amplitude α k (eq. (6.8)) Normalize to α 2 l = 1 and possibly

end while Update displacement field u k , or residual fields k ← k + 1 end while

Application of the method

The method has been tested with two virtual examples. A numerically generated specimen of size 384×384×386 voxels is used. The truncated cone shape is defined from a sample scanned for another application in the LMT lab-CT. A 3D view of the sample is presented figure 6.1. Its microstructure is a light matrix on which absorbing particles are included. The particles are described by 3D Gaussian profile. The volume was generated using an efficient matlab function provided by Jan Neggers. The numerical test case is generated with 3 projection angles at [0 -120 -270] • . Using 3 projections is indeed interesting because at least 2 will be sensitive. With a regular angular sampling and a prime number of projections (min 3) the displacement can always be correctly captured.

The set-up geometry is a cone beam (with parameters corresponding to a real experiment from which the shape is extracted) and projections of the volume with the ASTRA toolbox [13]. The volume is deformed considering mode shapes, with random amplitudes The motion of a sample while being scanned in a tomograph prevents its proper volume reconstruction. In the present study, a procedure is proposed that aims at estimating both the kinematics of the sample and its standard 3D imaging from a standard acquisition protocole (no more projection than for a rigid specimen). The proposed procedure is a staggered two-step algorithm where the volume is first reconstructed using a "Dynamic Reconstruction" technique, a variant of Algebraic Reconstruction Technique (ART) compensating for a "frozen" determination of the motion, followed by a Projection-based Digital Volume Correlation (P-DVC) algorithm that estimates the space/time displacement field, with a "frozen" microstructure and shape of the sample. Additionally, this procedure is combined with a multi-scale approach that is essential for a proper separation between motion and microstructure. A proof-of-concept of the validity and performance of this approach is proposed based on two virtual examples. The studied cases involve a small number of projections, large strains, up to 25%, and noise.

Discussions 10.2.1 Technical aspects

Unified reconstruction and measurement framework: The correct procedure when dealing with a sequence of successive inverse problems is to transport the covariance matrix, from the raw data to the final parameters. However, it is not easy to evaluate the uncertainty for reconstructed volumes (corrupted with reconstruction artifacts e.g., rings). Being able to work directly with the projections provides a more easily accessible description of the noise. The dynamic procedure is indeed based on the exact same functional split in two for the displacement and microstructure identification.

Computation time: Computation time can be an issue in P-DVC procedures. The heaviest operations are: (i) computation of the gray level sensitivity fields (i.e., weighted gradient), Π θ [∇f (x)ψ i (x)], (ii) update of the volume (i.e., 1 full 3D registration per loading step: 360 times per measurement iteration in the dynamic tomography of the lungs) and (iii) projections (and reconstructions) performed with the ASTRA toolbox. Different approaches can be implemented such as working with GPU (already performed with ASTRA), update of the sensitivities and volumes at specific times, parallel computation, smart initializations based on RBM measurements (computed in Fourier space), etc.The use of PGD algorithm may be very useful when dealing with a large number of degrees of freedom that can easily be computed (the fluid tracking test case for example).

Measurement accuracy

The measurement accuracy depends on various factors, a large part of which can be tuned and optimized prior or during the experiment. Figure 10.2 shows 4 possible levers to enhance the measurement accuracy: the choice of the kinematic regularization, the sample shape and texture optimization, the imaging procedure (angles, frame rate, rotation speed, frame average,etc.) and the mechanical test procedure (type of testing machine, loads). Some other parameters are not considered in this figure such as the noise and artifacts because they are associated to the acquisition itself and they can be reduced when dealing with the frame averaging for example. A measurement method, called Projection-based DVC (P-DVC), shortens the previous sequence and identifies the kinematics directly from the projections. The number of radiographs needed for tracking the time evolution of the test is thereby reduced from 500 to 1000 down to 2. This thesis extends this projection-based approach to further reduce the required data, letting faster phenomena be captured and pushing the limits of time resolution. Two main axes were developed: -On the one hand, the use of different spatial and temporal regularizations of the 4D fields (space / time) generalizes the P-DVC approach (with a known reference volume) to the exploitation of a single radiograph per loading step. Thus, the test can be carried out with no interruptions, in a few minutes instead of several days. -On the other hand, the measured motion can be used to correct the reconstructed volume itself. This observation leads to the proposition of a novel procedure for the joint determination of the volume and its kinematics (without prior knowledge) opening up new perspectives for material and medical imaging where sometimes motion cannot be interrupted. The development of these two axes opens up new ways of performing tests, faster and driven to the identification of key quantities of interest. These methods are compatible with the recent "hardware" developments of fast tomography, both at synchrotron beamlines or laboratory, and save several orders of magnitude in acquisition time and radiation dose.