1.

INTRODUCTION

X-ray imaging is an invaluable technique in medical imaging that enables fast measurements of a sample’s attenuation coefficient. However, it is limited to attenuation based contrast and cannot exploit additional information of the X-ray wavefront like its phase and small-angle scattering. Grating-based X-ray dark-field imaging is capable of measuring all three effects by introducing a Talbot-Lau interferometer in the beam path.^{1, 2}

Due to the added diffraction gratings the X-ray wavefront is modulated to create a reference pattern on the detector. As the sample interacts with the incident wavefront the observed pattern is distorted: Porous materials induce small-angle scattering that smears out the observed interference pattern. A proper analysis of this pattern thus yields additional information about the sample’s microstructure.

The so-called dark-field signal measured by this procedure has proven a promising new tool in the diagnosis of lung diseases as it can show micro-structural changes in the lung parenchyma.^{3, 4}

As a first step in transferring this technology to clinical usage a prototype human dark-field scanner is already in use in the university hospital München rechts der Isar. In a first clinical study it could show the potential of dark-field imaging for the diagnosis of COPD and COVID-19.⁵

However, this system is only capable of measuring chest radiographs. For an unobstructed 3D-view of a patients lung a tomographic setup is needed, as demonstrated for living mice.⁶ Therefore, in a next step we installed a Talbot-Lau interferometer in a clinical CT. The design of this first human sized dark-field CT prototype is presented in.⁷

Previous lab-based dark-field CT systems were designed to meet the stability requirements of the interferometric method. In contrast, our goal was to fit a Talbot-Lau interferometer into an existing clinical CT system with minimal hardware changes. This implies, in particular, a higher level of vibrations than desired, which we address by more complex data processing. Furthermore, the continuous rotation of the gantry prohibits the use of conventional phase demodulation techniques since the detector is positioned differently in each shot.

This work presents the second part of the data processing pipeline of the setup - the processing of sample scans covering the detectors full field of view. The prerequisite reference processing of air scans is discussed in.⁸ A schematic overview of the full processing pipeline is depicted in Fig. 1.

Figure 1.

Processing pipelines for reference and sample scans. a, Reference processing pipeline to extract scan-to-scan persistent system characteristics from an air scan. Here we introduce local intensity fluctuation corrections and information compression using principal component analysis (PCA). b, Interferometer performance and correction array results from reference processing. c, Sample processing pipeline, based on sliding window processing. As the high frequency oscillations differ from the reference scan, we use an optimization step to identify the optimal linear combination of the correction arrays to estimate the current sample free fringe parameters and thus to suppress vibration artefacts. d, Results of the correction coefficient optimization using prior knowledge from the reference scan. e, Coefficients of the three intensity correction arrays. The angular position dependent drift (orange) is scan-to-scan consistent and is utilized as prior knowledge in the correction optimization step. f, Three-dimensional scatter plot of the three intensity correction coefficients after subtraction of the low-frequency component shown in e. A correlation can be observed which also is used during the optimization step as prior knowledge. Figure and caption adapted from.⁷

2.

MODEL

The basic forward model of grating-based phase-contrast and dark-field X-ray imaging is given by¹

where the state of the interferometer is described by the parameters I, V and 𝜙 which denote the mean intensity, the fringe amplitude (visibility), and the interferometer phase, respectively. Once a sample is introduced it attenuates, small-angle scatters, and phase-shifts the incident wave which alters the measured fringe pattern. These changes can be modeled by three sample parameters: the sample transmission T, the dark-field D, and the (differential) phase Φ.

For a data acquisition scheme utilizing a continuously rotating gantry and in the presence of vibrations and other fluctuations, all of the introduced parameters generally depend on the individual detector pixel x and time t. For the presented dark-field CT prototype the effects of the observed fluctuations are described by linear combinations of per-pixel correction arrays R that describe the impact of the main modes of grating deformations on the measured fringe pattern and their respective per-shot strength coefficients c. These perturbations alter the flatfields I₀, V₀, and 𝜙₀ of the three interferometer parameters as follows:

While the flatfields I₀, V₀, and 𝜙₀ and the correction arrays R_i, R_v, and R_p are determined during the reference processing of an air scan (see Fig. 1a, b), the fluctuation strength coefficients of the sample scan c_i, c_v, and c_p have to be determined during sample processing. Subsequently, the sample parameters T, D and Φ can be extracted.

3.

FLUCTUATION ESTIMATION

The observed fluctuations described by c_i, c_v, and c_p are caused by periodic processes happening during the measurement. These include, most importantly, the gantry rotation, the change of the X-ray focal spot due to induced magnetic fields, and vibrations caused by motors on the gantry. While the effects of the gantry movement lead to a reproducible low frequency change in the coefficients, the other effects result in high frequency oscillations. It is therefore possible to extract the low frequency components of the coefficients LF from the reference scan by low-pass filtering the reference coefficients. Only high frequency effects are estimated separately during sample processing.

By analyzing multiple reference scans it was found that the high frequency components of the fluctuation coefficients can be approximated well by a sum of a few sinusoidals. These oscillate at characteristic frequencies of the setup which can be assigned to different components on the gantry. The most important frequencies are two dominant frequencies from the X-ray tube asynchronous motor which drives the anode rotation: the stator frequency is 189 and the rotor frequency approximately 176. Moreover, vibrations induced by the cooling unit at 117 are included in our model.

It was further found that the amplitudes of these oscillations stay constant from scan to scan and can therefore also be extracted from the respective reference coefficient. Additionally, the phase relationship of the oscillations between different coefficients stays constant. These two facts can be illustrated nicely by using pairs of coefficients to plot Lissajous curves which can then be compared between scans. An example is given in Fig. 1f: Here it can be seen that all values of the intensity fluctuation coefficients lie on a well defined Lissajous ellipse. This ellipse stays constant between scans, which indicates a constant amplitude and phase relation of the coefficients oscillations.

Based on this analysis each coefficient c_n is expressed as

where and are the reference amplitude and oscillation phase of the respective reference coefficient and frequency f_k which were retrieved during processing of an air scan. Employing this model has the advantage that only the global phase offsets of each frequency δ_k have to be found for the sample scan. They are extracted via Fourier-analysis and subsequent fitting of initial guesses of the total flux and interferometer phase. These are acquired by conventional processing algorithms. Then, all coefficients can be synthesized from parameters extracted during the reference processing.

4.

SIGNAL EXTRACTION AND RECONSTRUCTION

Once the fluctuation strengths have been estimated, per-shot references can be generated by applying equations (2), (3), and (4). Next, the measured interference pattern has to be demodulated to retrieve the three sample channels T, D, and Φ. This is achieved by first applying a change of variables to equation (1) which is defined by

The forward model is now linear in the new set of sample parameters A, B, and C:

Demodulation is performed by defining patches of pixels on the detector and consecutive shots in the measured sinogram y_meas. For each of these patches a least squares optimization is performed to minimize a cost function 𝒞_patch and find the sample parameters. The cost function is defined as

with weights w that can be used to introduce a distance weighting in the patch to weight pixels and shots less towards the borders of the patch.

To account for a continuously rotating gantry and spatial changes of the sample parameters, in each patch they are represented as linear combinations of known basis functions and their a priori unknown coefficients:

For simplicity the basis functions α, β, and γ were chosen to be polynomials of x and t. Consequently, the model function y_patch assumed in each patch is given by:

This model is linear in the basis function coefficients. Therefore, the least squares optimization is a linear regression and possesses an analytical solution that yields the optimal coefficients:

The full demodulation of the measured sinogram is achieved by iterating over all pixels and at each step processing the patch with the current pixel positioned in its center. The sample parameters at the current position are then given by evaluating equations (10) at the center of the patch using the least squares optimization results.

After performing this patch-wise sliding window phase retrieval the variables are transformed back using:

This yields the three sinograms T(x, t), D(x, t) and Φ(x, t). After transforming T and D into line integrals by applying the negative logarithm function they can be reconstructed using a filtered backprojection. For the dark-field a weighting has to be performed during reconstruction to account for sensitivity differences depending on the sample’s position within the interferometer.⁹ Lastly, post processing steps are applied to enhance the image quality and correct for dark-field beam hardening effects.

5.

RESULTS

The results of applying the processing pipeline to a measurement of a human thorax phantom are depicted in Fig. 2. The phantom was filled with a neoprene insert to simulate the scattering behavior of lung tissue. Additionally, three Falcon tubes filled with dry wool, wet wool, and powdered sugar as well as a POM cylinder were placed inside the neoprene insert.

Figure 2.

Processing results of a measurement of a human thorax phantom. Sub-figure a shows the reconstructed attenuation channel (-1000, 160) [HU] while sub-figure b depicts the dark-field. It can be observed that dry wool, powdered sugar, and the large porous neoprene insert simulating the lung show a strong dark-field signal, while the dense POM cylinder vanishes fully.

The attenuation channel reconstruction appears artefact free. In the dark-field, there are minor streak artefacts, and noise is corrupting the image quality slightly. However, this could be mitigated by employing further post-processing and filter algorithms.

It can be observed that the neoprene insert, dry wool, powdered sugar and bones, show a strong dark-field signal, while the POM cylinder and the soft tissue surrounding the ribs vanish fully in the dark-field reconstruction. This demonstrates the system’s capability to differentiate between porous and dense materials. Moreover, it shows the possibility of analyzing the scattering strengths of the different materials quantitatively.

6.

CONCLUSION

This work presented the sample processing pipeline of the first clinical dark-field CT scanner. The two main processing challenges in translating grating-based dark-field imaging to a clinical CT system are addressed and solutions are proposed.

The first challenge lies in describing the interferometer state accurately for each shot in spite of vibrations. This was addressed by analyzing and adapting measurements of an empty reference scan to the sample measurement.

The second challenge lies in the demodulation of the signal. A continuous tomographic data acquisition shows the sample under a different angle for each projection. This leads to movement artefacts when trying to demodulate the signal with a conventional sliding window approach. To mitigate the effects of sample movement, an advanced patch-wise sliding window algorithm is presented that estimates the sample movement locally using continuous basis functions. Since the proposed model leads to a linear optimization problem, demodulation can be performed fast, while only minor artefacts related to noise are present in the final reconstruction.