Modeling Vibrations of a Tiled Talbot-Lau Interferometer on a Clinical CT

X-ray computed tomography (CT) is an invaluable imaging technique for non-invasive medical diagnosis. However, for soft tissue in the human body the difference in attenuation is inherently small. Grating-based X-ray phase-contrast is a relatively novel imaging method which detects additional interaction mechanisms between photons and matter, namely refraction and small-angle scattering, to generate additional images with different contrast. The experimental setup involves a Talbot-Lau interferometer whose susceptibility to mechanical vibrations hindered acquisition schemes suitable for clinical routine in the past. We present a processing pipeline to identify spatially and temporally variable fluctuations occurring in an interferometer installed on a continuously rotating clinical CT gantry. The correlations of the vibrations in the modular grating setup are exploited to identify a small number of relevant fluctuation modes, allowing for a sample reconstruction free of vibration artifacts.


I. INTRODUCTION
G RATING-based X-ray differential phase-contrast [1], [2], [3] 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 [4] can be obtained. It has been shown that the latter provides signals from structures smaller than the system's direct resolution [5], [6], [7], [8]. A periodic intensity pattern is generated by a modulation grating G 1 , 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 2 is placed in front of the detector to sub-sample the wavefront [3]. One of the gratings is moved in small increments to obtain the convolution of the G 2 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" [3], [9].
The method was developed with highly coherent synchrotron radiation and brought to laboratory setups by including a third grating G 0 in the interferometer [9]. Placed between G 1 and a conventional X-ray source, it transforms the latter into many small slit sources which are mutually incoherent but produce individual G 1 interference patterns adding up constructively at the detector plane. The combination of G 0 , G 1 , and G 2 gratings is called a Talbot-Lau interferometer. Fig. 1 shows a sketch of the experimental setup of such an interferometer with inverse geometry [10], in which the sample is placed between G 1 and G 2 .
The first Talbot-Lau interferometer mounted in a continuously rotating clinical gantry is presented in [11]. It is a modified, commercial CT platform with 70 cm bore size, operated in a standard clinical scan protocol and with sufficient field-of-view to capture 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 spatial fluctuations in intensity, visibility, and fringe It is an inverse Talbot-Lau interferometer with the sample placed between G 1 and G 2 . G 0 and G 1 are close to the X-ray source and consist of a single grating respectively. G 2 consists of multiple smaller gratings which are tiled to cover the whole detector. The G 0 -G 1 combination vibration creates global phase fluctuations on the detector. The individual G 2 tile movements lead to tile-wise phase and visibility 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.
phase independently of the sample properties: vibrations and centrifugal forces from the rotation move the gratings and cause complex spatial phase variations. These also lead to a reduced visibility as the interferometer movement occurs during detector exposures, effectively blurring the modulation pattern. Grating G 2 is realized by tiling multiple smaller gratings to cover the whole detector, all of them vibrating potentially independently. Finally, small movements of the X-ray focal spot on the rotating anode cause defocusing of the interferometer and lead to spatial visibility and intensity fluctuations.
In this paper we propose a method to identify those variations in a Talbot-Lau interferometer. The fluctuations in visibility and phase are modeled per G 2 tile as linear combinations of two-dimensional polynomials. To enable reconstruction of a sample scan, these are reduced to a joint vibration model over all G 2 tiles by employing principal component analysis (PCA) and keeping only few components after factor analysis. In the intensity channel, the dominant fluctuations are determined by applying PCA on the model residual directly. The joint model is translated to a sample scan and enables the estimation of tile-wise vibrations behind a sample by determining the coefficients of the PCA model outside of the sample.

II. PREVIOUS WORK
A Talbot-Lau interferometer as described in section I is sensitive to mechanical instabilities due to the fine grating structures in the order of few micrometers. By including a parameterization of vibrations in the image reconstruction process, which is usually formulated as an optimization problem, artifacts can be avoided. Most common is the estimation of the position of the stepped grating modeled by a global shift of the observed phase on the detector. A popular algorithm based on the model likelihood, alternating between optimization of images and grating shift, has been presented in [12] for phase shifting interferometry and transferred to grating-based X-ray differential phase-contrast in [13] and [14]. Instead of optimizing the likelihood one can also maximize image smoothness [15], [16]. The vibration parameters are determined by a nested optimization, using a linearized least-squares algorithm for obtaining the projections during optimization. The errors resulting from incorrect stepping positions can also be removed in post-processing as proposed in [17]: The authors predict the artifacts from the reference phase image, assuming the errors are uncorrelated with the sample images.
Additional to a uniform phase offset, mechanical instabilities can also lead to relative rotations of the gratings resulting in spatial changes of the phase [18]. The authors propose an alternating algorithm, optimizing the model likelihood with respect to grating positions and pixel-wise flatfields, respectively. If only a transverse shift is considered, finding the vibration coefficients can be done in one step. For movement along the optical axis or rotation of the gratings the optimization is iterative.
While the previously discussed approaches are limited to phase vibrations, in [19] an algorithm is presented to correct for flux variations (i.e. global intensity offsets) together with stepping errors. It is computationally very cheap, but considers only global offsets in intensity and phase.
Concepts to compensate multiple sources of error become increasingly important when translating grating-based X-ray imaging to clinical practice. To scan human patients, a large field-of-view compared to lab-based setups is required. As cardiac and respiratory motion of the patient cause movement artifacts, a short acquisition time is desirable. Combining those requirements results in quickly moving interferometers in such systems, making them more prone to vibration artifacts than their laboratory predecessors. The scanning setup for chest radiography presented in [20] exhibits fluctuations in visibility caused by the interferometer moving during one exposure. A processing algorithm incorporating vibrations both during and between exposures, i.e. variations of visibility and phase, is presented in [21].
The work by [22] uses a tiled, but stationary interferometer and does not model tile-wise vibrations. We show in our comparisons that such a global model is not sufficient for our setup.
Concerning the translation from radiography to computed tomography, a small animal scanner using a rotating gantry equipped with gratings was discussed in [23] and [24]. The system can be operated in "step-and-shoot" mode where only one grating position is measured per rotation angle. The gantry is not rotating continuously, but is precisely moved to each angular position and kept there for phase stepping. Because of the limited scale and rotation speed of the setup the authors observed phase drifts only dependent on the orientation of the gantry.
A larger phase-contrast CT with clinical components was shown in [25]. It is a tabletop setup in which the sample is rotated and the authors do not discuss vibrations of any kind.
In summary, there exist methods for either spatial phase fluctuations or the combination of either phase and visibility or phase and flux offsets. None of the literature above considers spatial fluctuations in all channels simultaneously. Furthermore, we are not aware of previous work on modeling the coupled vibration of a tiled G 2 grating.

A. Vibrations in a Rotating Talbot-Lau Interferometer
The canonical forward model of a stepped Talbot-Lau interferometer without a sample is [9] y simple pt with the expected intensity y simple pt in detector pixel p ∈ {1, . . . , P} at stepping position index t ∈ {1, . . . , T }, flatfield intensity I p , flat-field visibility V p , flat-field phase φ p , 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. The only parameter which is changed between exposures is the global phase offset γ t by controlling one of the grating positions. In laboratory setups γ t can be deliberately chosen. We call this "explicit" phase stepping.
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. The system design is presented in [11]. The gantry is operated in a continuously rotating manner and, given the fact that the tube is not pulsed and the detector is read out at 2 to 4 kHz, it is impractical to implement explicit phase stepping. Instead, the acquisition relies on the vibrations intrinsic to the system which generate sufficient sampling of γ t to perform phase retrieval both for air and sample scans. Therefore the index t refers to both the stepping position as well as the gantry angle as they are changing simultaneously. We use the term "implicit" phase stepping for exploiting the vibrations to sample the stepping curve at initially unknown points. A schematic of the interferometer and the expected modes of vibrations is shown in Fig. 1. We do not use a single grating for sampling of the stepping curve and all gratings are moving simultaneously. The observed phase variation is a compound effect from all gratings.
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 I vib pt , V vib pt , and φ vib pt representing spatial and temporal fluctuations over p and t.
We initially assume that all G 2 tiles vibrate independently and process the measurement behind each tile separately. The following subsections III-A.1 and III-A.2 consider the measurement behind one G 2 grating tile. In section III-A. 3 we discuss how to combine the vibrations of all G 2 gratings in a joint forward model. 1) Phase Fluctuations: As indicated in Fig. 1 we assume translations and rotations of the G 0 -G 1 assembly, the G 2 carrier, and the individual G 2 tiles. According to [18] the resulting changes of the phase can be accurately modeled by low-order two-dimensional polynomials over the detector behind a G 2 tile. We define a two-dimensional polynomial term P i j = (P i j p ) of order i along the width of the detector and order j along the height of the detector in pixel p as with w( p) the w coordinate along the width and h( p) the h coordinate along the height of the detector in pixel p. Here, the pixel index p refers to the area behind a single G 2 tile. Each term of the polynomial is multiplied with a coefficient γ i j t to calculate the amplitude of the local phase offset per exposure t, leading to the phase vibration A maximum order of one is chosen along the grating bars (in j here) and two perpendicular to them (in i here) [18]. The global phase shift γ t from (1) is represented in component P 00 p and coefficient γ 00t .
2) Visibility Fluctuations: Additional to the phase fluctuations we experimentally observe spatial variations in visibility. These originate from movement of the gratings during one exposure and movement of the focal spot. The latter causes defocusing, or more precisely a non-perpendicular incident angle of the radiation onto the grating, referred to as "shadowing".
A change in total phase (i.e. the terms inside the cosine) during an exposure leads to a drop in visibility for that exposure [26]. Here, we use the model introduced in [21] for visibility loss due to comparatively slow vibrations. The phase coefficients γ i j t are regarded as continuous variables in time τ and the vibration is approximated as Let t be the exposure time during which the phase changes. The cosine is averaged over the exposure time t as The factor t 2 /24 originates from the integration of the cosine. Apparently, the space of possible visibility reductions caused by the finite exposure time is given as a sum of pairwise products of the phase displacements, when using (4) for φ vib pt . 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 they 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 (4), such that the squared phase change in (7) can be modeled. The reduction of visibility varying over exposures leads to variations of the noise level and potentially artifacts in dark-field and differential phase, as they are sensitive to the system visibility [27]. However, this effect is not captured by the presented model.
3) Principal Vibration Components: It is desirable to approximate the fluctuations with few dominant shared modes to reduce the number of free parameters, instead of modeling each G 2 tile separately. The tile movement is expected to be correlated because of the shared G 2 carrier. This leads to similar fluctuations γ i j t over some or all G 2 tiles. It also implies correlations in the visibility drop caused by the finite exposure time as we expectγ i j t to be similar over G 2 tiles. It is beneficial to exploit these correlations and approximate the fluctuations in visibility and phase with a reduced model.
One can use principal component analysis (PCA) to find this new vibration representation. Let X ∈ R 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 ∈ R P×P , the diagonal matrix ∈ R P×T , and the orthogonal matrix 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 To apply PCA on the combined fluctuations from all G 2 tiles, their terms V vib pt and φ vib pt are concatenated along the width of the interferometer and the index p is changed from locally on a G 2 tile to globally on the detector.
With this definition of PCA, the joint principal components B = (B pt ) and C = (C pt ) (both ∈ R P×T ) over all G 2 tiles of the visibility and phase fluctuations are determined via The coefficients β kt , γ kt ∈ R T ×T 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 reduced set of vibration components is used for reconstruction in a sample scan. How to determine B and C empirically during air scan processing is discussed in section IV-A.

4) Intensity Fluctuations:
The defocusing of the interferometer caused by displacements of the focal spot introduced in section III-A.2 also leads to modulations in intensity. Instead of modeling the fluctuations with polynomials as in visibility and phase, the vibration components in intensity A = (A kp ) ∈ R P×T with magnitudes α = (α kt ) ∈ R T ×T can be directly determined. They are computed with PCA on the normalized residual with the measured valuesŷ pt and the predicted signal y pt including only visibility and phase fluctuations. Again we approximate I vib pt with a reduced number A of modes: How to determine A empirically is discussed in section IV-A.

B. Air Scan Processing
This section first defines linearized statistical phase retrieval (SPR) and extends it to include spatial variations in intensity, visibility, and phase. Then a minimization around SPR is presented to determine the temporal coefficients of the fluctuation model.
1) Linearized Statistical Phase Retrieval: We first discuss phase retrieval regarding (1), meaning no per-exposure changes except the global phase offset γ t . The flat-field terms (Î p ,V p ,φ p ) are determined by minimizing the pixel-wise weighted least-squares problem with the predicted signal y pt and the measured valuesŷ pt in pixel p at exposure t. The weights 1/ŷ pt are motivated by the measurement being a Poisson process and the variance of the signalŷ pt therefore approximated by the signal itself. The minimization in (15) can be performed directly by linearizing (1) with respect to the flat-field terms (I p , V p , φ p ) [27]. A trigonometric identity and change in variables is applied, and (15) is solved for (Î p ,â p ,b p ) as an overdetermined linear system as T > 3. The initial flat-field terms of (1) are obtained viaV For ease of notation the SPR can be written in matrix form. Let y p = (ŷ pt ) ∈ R T be the vector of T measurements in pixel p, w p = (w pt ) ∈ R T the weight vector with w pt = 1/ŷ pt , and γ = (γ t ) ∈ R T the vector of global phase offsets. The system matrix A is defined as Note that A does not depend on the pixel p. The weighted least squares version of the linearized model from (16) is formulated in matrix notation as with the diagonal weights matrix W p ∈ R T ×T given as W p = diag w p and the solution vector Solving (20) for x p is performed via matrix factorization such as QR decomposition. The linearized statistical phase retrieval can be extended to include vibrations in all image channels as in (2). The columns of the system matrix A p ∈ R T ×3 in pixel p are modified as The weighted least-squares problem is formulated the same as in (20) except that A p now also depends on the pixel p: 2) Optimizing Vibration Coefficients for an Air Scan: To determine the coefficients (α, β, γ ) for a given vibration model (A, B, C), the pixel-wise cost function L p is defined as a weighted sum over the squared residual after phase retrieval and the sum over all pixels p gives the total cost C For a given fluctuation model (A, B, C) the function C(α, β, γ ) is minimized with respect to the vibration parameters giving the optimal coefficients (α,β,γ ). For the firstorder gradient descent algorithm used for optimization (L-BFGS [28]), the gradient of C with respect to the temporal vibration parameters has to be computed. Of special interest is the derivative of L p through the linearized phase retrieval. The gradient of L p with respect to α kt is computed by as the solution vector x p minimizes L p by definition and ∂ L p /∂ x p = 0. Therefore the derivative ∂ x p /∂α kt through the linearized solve for x p is not necessary. The remaining derivative ∂ A p /∂α kt is straightforward. This concept works equivalently for β kt and γ kt .
3) Optimization Strategy: To avoid running into local minima during the optimization, the cost function C (α, β, γ ) is not optimized with respect to all vibrations immediately: we ignore most vibration terms in the beginning and introduce the parameters successively as outlined below. All vibration coefficients are initialized with α = β = γ = 0. Before being introduced, a parameter is not included in the forward model.
We found the linearized phase retrieval to be most sensitive to a global offset in intensity and phase, as modeled by P 00 p with P i j p defined in (3). A convenient method to find approximationsα 00t andγ 00t for the corresponding coefficients is described in [29], involving a PCA on the normalized measured intensitiesŷ pt − E t (ŷ pt ) where E t is the pixel-wise mean over all exposures.
Generally, we first approximate the whole G 2 as one single interferometer and introduce tile-wise terms only in the end. We first find global offsets (including global fluctuations in intensity), then global polynomials, and finally tile-wise polynomials in visibility and phase. After global visibility and phase polynomials, we compute an intermediate approximation of the intensity PCA model A. It contains slight fringe artifacts and has to be blurred, but is required to determine accurate tile-wise vibrations in visibility and phase afterwards. In the final steps the PCA components on all channels are computed, including an un-blurred intensity PCA model.
In detail, the steps of the optimization pipeline are: 1) Computeα 00t andγ 00t as described in [29].

C. Sample Scan Processing
The vibration coefficients at the time of the sample scan are denoted α s , β s , and γ s . They are similar to the air scan but have to be re-determined, because of the high-frequency oscillation which is not correlated to the gantry position [11]. The PCA vibration model (A, B, C) is inherited unchanged from the air scan processing. For the measurement of a sample scan the following forward model is assumed: with y s pt the expected intensity, I s pt the sample's transmission, V s pt the sample's dark-field, and φ s pt the sample's differential phase projected under gantry angle index t on detector pixel p. Although the focal spot displacement causes fluctuations in the interferometer, it is small enough that it can be ignored with respect to these projections through the sample.
The values of α s , β s , and γ s are found with a sample-free region of the measured projections. There, the forward model simplifies to as I s pt = V s pt = 1 and φ s pt = 0. A least-squares optimization yields the optimal vibration coefficients (α s ,β s ,γ s ): with the sample mask m s being 1 in sample-free detector area, and 0 otherwise. With the vibrations known, the sample's projections I s , V s , and φ s are determined by applying SPR on a small window of consecutive exposures. The result is assigned to the central exposure index of the window. In the literature this has been called "sliding window phase retrieval" [30]. It involves the approximation that the sample does not move inside the window. Due to the high-frequency oscillation creating the sampling of the stepping curve inside the angular window [11], our method most closely resembles "sliding window zigzag". The volumes of attenuation coefficient, linear diffusion coefficient, and refractive index decrement are obtained via filtered back-projection [31] of I s , V s , and φ s , respectively.

IV. RESULTS
We show and discuss the fluctuation model (A, B, C) obtained from an axial air scan. It is then used for processing both another air scan and a measurement of a phantom. We compare the results to a simplified processing pipeline which only uses global (i.e. spanning all G 2 tiles), twodimensional polynomial fluctuations in all channels. Finally, we analyze the temporal stability of the PCA model as well as the impact of the size of the sample-free detector area on the reconstruction quality.

A. Principal Vibration Components
We first discuss the dominant spatial fluctuations (A, B, C) determined by PCA. The number of components A and B in intensity and visibility (as introduced in (11) and (13)) to keep for processing a sample scan is chosen by plotting the singular values determined during PCA sorted by their magnitude for each channel. These so-called scree plots [32] and corresponding vibration components are shown in Fig. 2. The vertical columns in the principal components correspond to gaps between G 2 tiles which are ignored during processing. The factors before the first "knee" of the scree plots are considered important and A and B are assigned to this number in the respective channel. With this method we determine The scree plot for the phase channel does not show a pronounced first knee. C is empirically chosen by first determining A and B and then increasing C until no visual difference in the reconstructions is perceived. Here, we determine C = 4. These numbers of components have been constant for the system under investigation.

B. Comparison With a Polynomial Vibration Model
We show the relevance of the proposed reference processing pipeline by comparing reconstruction results using the PCA model (A, B, C) from Fig. 2 to global polynomial fluctuations P (as defined in (3)) in all channels. The used processing pipeline for the sample scans is described in section III-C.
1) Air Scan: We look at the reconstruction of an air scan as a homogeneous, known "sample". Two directly subsequent scans are performed with no sample in the beam path. 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. The second scan is used with the proposed reference processing pipeline from section III-B to extract the PCA vibration model shown in Fig. 2. The first scan is processed and reconstructed with the pipeline in III-C as if it was a sample scan. We use eleven subsequent angular views per sliding-window. The duration between scans was 30 seconds.
The results are compared with a simplified pipeline in which fluctuations in all channels are only modeled with polynomials P i j p as in (3) up to second order along width and height of the whole detector. In terms of the air scan processing in section III-B.3, this equates to performing steps 1-6 with i = j = 0, 1, 2 in visibility and phase. Also in intensity, polynomials P replace the dominant components A from PCA. For sample scan processing and specifically equations (30) and (31), the same global polynomial vibration model is used in intensity, visibility, and phase. Therefore, this simplified pipeline does not handle tile-wise vibrations nor the intricate details of the intensity fluctuations. Even still, it covers all methods discussed in section II.
We show the central slices of the resulting volumes in attenuation, diffusion coefficient, and refractive index decrement in Fig. 3. 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 the PCA vibration model.
2) Phantom Scan: To show the benefit of the proposed PCA model for reconstructing a real sample, we conduct the same comparison to the simplified pipeline explained in IV-B.1 for the scan of a test object. It consists 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. Again a subsequent air scan is performed to obtain the flat-fields (Î ,V ,φ) and, in the case of the proposed PCA pipeline, the dominant vibration components (A, B, C). All scan parameters are identical to those in IV-B.1. The duration between scans was 30 seconds. The mask m s pt is manually determined based on the measured sample.
The central slices of the reconstructions of attenuation coefficient, linear diffusion coefficient, and refractive index decrement are shown in Fig. 4. The artifacts in all channels are similar to the air scan reconstruction in Fig. 3. The insight from looking at a sample scan is the magnitude of the artifacts compared to the sample contrast: in the attenuation channel no artifacts are visible for both pipelines unless a narrow window is chosen. For the linear diffusion coefficient the artifacts are also visible in a large window capturing the whole sample. In the reconstruction of the refractive index decrement the sample is barely visible when using polynomial vibrations in the simplified pipeline. The artifacts already shown in Fig. 3 are larger in magnitude than the sample contrast. When processed with the proposed PCA model, the concentric rings corresponding to the G 2 tiles vanish and all sample features are clearly visible.

C. Temporal Stability
We analyze the temporal stability of the PCA model by performing the sample processing on several air-scans with increasing time spacing to a sample-scan. Five air-scans are recorded and processed with the proposed pipeline. The temporal distance to the sample-scan is 30, 20, ten, five, and one minute(s), respectively. The resulting reconstructions are shown in Fig. 5. None of them show vibration artifacts (as seen in Fig. 4). The reconstructions of the refractive index decrement in Fig. 5(M)-(O) are affected by strong rings artifacts stemming from slow drift of the flat-field reference phase pattern φ p . This drift is in the order of minutes and much slower than the duration of an axial scan (1 s). Therefore the resulting artifacts are perfectly concentric.

D. Impact of Sample-Free Detector Area Size
We analyze the impact of the size of the sample-free detector area for estimating α s , β s , and γ s during the samplescan. It should be noted that the method of estimating vibration parameters using sample-free detector area is chosen deliberately for simplicity in this manuscript. See section III-C for the definition of the sample-scan processing algorithm. The full detector width along all eleven G 2 tiles as seen in Fig. 2 is 545 pixels and 77.2 cm, corresponding to a field-of-view of 42.4 cm. A large sample is simulated by reducing the area used for vibration estimation evenly on both sides of the sample. The resulting reconstruction results are shown in Fig. 6. The image quality in attenuation and diffusion coefficient is largely unchanged down to 10 % relative width, meaning 5 % are used on each side for vibration estimation. The refractive index seems to be most sensitive and deteriorates already at 20 % width in the outer regions and for 10 % in the vicinity of the sample. For 5 % relative width, all image channels suffer from artifacts. Condition number as a measure for orthogonality of the PCA model (A, B, C) inside the sample-free detector area used for estimating vibration coefficients α s , β s , and γ s during a sample-scan. The components are orthogonal (condition number 1) only for full width. Using less detector area, the condition number increases, decreasing the stability of the coefficient fit.
Less detector width leads not only to a reduced number of data points per exposure for the fit, but also reduces the orthogonality of the PCA model (A, B, C). Thus, the linear system of equations that is used to estimate the vibration coefficients becomes less stable. We estimate the stability via the respective condition number of the matrices (A, B, C) limited to the reduced detector area as described above. The result is shown in Fig. 7. A condition number of 1 means orthogonal components and is achieved only for 100 % detector width. With less detector width for the coefficient fit, the condition number increases in all channels. The phase model C is impacted the most. A higher condition number indicates a higher sensitivity to noise of the coefficient fit, matching the trend in Fig. 6.

V. DISCUSSION AND CONCLUSION
The main finding in this work is that by using the proposed PCA model, we can significantly reduce the number of parameters to model the vibration state of each projection compared to previous methods. The approach is generic and may be applied to other setups, although the number of parameters per projection might be adjusted according to the scree plots.
We further propose estimating the parameters on a per-projection basis from sample-free detector area. This simple method was used in order to validate the model and to investigate the number of PCA components needed for a reconstruction without vibration artifacts. We investigate the limits of this approach by using different amounts of the detector area to estimate the vibration coefficients in the samplescan. The smaller the detector size, the more the estimation suffers from noise. The estimation is further deteriorated by increased inherent instability shown by an increased condition number of the models, in particular for the phase coefficients.
The proposed correction works independently of the sample attenuation or phase shift. We determine the inherent interferometer vibrations, which are independent of the sample properties in the forward model.
The air-scan processing in section III-B, including tile-wise vibrations and the PCA method, is robust to noise and shortened rotation times: we can choose a large tube current and if the statistics are still limited by short rotation time, we can use several gantry rotations to increase the number of data points.
For the sample-scan processing in section III-C these options are not applicable, but since only a small number of free parameters has to be optimized we found that the proposed method is robust if the sample-free area is sufficiently large. As shown in the analysis in section IV-D and especially Fig. 6, the method eventually fails for very small samplefree areas. Because of the comparatively constant condition number (Fig. 7) from 10 % to 5 %, we think the reason is low statistics due to the small area used.
More setup-specific methods to determine fluctuation parameters at the time of the sample-scan using more prior knowledge about their temporal progression and correlation may be viable. The sample processing pipeline with sliding-window phase retrieval and filtered back-projection, as explained in section III-C, is only one possibility to solve equation (31) for the sample projections (I s , V s , φ s ) given the vibration model (A, B, C) from PCA. The crucial step is determining the fluctuation coefficients (α s , β s , γ s ) at the time of the sample scan.
The tiles of G 2 have to exhibit some joint vibration such that the number of model parameters can be reduced by PCA. Theoretically, the tiles could be completely independent. However, as the fluctuation is induced by an external vibration and the tiles are mounted on a common holder, similar vibrations can be expected in practice.
The system has to be temporally stable such that the PCA model from an air scan is applicable on a sample scan. An analysis shows that this is applicable for the presented system for a time delay of at least 30 minutes if only attenuation and dark-field images are desired. The system under investigation exhibits slow drift of the pixel-wise phase reference pattern, leading to concentric artifacts in the reconstruction of the refractive index decrement. This effect is independent of the presented method.
The exact order of the optimization steps in section III-B is empirically found and may be different for another setup. We show the detailed algorithm due to the novelty of the experimental platform and the lack of comparison in the literature.
Removing one grating may reduce system vibrations and their impact on image quality. Recent work uses structured anodes [33], [34] and could make the G 0 grating obsolete. However, the fact that the vibration is not the same on all G 2 tiles clearly shows that at least some vibrations cannot be avoided by omitting G 0 . Furthermore, the structured anodes can only be used for very small fan-angles. Reducing the number of G 2 tiles by implementing larger tiles would naturally reduce the complexity and should be a motivation for further advances in the grating fabrication field.
In summary, we propose a processing scheme to identify and correct for vibrations of a Talbot-Lau interferometer mounted inside a rotating clinical CT gantry. All channels show spatial fluctuations changing between exposures which have a specific shape over the detector. Polynomials in visibility and phase over the whole detector and per G 2 tile are used in the data model and an implicit optimization scheme with linearized phase retrieval as an internal layer is employed to determine their individual magnitude per exposure. The resulting tile-wise vibrations are coupled by applying principal component analysis (PCA) and keeping only the first few dominant components for processing a sample scan. In the intensity channel dominant fluctuation components are identified by PCA on the normalized residual. The vibration model for the intricate fluctuations allows an artifact-free reconstruction in presence of a sample with only a small number of parameters. A comparison with a vibration model using global polynomials shows that the latter is not sufficient to capture the system's dynamics and leads to artifacts in the reconstruction of a sample (see Fig. 3 and 4). The PCA model is temporally stable up to at least 30 minutes. The method of determining vibration coefficients with sample-free detector area limits the usable field-of-view of the scanner. The proposed algorithm can be used to identify setup-specific vibrations in a clinically relevant platform where mechanical instabilities can not be avoided.