1.

INTRODUCTION

Grating-based X-ray differential phase-contrast¹ uses the Talbot effect to retrieve additional information about the sample from the X-ray wavefront. Besides the conventional attenuation coefficient, the refractive index decrement and ultra-small-angle scattering as the linear diffusion coefficient² can be obtained. This is achieved by generating a periodic interference pattern by a modulation grating G₁, creating self-images at specific distances. The “intensity” is the local mean of the pattern, the relative magnitude of the modulation is called “visibility”, and the position of the pattern is the “phase”. These quantities are altered by the presence of a sample, where attenuation leads to an overall intensity reduction of the pattern, refraction shifts its lateral position and coherent small-angle scattering (diffusion) reduces the visibility. As the interference pattern is usually too small to be resolved directly, an analyzer grating G₂ is placed in front of the detector to sub-sample the wavefront.³ One of the gratings is moved in small increments to obtain the convolution of the G₂ modulation with the interference pattern at multiple positions. From these data points the three signals transmission, dark-field contrast, and differential phase shift can be retrieved. This procedure is called “phase stepping”. The method was developed with highly coherent synchrotron radiation and brought to laboratory setups by including a third grating G₀ in the interferometer.⁴ Placed between G₁ and a conventional X-ray source, it transforms the latter into many narrow slit sources which are mutually incoherent but produce individual G₁ interference patterns adding up constructively at the detector plane. The combination of G₀, G₁, and G₂ gratings is called a Talbot-Lau interferometer. Fig. 1 shows a sketch of the experimental setup of such an interferometer with inverse geometry,⁵ in which the sample is placed between G₁ and G₂.

Figure 1.

Schematic depiction of the setup. It is an inverse Talbot-Lau interferometer with the sample placed between G₁ and G₂. G₀ and G₁ are close to the X-ray source and consist of a single grating, respectively. The G₂ consists of multiple smaller gratings which are tiled to cover the whole detector. The G₀-G₁ combination vibration creates phase fluctuations globally on the detector. The individual G₂ tile vibration creates additional tile-wise variations. The focal spot movement creates global intensity and visibility fluctuations due to shadowing. The interference pattern amplitude is further reduced by grating movement during an exposure.

The first Talbot-Lau interferometer mounted in a continuously rotating clinical gantry is presented in.⁶ It is a modified, commercial CT platform with 80 cm bore size, operated in a standard clinical scan protocol and with sufficient field-of-view for imaging a human. The sampling of the stepping curve relies on vibrations intrinsic to the system, meaning no explicit phase stepping is performed. These vibrations also cause fluctuations of the interference pattern (intensity, visibility, and fringe phase) independently from the sample properties.

For artefact-free reconstruction, it is required to separate changes of the interference pattern caused by vibrations from those caused by the sample. The goal of the presented work is to model the vibration-induced changes in such a way that this separation is possible.

2.

METHODS

The experimental setup of interest in this work is a commercial clinical gantry platform (Brilliance iCT, Philips) which has been retrofitted with gratings to enable human CT scans giving dark-field contrast. A schematic of the interferometer and the expected modes of vibrations is shown in Fig. 1. The system design is presented in.⁶ The gantry is operated in a continuously rotating manner and the acquisition utilizes the vibrations intrinsic to the system which generate sufficient sampling of the stepping curve to perform phase retrieval.

The general approach for separating the changes of the interference pattern induced by vibrations from those induced by the sample is to derive a parametric model from an air scan, which describes the interference pattern generated by the gratings and its variability. The goal is to minimize the number of parameters related to the actual vibration state, as they need to be determined for each subsequent sample scan due to their limited reproducibility. The process is set up in two steps: In the first step, polynomial variations for phase and visibility are assumed. In a second step, the number of parameters is reduced using a principal component analysis.

2.1

Forward model

The canonical forward model for the measured intensity using a stepped Talbot-Lau interferometer is⁴

with the expected intensity in detector pixel p ∈ {1,…, P} at stepping position index t ∈ {1,…, T}, flat-field intensity I_p, flat-field visibility V_p, flat-field phase 𝜙_ρ, and the global phase shift γ_t induced by moving one of the gratings perpendicular to the grating bars. The flat-fields (I_p, V_p, 𝜙_p) are intrinsic to the interferometer setup and usually assumed to be constant during a scan and between scans. In laboratory setups it is assumed that γ_t is known and exactly the same for air scan and sample scan.

As indicated in Fig. 1 we expect various vibrations which will lead to a pixel- and time-dependent change of the intensity, visibility, and phase of the interference pattern. The forward model (1) is extended to

with , , and representing spatial and temporal fluctuations over p and t.

2.2

Phase fluctuations

As indicated in Fig. 1 we assume translations and rotations of the G₀ grating, G₁ grating, the G₂ carrier, and the individual G₂ tiles. We initially treat each G₂ tile as an independent interferometer, so p refers to a pixel behind one G₂ tile.

According to Dittmann et al.,⁷ the resulting changes of the phase can be accurately modeled by low-order two-dimensional polynomials over the detector. We define a two-dimensional polynomial term 𝓟 = (𝒫_ijp) of order i along the width and order j along the height of the interferometer in pixel p as 𝒫_ijp = w(p)ⁱh(p)^j. w(p) is the w coordinate along the width and h(p) the h coordinate along the height of the interferometer in pixel p. Each term of the polynomial is multiplied with a coefficient γ_ijt to give

A maximum order of one is chosen along the grating bars (in j here) and two perpendicular to them (in i here).⁷

2.3

Visibility fluctuations

According to Horn et al.⁸ change in total phase (i.e. the terms inside the cosine) during an exposure leads to a visibility drop proportional to the first time-derivative of the total phase. Besides movement of the gratings, changes of the X-ray focal spot location can also cause visibility fluctuations as the gratings are bent to a cylindrical surface and carefully focused onto the intended source position. It is assumed that both can be approximated by two-dimensional polynomials. We formulate a general model for the visibility fluctuation

with the maximum polynomial order in i and j doubled in comparison to (3), such that the first time-derivative of the phase vibrations can be modeled.

2.4

Principal vibration components

Instead of modeling each G₂ tile independently, it is desirable to approximate the fluctuations with a small number of dominant shared modes to reduce the amount of free parameters per exposure.

We use principal component analysis (PCA) to reduce the number of parameters. Let X ∈ ℝ^P×T be a data matrix with P variables in the rows and T observations in the columns. The singular value decomposition on X is given as X = UΣV^T, with the orthogonal matrix U ∈ ℝ^P×P, the diagonal matrix Σ ∈ ℝ^P×T, and the orthogonal matrix V ∈ ℝ^T×T. Σ contains the singular values of X on its diagonal which are defined to be in descending order. The rows of ΣV^T are the “principal components” of X and the columns of U are the “principal directions” of X, i.e. the magnitude of each principal component per observation. Let the function PCA be defined as acting on a matrix X and returning ΣV^T and U^T: PCA(X) → (ΣV^T, U^T).

To apply PCA on the combined fluctuations from all G₂ tiles, their terms and are concatenated along the width of the interferometer and the index p is changed from locally on a G₂ tile to globally on the detector. With this definition of PCA, the joint principal components over all G₂ tiles of the visibility and phase fluctuations are determined via

The coefficients and correspond only to the magnitude of principal component k at exposure t, not to the original polynomial model.

The fluctuations can now be approximated with a reduced number B^⋆ and C^⋆ of modes:

The dominant modes of the intensity fluctuations 𝓐=(𝒜_kp)∈ ℝ^P×T are determined by applying PCA on the normalized residuum

with the measured values ŷ_pt and the predicted signal y_pt including only visibility and phase fluctuations. Again we approximate with a reduced number A^⋆ of modes:

3.

RESULTS

We show the results from processing an air scan with the proposed PCA method in Fig. 2. The number of eigenvalues of X^TX before the first “knee” in the scree plots of the principal components in intensity and visibility are used as A^⋆ and B^⋆, i.e. the number of vibration modes to keep for processing a sample scan. Because of the lack of a distinct knee in the phase channel, the number of modes C^⋆ are increased until there is no noticeable difference in the reconstructed refractive index decrement. We obtain A^⋆ = B^⋆ = 2 and C^⋆ = 4. Especially the phase vibrations show tile-wise behavior, while substrate structure is visible in the intensity. The vertical gray lines correspond to the gaps between G₂ tiles.

Figure 2.

The most dominant principal components of the spatial fluctuations in intensity, visibility, and phase (left). All images are scaled to [-1, 1]. Vertical lines correspond to G₂ tile boundaries. The size of each principal component is equal to the number of detector columns times rows. The number of principal components used for processing is determined by the scree plots of the respective (normalized) eigenvalues of X^TX obtained by the PCA (right). Two dominant components in intensity and visibility are identified based on the first “knee” in the scree plots. For the phase channel there is no distinct knee and a number of four principal components are chosen based on the impact on the reconstruction of the refractive index decrement. The vibrations in visibility and phase are a combination of global and tile-wise characteristics.

The main goal of the vibration modeling is to facilitate artefact-free reconstruction. This requires a simultaneous estimation of vibration parameters and sample parameters, which is in general a difficult task and described in a different paper.⁹ Here, we disentangle these two estimations by the following approach for demonstrating the accuracy of the model: The selected sample is small so that only the central part of the detector is covered. This allows us to fit the vibration parameters robustly by a least-squares fit to the outer parts of the detector for each exposure. Subsequently, object parameters are estimated using a sliding-window phase-retrieval.¹⁰

The vibration modes in Fig. 2 are used on the scan of an object consisting of a Polyoxymethylene (POM) cylinder of 5 cm diameter and six falcon tubes with 3 cm diameter each, filled with wool at three different levels of dampness, chocolate chips, water, and neoprene, respectively.

One scan consists of 2400 exposures over a full 360° gantry rotation which takes 1 s. The X-ray tube is operated at 80 kVp and 550 mA. A subsequent air scan is used with the proposed reference processing method to extract the PCA vibration model shown in Fig. 2.

The results are compared with a simplified pipeline in which fluctuations in all channels are only modeled with polynomials 𝒫_ijp up to second order along width and height of the whole detector. Also in intensity, polynomials 𝓟 replace the dominant components 𝓐 from PCA. This simplified pipeline does not handle tile-wise vibrations nor the intricate details of the intensity fluctuations shown in Fig. 2.

Fig. 3 shows the central slices of the resulting volumes in attenuation, diffusion coefficient,¹¹ and refractive index decrement.¹² The reconstruction is performed via filtered back-projection.¹³ The simplified pipeline using global polynomials suffers from artifacts in all channels. The attenuation and linear diffusion coefficient show roughly circular fringe artifacts. The refractive index decrement is dominated by concentric rings and a bright/dark blob structure in the center. None of these artifacts appear in the volumes reconstructed with vibrations from PCA.

Figure 3.

Reconstruction of a phantom in attenuation coefficient (left), linear diffusion coefficient (center), and refractive index decrement (right) with a global polynomial fluctuation model (top) in all channels and principal vibration components obtained via the proposed pipeline (bottom). All channels are free of vibration artifacts when processed with the PCA model as shown in (D), (E) and (F). The materials of the cylinders are (clock-wise starting with the large cylinder) Polyoxymethylene (POM), dry wool, moist wool, soaked wool, chocolate chips, water, and neoprene. No beam-hardening correction is applied which leads to POM and water showing a non-zero diffusion coefficient. The windowing is (A), (D): [-0.05, 2.2] × 10^-1 cm^-1; (B), (E): [-2, 8] × 10^-2 cm^-1; (C), (F): [-4, 2] × 10^-1 cm^-1. Small window in (A), (D): [-3, 3] × 10^-3 cm^-1.

4.

CONCLUSION

We proposed a processing scheme to identify and correct for vibrations of a Talbot-Lau interferometer mounted inside a rotating clinical CT gantry. The tile-wise vibrations are coupled by applying principal component analysis (PCA) and only the first few dominant components are used for processing a sample scan. In the intensity channel dominant fluctuation components are identified by PCA on the normalized residual. The resulting vibration model has few parameters per exposure for the intricate fluctuations, still allowing for an artifact-free reconstruction of a sample. A comparison with a vibration model using global polynomial fluctuations in all channels shows that they are not sufficient to capture the system’s dynamics and lead to artifacts in the reconstruction of a sample.