Measurement of 1–10 Hz 3D vibration modes with a CT-scanner

Abstract

A new in situ vibration mode measurement method within a tomograph is proposed based on Projection-based Digital Volume Correlation techniques. Several projection angles are selected and a large number of radiographs of the vibrating sample are acquired at random instants with a small exposure time in order to ‘freeze out’ the displacement and avoid motion blurring. Based on an initial reconstruction acquired in a static configuration, the displacement field measurement is performed using a Proper Generalized Decomposition technique. All projections are analyzed as being due to a few vibration modes deforming the known reference volume. The different projection directions are related to each other assuming that the modal amplitude probability distribution functions are statistically similar. A synthetic test case, mock-up of a liver, is used to illustrate and validate the approach. In this case, 5 projection angles were chosen, 300 radiographs per angle, and the first three vibration modes could be recovered with a good accuracy.

Introduction

In the last decades, full-field measurement methods have become a key element for the validation and identification of (bio-)mechanical models [1, 2]. Model identification is an inverse problem that consists in estimating the constitutive parameters of a mechanical model. The maturity of these approaches is such that one may argue that full-field based identification can be properly considered as a genuine measurement. Targeting identification or validation of a model (i.e., with a high signal to noise ratio), an experiment can be designed and optimized to be discriminant, i.e., highly sensitive to parameter changes [3]. Various levers can be activated to enhance sensitivity and hence mechanical model identification:

• selecting adequate measurement device, with surface measurements using e.g., visible camera(s), high speed devices, thermal measurements, or volume analyzes: X-ray or neutron Computed Tomography (CT) [4], etc.

• choosing the experiment procedure and loadings: e.g., multi-axial, with thermo-mechanical loadings [5], cyclic loads [3], controlled crack paths [6, 7], vibrations [8], etc.

• optimizing the sample shape (e.g., with holes or notches to enhance a plastic behavior as in [9] or more sophisticated methods based on topology optimization [10]).

Being able to measure the vibration displacement field (thus the modal basis) of a sample, in 3D, is both relevant for model discrimination [11] but also very challenging. And, even if the time evolution of the kinematics is not known, the determination of the modal basis of a sample (even limited to its first few modes) is extremely sensitive to the density and stiffness of the studied samples, and can therefore be used to (un)validate a model in a very discriminating way. Although stereo-DIC methods coupled with high speed setups appear to be very efficient and relevant measurement procedures [8, 12], they are limited to surface measurements [13]. Coupling vibration measurement and CT imaging is very appealing. However, vibrations with multiple modes (say from 1 to 10 Hz) are much too fast (or lack perfect periodicity) to be captured by standard tomography acquisition approaches.

Concerning general mechanical tests, coupled with 3D imaging devices such as X-ray CT scanners, in situ measurements of the displacement field using Digital Volume Correlation (DVC) [14] offer remarkable opportunities for identification. However, temporal resolution is a critical issue. To run standard 3D full field measurement techniques, with its classical 4D (space-time) implementation [15], several volumes are required, and the acquisition of each volume costs tens of minutes or hours in a lab-CT (although a much faster rate, up to 20 Hz, can be achieved at synchrotron facilities [16, 17], leading to 4D scan movies). A classical analysis would then require quasistatic acquisitions to avoid blur in the reconstructions, what would imply a period larger than several tens of hours. 10 Hz vibration would thus be out of reach by no less than a factor of $10^5$ in frequency, or 5 orders of magnitude.

Different methods have been developed to be able to acquire tomographic 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 [4].

• For periodic motion, the acquisition can be triggered by other sensors giving the phase of the motion [18]

• An extension of DVC where the measured quantities are identified on projections is the Projection-based DVC (P-DVC) [19,20,21]. Those developments allow a 4D (space/time) kinematics to be followed from a single projection per state. The developed P-DVC technique could also be applied to follow deformations radiograph after radiograph, while the sample is rotating [22, 23].

• Some dynamic tomography techniques [24], aim at imaging a moving sample by capturing the kinematics and reconstructing the sample at the same time.

Although those methods have allowed extending measurements to fast phenomenon, they are not able to deal with vibrations.

In [25], 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 out the displacement without motion blurring. Without time information, frequencies cannot be measured. However, images are related to each others by a kinematic field that is mainly composed of the first low frequency modes. From DIC analysis of those images, a weighted principal component analysis (PCA) [26, 27] could extract the first vibration modes. The latter were favorably compared with numerical results. Such an approach overcomes the obstacle of not being able to resolve finely the time evolution of the vibrations (because acquisition instants of time are random) and only requires fast image capture to avoid blurring. Here and in the following, the acronym PCA will be used, although a number of different acronyms exist with essentially the same meaning, such as Proper Order Decomposition (POD), Singular Value Decomposition (SVD), Karhunen-Loeve (KL), etc.

The previous analysis [25] motivates the development of a similar technique to extract 3D vibration modes with a CT-scanner, by-passing the severe limitation of the slow acquisition of 3D images. P-DVC is a way to reduce the time scale of tomographic acquisition to that of radiographs (leading to several orders of magnitude gain in time scale). However, a new challenge arises from the fact that each radiograph brings only a partial information on the kinematic field, and the random acquisition times forbid any regularization that would exploit a smooth temporal evolution. The present study addresses this challenge based on a statistical assumption of a proper sampling of the first vibration modes. It will be shown that exploiting such a property will allow to relate the different projection angles, and hence reassemble from the various observed projections the spatial vibration modes.

A standard static acquisition is first performed in a reference configuration. Five projection angles for the multi-view P-DVC are selected. A large series of radiographs (300 frames) is acquired at each angle, at random instants, using a short exposure time. To mention a realistic order of magnitude, 80 Hz is a frame rate (per single projection) that can be achieved in a lab-CT (whereas 1 to 10 kHz [16, 28] can be reached at synchrotron beamlines). The identification of the kinematics is then performed with a Proper Generalized Decomposition [29, 30] procedure computing sequentially the relevant modes using statistical assumption(s) on the displacement field. A synthetic application for the measurement of free vibration modes in a human liver shows the accuracy of the proposed method with a CT-scanned liver that is virtually deformed with 3 numerical modes.

In “Modal measurement method section”, after having introduced the data acquisition protocol, and recalled the principles of the Projection-based Digital Volume Correlation, the methodology of the PGD modal analysis is presented, together with the way to extract the spatial vibration modes from the PGD modes. “Application of the method” section presents an application to a numerical mock-up of the free vibration of a liver. Synthetic projections are computed based on a numerical model, trying to reproduce the actual complexity of a real application (namely, with added noise and no contrast within the organ, so that the projections contain only boundary information). Yet, those synthetic data will be shown to lead to a good reproduction of the first three vibration modes. “Conclusion and discussion” section concludes this paper with a summary of the major results and suggests some perspectives.

Modal measurement method

Problem definition and assumptions

The studied sample is first CT-scanned in a static configuration. A set of projections $p^{\text {ref}}$ is acquired when the source and detector are both rotated by an angle $\theta$ with respect to the sample to refer to a medical imaging set-up as schematized in Fig. 1 (or equivalently the sample is rotated in a static set-up including source and detector, as in lab CT scanners). ${{\varvec{r}}}$ denotes the 2D detector pixel position, and $p^{\text {ref}}({{\varvec{r}}},\theta )$ is the cologarithm of the ratio of X-ray intensity at pixel ${{\varvec{r}}}$ to the flat-field at the same position. From $p^{\text {ref}}$, one may compute a 3D image where each voxel ${{\varvec{x}}}$ is characterized by its coefficient of X-ray absorption, $f({{\varvec{x}}})$. This image is computed (or reconstructed) from the inversion of the linear relation

$$\begin{aligned} \Pi _{\theta } [f({{\varvec{x}}})]=p^{\text {ref}}({{\varvec{r}}},\theta ) \end{aligned}$$

(1)

where $\Pi _\theta$ is the projection operator [31] that sums the coefficient of absorption $f({{\varvec{x}}})$ along the beam that hits the detector at position ${{\varvec{r}}}$. In practice, $f({{\varvec{x}}})$ is computed from the minimization of

$$\begin{aligned} \Gamma _{\text {ref}}(f) = \dfrac{1}{\beta _p^2} \sum _{{{\varvec{r}}},\theta } \left\| \Pi _{\theta } [f({{\varvec{x}}})]-p^{\text {ref}}({{\varvec{r}}},\theta )\right\| ^2 \end{aligned}$$

(2)

where $\beta _p^2$ is the detector noise variance.

Fig. 1

Scheme of the acquisition procedure. The point-like source (shown as a $\bullet$) and the detector may rotate around the sample by an angle $\theta$

The sample is then assumed to be excited to a steady vibrating motion. After having scanned the reference configuration, the experiment consists in recording several sets of projections, $p^{\text {vib}}({{\varvec{r}}},t)$. These radiographs of the vibrating sample are acquired at random times using the shortest possible exposure time to avoid motion blur. Exposure times of order 0.0125 s are available in some micro-lab-CT (as the one at the LMT Paris-Saclay). Few, $N_l$, projection angles are chosen $\theta _l$, and for each of them, $N_t$ projections are recorded. The projection angle attached to a label t is denoted $\theta (t)$. These projections are indexed by a time label t, however their acquisitions are random instants does not allow t to be fully interpreted as a time, as it would if the experiment were performed with a very fine temporal sampling, i.e., much smaller than the inverse of the highest observed frequency.

Otherwise, any additional sensor information (i.e., laser interferometry, acoustic measurements) could potentially be used in a temporal regularization if available. In the sequel, no such supplementary information is assumed to be available.

Its 3D motion ${{\varvec{u}}}({{\varvec{x}}},t)$ is described by a displacement field that can be decomposed on its vibration modal basis as

$$\begin{aligned} {{\varvec{u}}}({{\varvec{x}}},t) = \sum _{\text {mode~}i} \alpha ^{\text {vib}}_i(t) {{\varvec{\Phi }}}^{\text {vib}}_i({{\varvec{x}}}) \end{aligned}$$

(3)

with ${{\varvec{\Phi }}}^{\text {vib}}_i({{\varvec{x}}})$, the ith spatial vibration mode, and $\alpha ^{\text {vib}}_i(t)$ its temporal amplitude. Vibration displacements are assumed to be of small amplitude (small to be compared with the characteristic length of the texture) and centered around the reference configuration (i.e., null mean displacement).

P-DVC

Projection-based DVC (P-DVC) aims at evaluating a displacement field, ${{\varvec{u}}}({{\varvec{x}}},t)$, from projections based on the consistency of the measured projection and the one computed from the warping of the reference volume by ${{\varvec{u}}}$. In the present case, considering vibrations, the functional to minimize is similar to the one used for reconstruction, at the exception of the warping. It reads

$$\begin{aligned} \Gamma _{\text {vib}}({{\varvec{u}}}) = \dfrac{1}{\beta _p^2} \sum _{{{\varvec{r}}},t} \left\| \Pi _{\theta (t)} [f({{\varvec{x}}} + {{\varvec{u}}}({{\varvec{x}}},t))]-p^{\text {vib}}({{\varvec{r}}},t)\right\| ^2 \end{aligned}$$

(4)

In the above expression, the functional $\Gamma _{\text {vib}}$ is normalized using the variance of the detector noise, $\beta _p^2$.

P-DVC is a non-linear problem because the unknown displacement field comes in as an argument of the sample texture f(.). However using multiscale approaches, the problem is usually addressed from successive linearizations, treating the displacement correction as a small perturbation, for which a Taylor expansion to first order is sufficient. In the present case, it is assumed that the displacement magnitude is small enough so that a multiscale approach is not needed and a mere linearization about the reference state is sufficient and does not necessitate successive further updates. In practice, this means that at each point of the sample the displacement amplitude is smaller than the microstructure correlation length. Thus, the above functional can be linearized to

$$\begin{aligned} \Gamma _{\text {vib-lin}}({{\varvec{u}}})= & {} \dfrac{1}{\beta _p^2} \sum _{{{\varvec{r}}},t} \left\| \Pi _{\theta (t)} [f({{\varvec{x}}})]+ \Pi _{\theta (t)} [{\varvec{\nabla }}f({{\varvec{x}}})\cdot {{\varvec{u}}}({{\varvec{x}}},t)]-p^{\text {vib}}({{\varvec{r}}},t)\right\| ^2\nonumber \\= & {} \dfrac{1}{\beta _p^2} \sum _{{{\varvec{r}}},t} \left\| \Pi _{\theta (t)} [{\varvec{\nabla }}f({{\varvec{x}}})\cdot {{\varvec{u}}}({{\varvec{x}}},t)]-(p^{\text {vib}}({{\varvec{r}}},t)-p^{\text {ref}}({{\varvec{r}}},\theta (t)))\right\| ^2 \end{aligned}$$

(5)

It is convenient to introduce $\rho ^0({{\varvec{r}}},t)=p^{\text {vib}}({{\varvec{r}}},t)-p^{\text {ref}}({{\varvec{r}}},\theta (t))$, the projection difference between the vibrating state at time t, and the reference state with the same projection angle. This mere difference has the additional merit of cancelling some spurious artefacts due to the X-ray detectors, and to avoid possible reconstruction inaccuracies (for instance due to approximate acquisition geometry parameters).

The above formulation is a very ill-posed problem and it calls for regularizations to be solvable. This holds when dealing with a motion that is smooth in time, or when prior knowledge can be incorporated from say a parametrization based on additional measurement such as a load. However, in the present case, the randomness of the time acquisition prevents such a strategy.

Instead, a statistical approach is proposed, based on a modal decomposition.

PGD modal analysis

The basic idea is to introduce modes as separated time functions and space fields, in the same spirit as Eq. 3.

$$\begin{aligned} {{\varvec{u}}}({{\varvec{x}}},t) = \sum _{\text {mode~}i} \alpha ^{\text {PGD}}_i(t) {{\varvec{\Phi }}}^{\text {PGD}}_i({{\varvec{x}}}) \end{aligned}$$

(6)

However, it is important to note that nothing guaranties at this stage that the modes which are identified from the measurements ${{\varvec{\Phi }}}^{\text {PGD}}_i$ coincide with the vibration modes ${{\varvec{\Phi }}}^{\text {vib}}_i$.

To render the problem solvable, one needs another assumption, that is similar to the one used to associate PCA modes and vibration modes [26, 27]. The additional assumption is that the displacement fields have a hierarchical structure. Namely, it is assumed that the first mode dominates over all others, so that—to dominant order—the above sum can be truncated to the first mode only. However, when the first mode is understood, and both $\alpha ^{\text {PGD}}_1(t)$ and ${{\varvec{\Phi }}}^{\text {PGD}}_1({{\varvec{x}}})$ are identified, they can be subtracted off from the data, $\rho$, and more modes appear, and again, it is assumed that the second mode dominates over all the other ones, and similarly order after order. This assumption is essentially based on a clear separation of eigen-frequencies (large ratio of consecutive resonant frequencies). Such an assumption leads naturally to a greedy approach where only one mode is identified at a time. Modes are extracted iteratively from the remaining projection data, $\rho$, from which the contribution of the previous modes has been erased. Thus, the elementary procedure is to determine the displacement field as if it were composed of a single mode. In this section, the superscript PGD is omitted for the sake of conciseness.

Displacements are discretized in space using a finite element mesh, fitted to the sample shape, as usually done in global DIC or DVC procedures [32, 33]. To specify our notations, mesh shape functions ${\varvec{\psi }}_j({{\varvec{x}}})$ are introduced, so that mode i, ${{\varvec{\Phi }}}_i$ reads

$$\begin{aligned} {{\varvec{\Phi }}}_i({{\varvec{x}}}) = \sum _j \phi ^j_i {\varvec{\psi }}^j({{\varvec{x}}}) \end{aligned}$$

(7)

With this FE discretization, mode i reduces to a vector ${\varvec{\phi }}_i$, that gathers all components $\phi _i^j$.

To determine a single mode, the space-time separation of ${{\varvec{u}}}(x,t)=\alpha (t){{\varvec{\Phi }}}(x)$ is inserted in the functional $\Gamma _{\text {vib-lin}}$ (since only one mode is considered at a time, the index of the mode is dropped in the following)

$$\begin{aligned} \Gamma _{\text {vib-lin}}({{\varvec{\Phi }}}(x),\alpha (t))= & {} \dfrac{1}{\beta _p^2} \sum _{{{\varvec{r}}},t} \left\| \alpha (t)\Pi _{\theta (t)} [{\varvec{\nabla }}f({{\varvec{x}}})\cdot {{\varvec{\Phi }}}({{\varvec{x}}})]-\rho ^{0}({{\varvec{r}}},t)\right\| ^2\nonumber \\= & {} \dfrac{1}{\beta _p^2} \sum _{{{\varvec{r}}},t} \left\| \alpha (t)\phi ^j \Pi _{\theta (t)} [{\varvec{\nabla }}f({{\varvec{x}}})\cdot {\varvec{\psi }}^j({{\varvec{x}}})]-\rho ^{0}({{\varvec{r}}},t)\right\| ^2 \end{aligned}$$

(8)

Because the above functional couples space and time, it is proposed to use a staggered minimization. The following two steps are to be performed alternately

$$\begin{aligned} {\varvec{\phi }}= & {} \underset{{\varvec{v}}}{\mathop {}\!\mathrm {Argmin}}~\Gamma _{\text {vib-lin}}({\varvec{v}}, {\varvec{\alpha }}) \end{aligned}$$

(9)

$$\begin{aligned} {\varvec{\alpha }}= & {} \underset{{\varvec{a}}}{\mathop {}\!\mathrm {Argmin}}~\Gamma _{\text {vib-lin}}({\varvec{\phi }}, {\varvec{a}}) \end{aligned}$$

(10)

leading to the following linear system, where l indexes the different $N_l$ angles, each composed of $N_t$ projections

$$\begin{aligned} \phi _i = \left[ \sum _{l=1}^{N_l} {\varvec{H}}^l\right] _{ij}^{-1} \left( \sum _{l=1}^{N_l} h^l_j\right) \end{aligned}$$

(11)

where ${\varvec{H}}^l$, indexed by l for angle $\theta _l$, denotes the hessian of $\Gamma _{\text {vib-lin}}$ minimized with respect to ${\varvec{\phi }}$,

$$\begin{aligned} H^l_{ij}= & {} \sum _{{{\varvec{r}}},t_l}\alpha (t_l)^2 \left( \Pi _{\theta _l}\left[ {\varvec{\nabla }}f({{\varvec{x}}})\cdot {\varvec{\psi }}_i({{\varvec{x}}}) \right] \right) \left( \Pi _{\theta _l}\big [ {\varvec{\nabla }}f({{\varvec{x}}})\cdot {\varvec{\psi }}_j({{\varvec{x}}}) \big ] \right) \end{aligned}$$

(12)

and where the second member is

$$\begin{aligned} h_i^l=\sum _{{{\varvec{r}}},t_l} \rho ({{\varvec{r}}},t_l) \alpha (t_l)\Pi _{\theta _l}\left[ {\varvec{\nabla }}f({{\varvec{x}}})\cdot {\varvec{\psi }}_i({{\varvec{x}}}) \right] . \end{aligned}$$

(13)

As the acquisition times are randomly sampled, to solve this problem, it is proposed to resort to a statistical property, namely that, due to the large number of snapshots and steady sample vibration, the full statistics of 3D motion is assumed to be identical for any projection angle. Thus, conventionally, one may set $(1/N_t)\sum _{t_l} \alpha ^k(t_l)^2=1$. Let us emphasize that this condition holds for the 3D displacement but not for residual fields. Indeed, some displacement modes may not be visible on the $\rho$ field if the projection direction is colinear to the displacement, yet this displacement amplitude over time may be normalized to unity. Using this assumption the Hessian matrix relative to one specific projection direction, l, can be written

$$\begin{aligned} H^l_{ij} = N_t \sum _{{{\varvec{r}}}} \left( \Pi _{\theta _l}\left[ {\varvec{\nabla }}f({{\varvec{x}}})\cdot {\varvec{\psi }}_i({{\varvec{x}}}) \right] \right) \left( \Pi _{\theta _l}\big [{\varvec{\nabla }}f({{\varvec{x}}})\cdot {\varvec{\psi }}_j({{\varvec{x}}}) \big ] \right) \end{aligned}$$

(14)

It is interesting to note that these Hessian matrices are exactly the same as static and fully defined P-DVC procedure, as presented e.g., in [34]. Indeed the assumption on the stationary amplitude allows one to treat this measurement as if it were static.

Once the spatial mode ${{\varvec{\Phi }}}({{\varvec{x}}})$ has been computed, the determination of time amplitudes becomes a small and elementary step,

$$\begin{aligned} \alpha (t) = \left[ \sum _{l=1}^{N_l} H^l \right] ^{-1} \sum _{l=1}^{N_l} h^l(t) \end{aligned}$$

(15)

where $H^l$ is the scalar Hessians of $\Gamma _{\text {vib-lin}}$ minimized with respect to ${\varvec{\alpha }}$ for angle l

$$\begin{aligned} H^l= & {} \sum _{{{\varvec{r}}}} \left( \Pi _{\theta _l}\left[ {\varvec{\nabla }}f({{\varvec{x}}})\cdot {{\varvec{\Phi }}}({{\varvec{x}}}) \right] \right) \left( \Pi _{\theta _l}\left[ {\varvec{\nabla }}f({{\varvec{x}}})\cdot {{\varvec{\Phi }}}({{\varvec{x}}}) \right] \right) \end{aligned}$$

(16)

and the second member

$$\begin{aligned} h^l(t)=\sum _{{{\varvec{r}}}} \rho ({{\varvec{r}}},t) \left( \Pi _{\theta _l}\left[ {\varvec{\nabla }}f({{\varvec{x}}})\cdot {\varvec{\Phi }}({{\varvec{x}}}) \right] \right) \end{aligned}$$

(17)

If the temporal amplitude $\alpha (t)$ differs from the one that allowed to compute ${{\varvec{\Phi }}}({{\varvec{x}}})$, both spatial mode and temporal amplitudes are recomputed, until a fixed point is reached.

The modal measurement method is summarized in the Algorithm 1.

Relating vibration modes to PGD modes

Up to this point, the vibration modes, ${{{\varvec{\Phi }}}}^{\text {vib}}$, and the above PGD modes, ${{\varvec{\Phi }}}^{\text {PGD}}$ have no reason to coincide. However, as explained in [35,36,37], 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 last requirement comes from the fact that vibrations modes are orthogonal with respect to a “mass” scalar product. In the FE setting, if ${{\varvec{M}}}$ denotes the mass matrix, then

$$\begin{aligned} {{{\varvec{\Phi }}}}^{\text {vib}}_i\cdot {{\varvec{M}}}\cdot {{{\varvec{\Phi }}}}^{\text {vib}}_j=\delta _{ij} \end{aligned}$$

(18)

Using a PCA decomposition, the metric is just the Euclidean one, and hence PCA modes ${{\varvec{\Phi }}}^{\text {PCA}}$ are orthogonal using the identity matrix instead of ${{\varvec{M}}}$. As a consequence, in general, there is no chance that PCA and vibration modes coincide. However, if one assumes that the resonant frequencies of the vibration modes are well separated, one may use a natural hierarchy of modes: The first one dominates over the second, which itself dominates over the third and so on. Because the first mode dominates in terms of amplitude, both PCA and vibration modes coincide up to a scale factor. This is simply because the displacement field is essentially unidirectional, and hence any vector is co-linear to the basis vector. The metric will only affect the norm of the mode, thus, knowing the mass matrix, it is easy to rescale the PCA mode so that it is of unit norm according to the mass matrix. Then (up to a possible difference of sign) the two PCA and vibration modes will coincide. The second (resp. the nth) mode is also assumed to dominate over higher order ones, which means that the first two (resp. n) modes, be they PCA or vibrations, generate identical vector spaces. However, because the scalar product is different for vibration and PCA, these second (nth) modes cannot be the same, but it is straightforward to orthogonalize the basis vectors ${{\varvec{\Phi }}}^{\text {PCA}}$ with the mass metric. Thus because the sub-spaces of displacements obey the hierarchy property, one may identify the vibration modes from the PCA ones (up to undetermined signs) after a mere Gram-Schmidt orthonormalization process (using the mass matrix).

In our case, the modes that we generate are PGD modes. Even though PGD shares many features with PCA, a natural question to ask is whether the above result may be extended to PGD modes. The answer comes from the basis assumption about mode separation: the vector spaces of vibration modes are again to match those from PGD simply because of the postulated hierarchy. Henceforth, PGD modes can be orthonormalized using the mass matrix and hence coincide with vibration modes. Additionally, one can show that the PGD modes are also orthogonal with respect to the natural scalar product issued from P-DVC, namely with a metric given by the Hessian ${\varvec{H}}=\textstyle \sum _l H^l$ (see Eq. 16).

Let us underline that the assumption on the hierarchy property can be probed without the orthonormalization of modes: the coincidence of sub-spaces truncated at order k can be written as

$$\begin{aligned} \Vert {{{\varvec{\Phi }}}}_k^{\text {vib}}-\sum _{i=1}^k ({{{\varvec{\Phi }}}}_k^{\text {vib}}\cdot {{\varvec{\Phi }}}_i^{\text {PGD}}){{\varvec{\Phi }}}_i^{\text {PGD}}\Vert _{{\varvec{H}}}\le \eta \Vert {{{\varvec{\Phi }}}}_k^{\text {vib}} \Vert _{{\varvec{H}}} \end{aligned}$$

(19)

where $\eta$ is a small parameter $\eta \ll 1$ introduced as a small tolerance, (otherwise the left hand side should be null). $\Vert \ldots \Vert _{{\varvec{H}}}$ refers to the ${\varvec{H}}$ norm: $\Vert {\varvec{A}} \Vert _{{\varvec{H}}}={\varvec{A}}^t{\varvec{H}} {\varvec{A}}$.

It is also noteworthy to observe that the assumption about the complete hierarchical structure of the vibration space, mode after mode, may be relaxed. It suffices to have a large frequency band gap separating the first n modes from the following ones, to expect that the nth dimensional space generated by the first n vibration modes,and n PGD modes, coincide. In the case, the above condition can be relaxed, for all $k\le n$

$$\begin{aligned} \Vert {{{\varvec{\Phi }}}}_k^{\text {vib}}-\sum _{i=1}^n ({ {{\varvec{\Phi }}}}_k^{\text {vib}}\cdot {{\varvec{\Phi }}}_i^{\text {PGD}}){{\varvec{\Phi }}}_i^{\text {PGD}}\Vert _{{\varvec{H}}}\le \eta \Vert {{{\varvec{\Phi }}}}_k^{\text {vib}}\Vert _{{\varvec{H}}} \end{aligned}$$

(20)

(note that the summation over the index i has been extended to n instead of k.)

Application of the method

Synthetic test case: free vibrations of a liver

The method is tested with the identification of 3 vibration modes of a human liver at [3.2–4.1–6.1] Hz. Although artificial, this test case has been designed to mimic the free vibration modes and standard elastography procedures. The test case is generated from the following steps:

• A CT-scan of a liver and the corresponding surface mesh is downloaded from the IRCAD database [38].

  The volume is composed of 512 $\times$ 512 $\times$ 124 voxels with a physical size of [0.56, 0.56, 1.6] mm. Figure 2(a, b) shows a slice of the scan, the 3D numerical mesh and mode 1 amplitude.

• The numerical vibration modes ${{{\varvec{\Phi }}}}_i^{\text {vib}}$ are computed using the meshed volume, a finite element code (Abaqus), and state of the art parameter values taken from the literature. As such these modes are considered as the ground truth, and should be distinguished from the PGD modes ${{{\varvec{\Phi }}}}_i^{\text {PGD}}$, and from the estimated vibration modes ${{\varvec{\Phi }}}_i^{\text {est}}$ computed from the latter.

• The CT image is warped with three numerical modes with realistic random amplitudes, corrupted by noise and projected to generate the experimental sinogram.

The commercial software Abaqus is used to compute the numerical vibration modes ${{{\varvec{\Phi }}}}_i^{\text {vib}}$. 201,627 degrees of freedom are used in the FE model to obtain the modal convergence and linear tetrahedron elements (type C3D4) are used to perform the frequency analysis. The material parameters are set to representative values of soft tissues [39,40,41]. Because only a small vibration amplitude is considered, a simple linear elastic constitutive law is chosen, and should be considered as the tangent behavior in a small strain regime. The elastic modulus is set to 30 kPa, the Poisson’s ratio is set to 0.45 and the mass density is set to 1.59 g/cm$^3$. The first three vibration frequencies are obtained at [3.2–4.1–6.1] Hz.

Fig. 2

a Slice of the liver CT scan and b 3D view of the scan with the fine numerical mesh. The color of the liver is chosen according to the amplitude of the first computed vibration mode ${{{\varvec{\Phi }}}}_1^{\text {vib}}$

Each frame is then generated with the projection of the static volume warped with a combination of the first three computed vibration modes: $p^{\text {est}}({{\varvec{r}}},t_l)=\Pi _{\theta _l}f({{\varvec{x}}} + {{\varvec{u}}}^{\text {vib}}({{\varvec{x}}},t_l))$ with ${{\varvec{u}}}^{\text {vib}}({{\varvec{x}}},t) = \sum _m \sigma _m(t){{\varvec{\Phi }}}_m^{\text {vib}}({{\varvec{x}}})$, ${{{\varvec{\Phi }}}}_m^{\text {vib}}$ being the computed vibration modes and $\sigma _m(t)$ random amplitudes following a normal distribution whose standard deviation is chosen to have a constant energy per mode. In the considered test cases, the vibration amplitude was such that the local strain never exceeded 1.5%, justifying the assumption of a linear elastic model.

The projection is performed using the ASTRA toolbox [31] and a parallel beam assumption (the beam model does not affect the procedure). Five projection angles equally distributed (every $2\pi /5$) are selected. The choice of five projections is motivated by the fact that sensitivity in the different space direction is needed (hence more than one projection angle) and a small number is more challenging for the proposed method. 300 projections of 512$\times$124 pixels are generated for each angle (i.e., total of 1500 projections). Finally, a white Gaussian noise with a standard deviation of 1% of the gray level dynamic is added to each projection.

PGD modes measurements

The kinematics is identified with a coarse boxed mesh: a circumscribed [$3\times 3\times 3$] nodded box with tri-cubic interpolations. This simple mesh was chosen to reduce the number of space unknowns (81 space degrees of freedom) and the spline interpolations are well adapted to low frequency mode shapes (and ensure long-range sensitivity for each node). The two meshes are compared in Fig. 3. In the following, when a field ${{\varvec{\Phi }}}$ is decomposed over the fine (Abaqus) mesh, it will be denoted as $\widehat{{\varvec{\Phi }}}$, whereas the coarse mesh interpolation will be written $\widetilde{{{\varvec{\Phi }}}}$.

Fig. 3

Fine numerical and circumscribed mesh composed of respectively 201,627 and 81 space degrees of freedom (with tri-cubic $C_1$ interpolation). Note that the bottom left node has a non-zero sensitivity as a result of the cubic interpolation

In order to compare numerical and measured modes, the initial Abaqus modes $\widehat{{\varvec{\Phi }}}_m^{\text {vib}}$ are projected onto this coarse boxed mesh and written $\widetilde{{{\varvec{\Phi }}}}_m^{\text {vib}}$. For the different spatial modes considered in this study, the difference between the initial and projected displacement fields was observed to be negligible. A summary of all mode notations is proposed in Table 1.

Table 1 Summary of all measured and numerical mode notations

Two initial residual fields acquired at the same angle are presented in Fig. 4. It can be seen that the residual fields (b and c) show different patterns because the displacement amplitudes are different (they have been acquired at random instants of time, and in the two chosen cases, in most parts of the projected liver, the two displacements point to opposite directions).

Fig. 4

a Projection of the liver at angle $8\pi /5$ and b, c two initial projections corresponding to two different vibration states. The ground truth amplitudes for mode 1–2–3 are respectively for b and c [0.40; 0.76; − 0.19] and [0.32; − 0.20; − 0.23]. Note that projections are shown with the acquisition scale factor [0.56, 1.6]

Two analyses can be performed with the results of this test case:

• The evolution of the projected residual field that shows how the procedure works and converges. This is the true minimized quantity of this procedure.

• The comparison between the measured kinematics and the numerical one used to generate the test. It is interesting to note that the measured modes can be different from the numerical modes yet generate the correct kinematic subspace (e.g., because of the structure symmetry, two modes with the same frequency may be identified up to an arbitrary linear combination of each other.) More than a mode to mode comparison, considering the space generated by the kinematics will allow to validate the procedure.

The identification of the first modes converges in 5–10 iterations. The measured 3D space displacement fields, expressed on the numerical finite element mesh, are shown for the [x, y, z] components Fig. 5. The norm of the residual fields is decreasing with the iterations [3.48–2.04–1.77–1.47–1.38]. The associated amplitudes will be discussed in the next session.

Fig. 5

Measured displacement fields of the five first modes ${{\varvec{\Phi }}}_i^{\text {PGD}}$ in the x, y, z directions (here shown on the numerical mesh yet computed on the boxed mesh composed of 81 space degrees of freedom). The amplitude is to be associated with the 1500 time amplitudes normalized angle-wise

With the kinematics correction, the motion-induced patterns of the residual fields are erased (see Fig. 6). The standard deviation of the residual is presented in Fig. 7. It can be seen that the residual is reduced down to 1.3% (1% is the noise standard deviation), no additional modes are required.

Fig. 6

Residual fields projected at $6\pi /5$ (a) before kinematics correction and (d) after the 5 mode identification. The dashed line corresponds to the envelop of the mesh projection and defines the region of interest. The residual on the top left part corresponds to large motion that may not be correctly captured with the boxed mesh

At this stage, the comparison with the ground truth (space modes and amplitudes) cannot be performed since, as earlier discussed, the PGD modes are not coincident with vibration modes.

Fig. 7

Evolution of the residual norm at each iteration. The residual converges after 4–5 modes meaning that no additional mode are required

Vibration modes analysis

The measured vibration modes can be orthogonalized with respect to the scalar product in order to be compared with the ground truth vibration modes. A Gram-Schmidt procedure is thus applied: the first mode remains unchanged, the following modes are orthogonalized with respect to the mass scalar product with the previous ones. Figure 8 hows the components of ${{\varvec{\Phi }}}_i^{\text {est}}$ on the basis of ${{\varvec{\Phi }}}_i^{\text {PGD}}$. The rather large absolute value of the off-diagonal component do show that the PGD basis is far from being orthogonal with respect to the mass matrix, and hence, even if the vector space generated by the first $i=1,\ldots ,n$ modes, either ${{\varvec{\Phi }}}_i^{\text {est}}$ or ${{\varvec{\Phi }}}_i^{\text {PGD}}$ coincide, it is imperative to carry out the orthogonalization in order to compare ${{\varvec{\Phi }}}_i^{\text {est}}$ with the actual ground truth modes.

Fig. 8

The estimated vibration modes ${{\varvec{\Phi }}}_i^{\text {est}}$ are computed with a Gram-Schmidt orthonormalization procedure from ${{\varvec{\Phi }}}_i^{\text {PGD}}$

In order to appreciate the obtained results, the 5 normalized vibration modes can be projected onto the initial numerical modes used to generate this test case (in fact their projections onto the boxed mesh ${{\varvec{\Phi }}}_i^{\text {vib}}$). To compare the measured PGD modes and the numerical ones, Fig. 9 shows with the same color code as in the previous figure, the components of ${{\varvec{\Phi }}}_i^{\text {vib}}$ on the set of vectors ${{\varvec{\Phi }}}_j^{\text {est}}$. Now it is clear that the three true vibration modes almost coincide with the estimated ones (possibly with an opposite sign, which is purely conventional). The measured 4th and 5th modes are computed to correct for the slight misalignment of the two spaces. It is seen that these two extra estimated modes are orthogonal to the three true vibration modes.

Fig. 9

Projection of the measured vibration modes $\widetilde{{{\varvec{\Phi }}}}_i^{\text {est}}$ with the ground truth vibration modes $\widetilde{{{\varvec{\Phi }}}}_i^{\text {vib}}$

The (x-component of the) first three estimated vibration modes are displayed together with their imposed counterpart in Fig. 10 (first and third column). The central column shows the comparison of the corresponding estimated and prescribed amplitudes for all the 1500 generated cases.

Fig. 10

a Estimated space vibration modes $\widetilde{{{\varvec{\Phi }}}}^{\text {est}}_i$ (1–2–3) shown here projected along the x direction, c Imposed space vibration modes (actually its boxed mesh approximation) $\widetilde{{{\varvec{\Phi }}}_i}^{\text {vib}}$ (1–2–3) also projected along the x direction and b comparison between the 1500 measured and ground truth amplitudes. For this comparison, the space functions were normalized and plotted with the same color scale and the sign of the space functions was chosen to be the same. The x-component was arbitrarily selected to illustrate this comparison, but the y or z components show a similar proximity

Conclusions and discussion

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 set 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 cannot 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 at rest (known) and the modal basis (unknown). Based on a PGD approach, each mode is identified successively using the statistical assumption that the statistical distributions of the 3D displacement amplitudes are identical for each angle of acquisition (valid for a statistically steady excitation, and a sufficiently abundant sampling).

Let us stress here that the present paper does not pretend that any vibrating system in the considered 1–10 Hz frequency range can be studied in a tomograph. The size of the studied object/organ, its visibility in the scanner, its inner contrast, the amplitude of displacements, and thus in turn the excitation and the mechanical properties, etc. are key factors to make such an experiment feasible. Above all, the proposed methodology relies heavily on linear elasticity, so that displacement superposition holds, and no frequency shift with amplitude is present. This restriction is not limiting with respect to biological tissues (or more generally constitutive materials), but rather it requires imperatively that the vibration excitation be small.

The method has been applied to a numerical test case, which is a mock-up of a low frequency excitation of a liver that would be analyzed with a medical CT scanner. Synthetic projection data are generated by deforming a volume with vibration modes associated with random amplitudes. The numerical experiment is performed with 5 angles and 300 projections per angle. The kinematics is composed of three different vibration modes associated to realistic amplitudes. The estimated displacement modes (i.e., spatial field and amplitude), after its orthonormalization with the mass, are compared with the ground truth and allow validating the modal identification procedure.

This approach requires the acquisition of a single reference static image and a collection of projections (here 1500). It is a huge gain in data/acquisition time compared with standard dynamic acquisition methods based for example on the motion phase [18] (that would be complex when dealing with multiple mode identification). This reduction is of utmost importance in material and bio/medical imaging [42].

This proposed procedure is very light to compute as the volume is never updated and the sensitivity fields computed only once. A finer mesh could easily be designed and is not a limit (and could be associated with an elastic regularization [43]). Implementing a multi-scale approach may be crucial when dealing with large displacements, especially because the volume is never advected in deformed configurations. The procedure would be: (i) low-pass filter to smooth the functional shape (convolution with a Gaussian kernel), (ii) downsizing to reduce computation time. This pyramidal identification would be performed from the coarsest to the finest scales.

It is noteworthy that from the Rayleigh-Ritz ratio, a coupling of the obtained results with assumptions on the stiffness $[{{\varvec{K}}}]$ and mass $[{{\varvec{M}}}]$ (given by the tomography though it has to be correctly weighted) matrices enables the frequencies of the measured modes to be identified through

$$\begin{aligned} \left( \omega ^k\right) ^2 = \dfrac{{{\varvec{\Phi }}}^k [{{\varvec{K}}}] {{\varvec{\Phi }}}^k}{{{\varvec{\Phi }}}^k [{{\varvec{M}}}] {{\varvec{\Phi }}}^k} \end{aligned}$$

(21)

Optimizing such a measurement, i.e., choice of angles, surface texture, number of projections per angle, is an interesting perspective. This prior optimization could be performed with respect to an initial numerical model (e.g., an elastic model) from which the kinematics could be identified. Any prior knowledge (for example learnt from other samples as proposed in [44] for the liver shapes) could be combined for a better suited identification.

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 has 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 (minimal exposure time of 0.0125 s). The sample could be soft, made out of gel or silicon to have low frequency modes and the excitation (ideally white) performed with, for example, acoustic excitation, air puffs, mechanical tapping etc. Moreover, with the use of a low absorbing sample, an X-ray beam chopper and an intense beam (i.e., a W-target source) the acquisition time could be significantly reduced, thereby making much faster phenomena accessible.

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 multi-view DIC setup composed of a single camera positioned at different places for the acquisition of snapshots could be designed and is an interesting perspective.

Finally, the identification of material or geometrical parameters could be performed with this modal P-DVC framework, giving access to a novel non-destructive testing procedure.