Open Access
Add to BibSonomy
Abstract
X-ray grating interferometry CT (GI-CT) is an emerging imaging modality which provides three complementary contrasts that could increase the diagnostic content of clinical breast CT: absorption, phase, and dark-field. Yet, reconstructing the three image channels under clinically compatible conditions is challenging because of severe ill-conditioning of the tomographic reconstruction problem. In this work we propose to solve this problem with a novel reconstruction algorithm that assumes a fixed relation between the absorption and the phase-contrast channel to reconstruct a single image by automatically fusing the absorption and phase channels. The results on both simulations and real data show that, enabled by the proposed algorithm, GI-CT outperforms conventional CT at a clinical dose.
© 2023 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
1.1 Grating interferometry
X-ray grating interferometry (GI) is a promising imaging technique which could extend the limits of existing clinical X-ray imaging by sensing the X-ray beam refraction and small-angle scattering in addition to the beam attenuation [1].
The interaction of X-rays with matter is described by the complex index of refraction $n = 1 - \delta + i\beta$. The real part $\delta$ describes the beam refraction, whereas the imaginary part $\beta$, the attenuation. The small-angle scattering originates from the refraction at many interfaces of fine unresolvable features and is determined by the phenomenological dark-field extinction coefficient $\epsilon$. Conventional attenuation-based CT utilizes only a part of the available information: the beam attenuation. On the contrary, a technique such as GI, which is able to extract more information per measured photon, bears great potential to play an ever-increasing role in the field of X-ray imaging. Refraction and small-angle scattering add complementary and potentially diagnostically relevant information that cannot be accessed with conventional X-ray imaging. Higher soft-tissue contrast can be obtained in phase compared to absorption [2,3], whereas the dark-field channel can highlight strongly scattering structures which often are not visible in the other two channels [4].
Breast imaging is one important application which is speculated to become the first to benefit from the additional information GI provides, since the phase-contrast channel could improve image quality for high spatial frequency details [5] and the dark-field channel could allow for non-invasive microcalcification classification [6]. In fact, to date all available breast imaging modalities come with some limitations [7]. Mammography involves painful breast compression, is intrinsically limited to planar imaging and offers little contrast, especially in highly dense breasts due to tissue overlap [7]. Breast MRI offers superior contrast but has poor resolution and is expensive [7]. Recently, it has been shown that breast tomosynthesis [8] and attenuation-based breast CT [9] can provide high resolution with little or no tissue overlap, respectively. Unfortunately though, the images are characterized by relatively modest contrast-to-noise ratio (CNR), especially for high-frequency details [7]. To tackle these limitations, GI has found increasing application in the field of breast imaging in the last decade [10–12]. Notably, a first-of-its-kind ex-vivo reader study has shown that planar GI imaging yields complementary and diagnostically relevant information [10]. Given the promising results in 2D GI-based breast imaging, our group has embarked upon a long-term effort to build a GI breast CT (GI-BCT) prototype scanner to extend the technology to the third dimension to combine the advantages of improved contrast with eliminating tissue overlap and the need for breast compression [13].
X-ray detectors only measure the intensity of the beam, which is directly related to the beam attenuation. The beam refraction and the small-angle scattering cannot directly be sensed and hence need to be translated into intensity variations that can be recorded with X-ray imaging detectors. To retrieve these additional signals with maximal sensitivity, Talbot(-Lau) GI is utilized with two (or three) gratings, depending on the X-ray source being used. With highly coherent sources such as synchrotrons or microfocal X-ray tubes [14], a phase-shifting grating G1 and an analyzer grating G2 suffice [1], whereas with standard X-ray tubes a third grating, called source grating or G0, must be placed in front of the source to achieve sufficiently high beam coherence [15] (see Fig. 1). The G1 grating, which is designed to diffract with as little absorption as possible, imprints a periodic pattern into the X-ray beam that leads to a self-imaging effect called Talbot effect, which allows to measure the sample-induced distortions of the pattern to analyse the attenuation, refraction and small-angle scattering properties of the sample [16].
Fig. 1. Principle of grating interferometry. The source illuminates G0 which leads to beamlets. G1 creates an interference pattern downstream. The detector pixels are significantly larger than the pitch of the interference pattern and hence an analyzer grating G2 (with a pitch as big as the interference pattern) is needed to translate the pattern into intensity variation on the detector. When the sample is in the beam the pattern is distorted. Refraction by an angle $\alpha$ causes a displacement of the fringes in the direction perpendicular to the grating lines. Attenuation corresponds to a change in intensity, dark-field to a reduction in the pattern amplitude. The G0 and G2 gratings are made from highly absorbing gold structures displayed in yellow. The G1 grating does not significantly reduce the photon flux. The refraction angle magnitude is unrealistically large for clarity. Moreover, while the X-rays would be bent outwards for objects with $n<1$, we displayed inward bending ($n>1$) for clarity.
When a sample refracts the beam the wave’s velocity changes according to $\delta$, which results in a phase shift of that particular wave. Thereby, the wavefront changes direction, i.e., is refracted according to
Fig. 2. Data acquisition in grating interferometry. The five blue dots indicate the photon counts at the five phase steps which are fitted with the blue sine curve from which the intensity map $I_0$, the visibility map $V_0$ and the phase map $\varphi _0$ of the reference scan without sample can be obtained for each detector pixel. The corresponding phase-stepping curve with the sample in the beam is shown in orange. The amplitude of the reference curve (blue) is given by the visibility of the system (9.4$\%$ in our case). By analysing the reduction in mean intensity, the displacement and the loss in visibility of the sample compared to the reference phase-stepping curve, the three image signals can be extracted.
Small-angle scattering, also known as the dark-field signal, reduces the amplitude of the interference pattern, and thus the amplitude of the measured phase-stepping curve (see Fig. 2). Analogous to attenuation, the dark-field signal follows the Beer-Lambert law [17]:
The tomographic forward operators related to the three GI signals are thus given by:
$\mu, \epsilon$ and $\delta$ are the attenuation, dark-field and phase images, $A_\mu, A_\epsilon$ and $A_\delta$ are the system matrices, $T$ denotes the transmission sinogram, $D$ the dark-field sinogram and $\varphi$ the differential phase-contrast (DPC) sinogram. Please note that the system matrix $A_\delta$ includes a differential term [18].As in conventional CT, we are interested in the change in the signals induced by the sample, which requires the projection data to be flat-field corrected [19]. The flat-field scan yields phase-stepping curves for every pixel, which can be characterized by the intensity map $I_0$, i.e., the average of the curves, the visibility map $V_0$, i.e., the amplitude of the curves, and the phase map $\varphi _0$, i.e., the phase of the curves. Arguably, the most critical parameter in GI is the visibility, which characterizes the contrast of the interference pattern in the Talbot carpet. With higher visibilities, the sensitivity to the phase and dark-field signals increases, and therefore the easier the separation of the three signals becomes.
2. Methods
2.1 GI-BCT prototype scanner and measurement protocol
Our GI-BCT setup [13] consists of a symmetric Talbot-Lau interferometer set to the 5th Talbot order. The gratings were bent inwards with in-house holders. The gratings have a pitch of 4.2 $\mathrm {\mu }$m with a duty cycle of 0.5, and G1 was designed to induce a $\pi$ phase shift at 46 keV. We used a single G0 and G1 grating, and an array of three G2 gratings. The G0 and G2 gratings had gold lamellae electroplated onto a graphite substrate. The G1 phase grating was manufactured on a double side polished 8-inch silicon wafer by deep reactive ion etching. The grating lines were 59 $\mathrm {\mu }$m thick and the silicon substrate thickness was 240 $\mathrm {\mu }$m [13]. The source-G0 distance is 100 mm, the source-detector distance 1756 mm and the rotation centre is 1003 mm away from the source. The setup is shown in Fig. 3.
Fig. 3. GI-BCT prototype scanner. From left to right: the source, the source grating, the phase-shifting grating, the scanned breast sample, the three analyzer gratings and the detector.
A Comet MXR-225HP/11 tungsten anode tube with a focal spot size of 250 $\mathrm {\mu }$m was operated at 70 kVp with an Aluminium filter of 3 mm thickness. The tube current was 10 mA or 2.5 mA, depending on the radiation dose of the scan. The projections were acquired with a photon-counting detector with 750 $\mathrm {\mu }$m-thick CdTe sensor and 75 $\mathrm {\mu }$m pixel size, manufactured by Dectris AG, Switzerland. We acquired 1200 projections with a frame-rate of 20Hz for all scans, while the sample was rotated 360$^{\circ }$ at 1 rpm. To reduce the computational cost of the reconstructions, a binning of 2$\times$2 was performed in the projection domain for all scans, thus yielding a pixel size of 150 $\mathrm {\mu }$m.
The simplest way to reconstruct the three tomograms is to acquire multiple exposures at every projection angle (at least three) to build an explicit phase-stepping curve, as shown in Fig. 2 for each pixel in each projection. Fourier analysis or sine curve fitting can then be used together with the flat-field maps to retrieve the three sinograms $T, D$ and $\varphi$ followed by tomographic reconstruction of the three separate volumes [20]. Unfortunately, this approach is unlikely to be compatible with clinical dose and scanning time constraints, requiring the scan to be performed as fast as possible with as few projections as possible. In a clinical environment, this could be achieved with a continuously rotating gantry without phase-stepping.
To address some of these challenges and enable the reconstruction of the three signals from a single measurement, which does not require phase-stepping, the authors in [21,22] proposed the following model-based optimization problem that incorporates the full GI forward model, thus operating directly on the detected photon counts $I$:
Note that according to the central limit theorem the Poisson distribution, which governs counting noise statistics, converges to the normal distribution when large numbers of photons are being detected. In our case, this assumption is satisfied and therefore this cost function assumes a Gaussian noise distribution. Unfortunately, this optimization problem only converges under very favourable conditions with little noise and a very high visibility [21]. In particular, the convergence of $\delta$ is challenging due to the highly ill-conditioned nature of $A_\delta$ [18] and the non-convexity of the cosine term. The convergence of $\epsilon$ is cumbersome as there is little signal compared to the other two channels. On the contrary, $\mu$ tends to converge even under highly noisy conditions since the attenuation signal is strong and $A_\mu$ is significantly less ill-conditioned.
The approach introduced by [21] has recently been extended to helical geometries and a static grating configuration [23]. While the authors could show convergence on simulated data with high counting statistics, the applicability of the approach to a clinically realistic dose regime has not been demonstrated to date. It can be expected that it will be challenging to make this approach compatible with clinical conditions at state-of-the-art GI setups due to the current limitations in visibility and the highly ill-conditioned nature of the problem.
It is worth noting that alternative approaches for solving this highly ill-conditioned problem have been proposed in the literature, most notably the reverse-projection approach [24] and its generalization [25]. However, they assume a linear dependence on the differential phase-contrast signal instead of the cosine relation, which is only valid for small refraction angles and small FOVs.
It is crucial to note that in GI we are sensing the first derivative of the phase, while we aim to reconstruct the refraction coefficient, which is directly proportional to the phase itself. This has important consequences for the noise properties of the reconstructed tomograms. As shown in [5], the noise power spectrum (NPS) in phase-contrast CT behaves differently than in absorption. In absorption the noise is concentrated in the high frequencies, whereas in phase the noise is predominantly in the the low frequencies. Moreover, it been has shown that the CNR in the two channels is a function of the spatial frequency. While the phase channel, which depends on the quality of the GI setup, tends to be superior at higher frequencies, the absorption signal excels in the lower frequencies [5]. In this work we propose to leverage on this prior knowledge to address the limitations of the existing reconstruction methods listed above. In particular, we propose a fused parameterization in which we reconstruct a single image representing both the absorption- and the phase-contrast channel.
If phase-stepping is to be avoided, and therefore a single projection per angle is to be acquired, the implicit PSC has to be sampled at different positions, i.e., with different flat-field phases, for different projections and/or pixels. By sampling the measurements more inhomogeneously, the ill-conditioning of the problem is reduced and the sensitivity to the phase signal is increased. There are essentially two ways an inhomogeneous phase map can be obtained. One can either introduce variations within the projection (from pixel to pixel) or from projection to projection. In the first case the gratings need to be misaligned to create fringes as it was done in [11,23]. This allows to keep the gratings fixed during image acquisition and is thus favourable from a mechanical point of view. In the second case the gratings are aligned to yield a homogeneous phase map and are then moved while the sample is rotating. Since the latter is less favourable from a mechanical point of view, we employed the former case.
We performed a simulation study where we optimized for the fringe angle and period. We varied the two parameters to simulate projection data, reconstructed these data with the proposed algorithm and quantified the mean squared error (MSE) of the reconstruction iterates to the clean ground truth. We found the period of the fringe in the horizontal direction to be the critical parameter (see Fig. 4). Consequently, various combinations of fringe angle and period gave rise to optimal reconstructions (see Fig. 4): in particular, larger fringe angles coupled with larger periods and small fringe angles coupled with smaller periods.
Fig. 4. Fringe analysis. The three plots show experiments with different fringe periods. The fringe angle is plotted along the y-axis, the iteration number along the x-axis. We varied the fringe angle from 0$^{\circ }$ to 45$^{\circ }$, where 0$^{\circ }$ refers to the horizontal axis of the detector, and the period from 20 to 180 pixels. The lower the MSE, the better the result. As expected, the MSEs rise again after reaching the minimum once the algorithm starts fitting the noise.
We would like to mention that decreasing the horizontal period of $\varphi _0$ to an arbitrarily small value does not automatically lead to an increasing image quality, as can be noted by looking at the 20-pixel-period/45$^{\circ }$-angle simulation, which leads to the worst results. In fact, when the horizontal period becomes too small, the reconstructions become affected by strong ring artefacts because the algorithm overfits to the fringes themselves by not being able to distinguish the contribution of the absorption and the phase signal. Still, an optimal horizontal period for which the DPC sensitivity is sufficiently strong, but artefacts are not yet significant, could be found. Since it was empirically easier to obtain a regular fringe pattern with small periods, we opted for a small fringe angle of $\sim$5$^{\circ }$ and a period of $\sim$20 pixels as shown in Fig. 5.
Fig. 5. Flat-field maps measured on our prototype scanner. From top to bottom: intensity map $I_0$, phase map $\varphi _0$, visibility map $V_0$ of the reference scan without sample and profile through central line of phase map. The visibility map has an average of 9.4$\%$. In the phase profile one can observe the horizontal gradient which is the critical parameter for determining the fringe orientation and period. The two vertical lines are caused by the inter-grating gaps.
2.2 Proposed algorithm
We circumvent the ill-posedness of (6) by incorporating a strong prior on the signals, and in particular on their inter-relation. Since there is an underlying structure (the distribution of the sample’s complex index of refraction $n$) that gives rise to the distributions of both $\mu$ and $\delta$, it is realistic to assume a function $g_\theta$ exists that describes the relation between the two channels, i.e.,
where $\theta$ is a set of parameters that characterize the function $g$ and which could either be learned from data or hard-coded based on prior knowledge. In this work we chose to parameterize $g_\theta$ as a piece-wise linear function based on tabulated $\mu$ and $\delta$ values [26–28]:We would like to mention that we could have inverted the direction of the $g_\theta$ mapping, i.e., from $\delta$ to $\mu$. However, this parameterization leads to the reconstruction of an image with phase-like contrast. To make the comparison to attenuation-based CT fairer, we opted for a mapping from $\mu$ to $\delta$ in which an image with attenuation-like contrast is reconstructed. Please note that we assume that a first-order approximation for the X-ray spectrum, i.e., an effective energy, can be found for both attenuation and phase, and that this effective energy is the same for both channels. In its current implementation, we did not model beam hardening in our algorithm. These two approximations constitute a limitation of our method. In fact, ideally one would like the forward operators to model distinct attenuation and refraction behavior depending on the X-ray energy spectrum and to let the mapping become a function of both $\mu$ and energy, i.e., $\delta (E) = g_\theta (\mu (E),E)$. A comprehensive modelling of the energy dependence would result in a significant increase in computational complexity and related reconstruction time which we considered unnecessary at this point.
With the assumption of (7) our optimization problem can be expressed as
Finally, we reformulate the loss function by eliminating the negative exponential term, as this alternative formulation improves convergence speed:
This optimization problem is significantly easier to solve compared to (6) since the number of variables is smaller and the problems caused by the bad conditioning of $A_\delta$ are avoided as $\delta$ is now directly tied to $\mu$, effectively reducing the number of local optima in (6) induced by the cosine term. Moreover, assuming this fixed relation practically means that we aim to reconstruct an image based on more information compared to conventional CT, i.e., not only based on its radon transform but also on a 1-dimensional derivative of its radon transform.
To demonstrate the practical value of our reconstruction approach we will show that the addition of gratings into the beam, combined with our prior knowledge of the relation between $\delta$ and $\mu$, can lead to superior image quality than classical attenuation-based CT. We will thus compare the optimization problem in (10) to
To compare the performance of the two methods in the fairest possible manner both optimization problems were solved with the same solver and with the same tomographic operators. We used the L-BFGS algorithm with line search [30], as it led to rapid and robust reconstructions. The system matrix $A_\delta$ has been implemented based on the derivative of spherically symmetric blob functions as in [18]. $A_\mu$ and $A_\epsilon$ were identical and were also based on spherically symmetric blob functions. However, since these forward operators do not contain a derivative, we used the blob function itself, instead of its derivative. Finally, to avoid confounding effects, no regularization has been used other than early stopping (and other than the $g_\theta$ mapping in the proposed algorithm). All codes have been implemented in Tensorflow [31] or CUDA [32].
2.3 In-silico experiments and measured specimen
An in-silico breast phantom described in [33] was used to validate and compare the performance of our algorithm against conventional CT in a controlled environment. While the grey values are different between the absorption- and phase-contrast phantoms, since attenuation and phase shift coefficients have been assigned based on tabulated $\mu$ and $\delta$ values [26–28], the underlying image content is identical.
We simulated the acquisition protocol described above making use of the experimental flat-field data obtained on our scanner (see Fig. 5) by applying the GI forward operator (see right part of (6)) and then adding Poisson noise. Our simulation did not model the point-spread-function of the detector, nor the electronic readout noise. To avoid using the same operators for the simulation and the reconstruction, the forward projections have been implemented with the ASTRA toolbox [34]. We varied the experimental conditions along two main parameters: dose and visibility.
To simulate a given dose, we matched the photon counts in our simulation to the counts of an experimentally obtained scan, for which dose calculations have been performed [13]. This was possible since our phantom matched the sample of the real scan both in terms of dimensions as well as in terms of tissue composition. Based on these calculations we then simulated acquisitions at $\sim$4mGy, $\sim$8mGy, $\sim$16mGy, $\sim$24mGy, $\sim$48mGy and $\sim$96mGy. As mentioned above, we used twice the amount of photons for the CT-equivalent reconstruction case as the G2 grating absorbs half of the photons [13].
In terms of visibility, we investigated three cases. First, we used the measured visibility map on our prototype (average visibility of 9.4$\%$). We then scaled the map by a constant in order to obtain an average visibility of 17.8$\%$, which, according to wave-propagation simulations, corresponds to the maximum visibility achievable with the current gratings on our scanner. This upper limit is given by grating fabrication accuracy. Finally, we simulated a case in-between with an average visibility of 13$\%$.
Subsequently, to validate the performance of the proposed method on real data, we scanned a formalin-fixed breast specimen provided by the Department for Pathology and Molecular Pathology at the University Hospital Zürich without any significant pathological findings. The specimen has been obtained with written informed consent under the ethical approval KEK-2012_554 obtained from the Cantonal Ethics Commission of canton Zürich. The specimen was scanned with an average visibility of 9.4$\%$ and different doses (5.5mGy, 11mGy, 16.5mGy, 22mGy, 44mGy and 66mGy).
3. Results
We will first focus on the $\sim$24mGy in-silico experiment to comment on some general properties of the FIBIR algorithm. Subsequently, we will include all dose simulations and draw conclusions on the conditions under which the proposed method outperforms the CT-equivalent. Finally, we will show that the reconstruction results of a real breast specimen confirm the results of the simulation study.
A crucial characteristic of the FIBIR algorithm is that the absorption and DPC contributions add complementary information and work symbiotically towards the reconstruction of the single fused image. We can visualize this by looking at the gradients coming from the differential phase part and the attenuation part. Let us define the loss function in (10) as $\mathcal {L}(\mu, \epsilon )$. We can then compute its derivative with respect to $\mu$:
Fig. 6. Gradients of FIBIR and the CT-equivalent on an in-silico phantom. From left to right: gradient of the attenuation term, gradient of the (differential) phase term, total gradient of the fused method, and gradient of the CT-equivalent. From top to bottom, the visibilities were 9.4$\%$, 13 $\%$ and 17.8$\%$. A dose of $\sim$24mGy was used for this experiment. The gradients are shown for iteration 7. This iteration number was chosen since at the first iterations the gradients only carry very coarse signal information while for the later iterations, when the image has almost converged, the gradients do not contain much signal anymore.
This complementarity of information leads to faster algorithm convergence as shown in Fig. 7. In fact, it can be seen that while the low frequencies appear both in classical CT reconstruction and in FIBIR from the beginning, the high frequencies appear early on only in the latter method. Therefore, the proposed method delivers the full range of spatial frequencies with significantly less iterations than the classical CT approach. Furthermore, in agreement with Fig. 6, a higher visibility amplifies the superiority of the proposed approach. The rings which are present during the first iterations of the reconstruction are caused by the impossibility to discern the contribution of the absorption and the phase signal to the measured data during the initial part of the reconstruction. We would like to note though that they disappear once the algorithm moves towards convergence.
Fig. 7. Convergence analysis on an in-silico phantom. On the top, from left to right: iterates for the CT-equivalent reconstruction, for FIBIR with 9.4$\%$ average visibility, 13$\%$ average visibility and 17.8$\%$ average visibility. The proposed algorithm leads to faster convergence and the high frequencies appear earlier with higher visibilities. Note that the noise in FIBIR is higher compared to classical CT since half of the flux reaches the detector. On the bottom: reconstruction error. The CT-equivalent is plotted in orange. FIBIR results are plotted in blue for the measured visibility map (average visibility of 9.4$\%$), in purple for 13$\%$ visibility and in green for the maximum achievable optimal visibility average of 17.8$\%$. A dose of $\sim$24mGy was used for this experiment.
The lower part of Fig. 7 shows the mean squared error (MSE) of the iterates compared to the clean ground truth. We see that FIBIR indeed reaches its optimum after less iterations compared to the CT-equivalent reconstruction. Since the FIBIR algorithm involves three forward operators, whereas classical CT reconstruction only involves one operator, the runtime of the FIBIR algorithm is higher compared to the iterative CT-equivalent. Currently, the computation time for one FIBIR iteration is around 45 s, whereas its classical counterpart employs roughly 15 s. The fact that the MSEs rise again after reaching the minimum is explained by the fact that the algorithms first fit the signal, but later start to overfit to the noise. The CT-equivalent reaches a slightly lower MSE error as it operates on twice the amount of photons.
Since GI-CT has the biggest potential in surpassing conventional CT in the high spatial frequencies [5], we compared the proposed FIBIR algorithm with the CT-equivalent reconstruction in terms of spatial resolution at a given image quality, i.e., at a given CNR. In particular, we smoothed all reconstructed images with a Gaussian kernel until we reached a CNR of 5 to satisfy the Rose criterion [35]. On these images we then performed a Fourier ring correlation (FRC) analysis [36] to estimate the spatial resolution of the images based on the commonly used threshold of 0.143. The reconstruction at the highest dose of $\sim$96mGy has been used as a reference for the FRC.
Figure 8 shows the reconstruction results at a dose of $\sim$24mGy. We can confirm that the proposed parameterization leads to superior spatial resolution (see also Fig. 11) at a CNR of 5 and that higher visibilities positively correlate with higher resolution. Also, from a qualitative point of view, we see that FIBIR leads to sharper high-frequency details (see red arrows in Fig. 8), which is further confirmed by the line profiles in Fig. 9 taken in the region where the red arrows in Fig. 8 point.
Fig. 8. Reconstruction results at $\sim$24mGy on an in-silico phantom. On the top, from left to right: CT-equivalent: FBP, iterative and ground truth. On the bottom: FIBIR with 3 different visibilities (9.4$\%$, 13$\%$ and 17.8$\%$). The spatial resolution was calculated using Fourier ring correlation. Stripe- [37] and ring-removal [38] have been applied to all data sets in the sinogram space and in the image space, respectively. This filtering step did not degrade image quality.
Fig. 9. Line profiles of the reconstruction results in Fig. 8. On the left: line profiles for the FIBIR algorithm with different visibilities. On the right: line profiles for the CT-equivalent reconstruction, using FBP and iterative reconstruction. The line profile of the ground truth image is shown in red. We see that FIBIR yields a sharper profile compared to the CT-equivalent.
So far, we have only considered the case of static gratings. However, we would like to point out that our algorithm is also compatible with phase-stepping GI measurements. When dose-equivalent phase-stepping data is simulated, the gradients contain more high frequency details compared to one-shot acquisitions as can be appreciated by comparing Fig. 10 to Fig. 6. Nonetheless, one can also note that the gradients are noisier because with phase-stepping data each measurement has a significantly reduced number of photon counts for a scan at an equivalent dose. At very low counts, this will increase statistical phase wrapping [5] which degrades the quality of the FIBIR gradients. An advantage of using phase-stepping data is that no rings appear in the reconstruction as it becomes easier to discern the two signals and hence overfitting to the wrong component is eliminated and thus no rings appear.
Fig. 10. Gradients of FIBIR and CT-equivalent with phase-stepping data on an in-silico phantom. From left to right: gradient of the attenuation term, gradient of the (differential) phase term, total gradient of the fused method, and gradient of the CT-equivalent. From top to bottom, the visibilites were 9.4$\%$, 13$\%$, and 17.8$\%$. A dose of $\sim$24mGy was used for this experiment. The gradients are shown for iteration 7.
Figure 11 summarizes the results of the simulation study (left) and the real data (right). It shows which spatial resolution can be achieved at a given dose with the FIBIR algorithm and the CT-equivalent, at a CNR of 5. The results for static, fringe-based scans are shown as solid lines, the results obtained with phase-stepping data are shown as dashed lines. The shade of the line is proportional to the visibility. As expected, the curves show that the iterative CT-equivalent reconstruction performs better compared to FBP at low doses and that the two converge at very high doses. If we first focus on the static scans we see that FIBIR performs better with increasing visibility. Moreover, it reveals that the advantage of FIBIR compared to the CT-equivalent increases with increasing dose. In fact, we see that for higher doses the grey line diverges from the blue lines. We find the intersection points of the three FIBIR lines with the iterative CT equivalent line at (289$\mathrm {\mu }$m, 6.3mGy) for a visibility of 9.4$\%$, at (314$\mathrm {\mu }$m, 2.5mGy) for a visibility of 13$\%$ and at (338$\mathrm {\mu }$m, 1.1mGy) for a visibility of 17.8$\%$. The intersection point can be regarded as the point beyond which GI-CT yields superior image quality compared to the CT-equivalent. FIBIR thus allows to either increase the spatial resolution for a given dose or to reduce the dose for a given spatial resolution. For example, if we wish to achieve a resolution of 250$\mathrm {\mu }$m, GI allows to reduce the dose from 31.5mGy to 16.2mGy with a visibility of 9.4$\%$, from 31.4mGy to 8.9mGy with a visibility of 13$\%$, and from 31.5mGy to 5.2mGy with a visibility of 17.8$\%$.
Fig. 11. Dose vs. spatial resolution at a CNR of 5. On the left: results on simulated data, on the right: results on real data. The CT-equivalent FBP reconstruction is plotted in black, the iterative CT-equivalent in grey. The FIBIR results with three different visibilities are plotted in blue, the results of the post-processing fusion are plotted in orange. Dashed lines indicate results which are based on phase stepping, solid lines indicate results coming from static scans. The shade of the lines varies with the visibility, with the brightest shades corresponding to 9.4$\%$, the intermediate one to 13$\%$ and the darkest one to 17.8$\%$ visibility. The solid horizontal line indicates the upper dose limit of clinical breast CT [9].
When comparing the results of the static scan with the results of the phase-stepping scan we see that the latter leads to slightly inferior results compared to the static scans across all visibilities. Sharper gradients and access to more measurement points can thus not compensate for the significantly higher amount of noise in the data and for the associated statistical phase wrapping [5]. The 9.4$\%$ visibility line now crosses with the iterative CT-equivalent line at (256$\mathrm {\mu }$m, 23.7mGy), the 13$\%$ visibility line at (311$\mathrm {\mu }$m, 2.85mGy) and the 17.8$\%$ line at (352$\mathrm {\mu }$m, 0.7mGy). With phase-stepping data the dose required to achieve 250$\mathrm {\mu }$m resolution can be reduced from 31.5mGy to 11.7mGy and to 6.1mGy, with visibilities of 13$\%$ and 17.8$\%$, respectively.
When phase-stepping is performed it is also possible to retrieve the signals in the sinogram space, then analytically reconstruct them and finally fuse the two channels in a post-processing step in which the absorption channel is low-pass filtered and the phase channel high-pass filtered [13,39]. The results of this approach are plotted in different shades of orange for the three different visibilities in Fig. 11. We can see that this post-processing fusion approach follows the same trend as the FIBIR reconstruction but it yields inferior results across all dose levels and visibilities. The results of the post-processing fusion are less favourable than in [13] since in our work we estimated the spatial resolution with Fourier ring correlation, whereas in [13] the results show the dose as a function of the Gaussian kernel size needed to achieve a CNR of 5, which constitutes a lower limit on the resolution [13]. The estimation of the resolution in this work does therefore not constitute a measure of the resolving power of our scanner. On the contrary, since it is based on Fourier ring correlation, it is a task-based assessment, i.e., it estimates the resolution based on the contrast in a given image. Importantly, since image resolution is a function of the contrast, samples with different contrast could yield superior resolution estimates.
The results on real data in Fig. 11 show a similar trend as the results of the simulation study. At low doses, conventional CT is superior to GI-CT, while from 19.4mGy upwards, GI-CT delivers higher spatial resolution at a CNR of 5. This data was obtained with the current visibility of 9.4$\%$ of our setup. The highest dose reconstruction, at 66mGy, has been used as a reference for the FRC-based resolution analysis. The post-processing fusion approach yields inferior results compared to both the FIBIR approach and to the CT-equivalent. Note that the dose at the intersection point between the CT-equivalent line and the FIBIR line is comparable to the dose of the intersection point in the phase-stepping simulation study (23.7mGy). Since the visibility of our setup is expected to increase in the coming years owing to improvements in grating fabrication methodologies, the dose beyond which GI-CT beats conventional CT will likely be lowered. Notably, while the trends in the simulation study and on real data are very similar, the achievable resolution values are significantly reduced for the real data. This is to be expected since the in-silico study simulates an ideal scenario with a monochromatic beam, an infinitely small source size and no scattering.
Finally, a look at Fig. 12, which plots the reconstruction results at 22mGy (i.e., just beyond the crossing point), confirms that even with a low visibility of 9.4$\%$ GI-CT in conjunction with the proposed FIBIR algorithm is able to yield visually clearly sharper images (see red arrows in Fig. 12).
Fig. 12. Reconstruction results of a fixed breast specimen measured with phase-stepping at our GI scanner. On the left, from top to bottom: analytically reconstructed CT-equivalent scan and iteratively reconstructed CT-equivalent scan. On the right, from top to bottom: post-processing fusion and FIBIR. Stripe- [37] and ring-removal [38] have been applied to all data sets in the sinogram space and in the image space, respectively. This filtering step did not degrade image quality.
4. Discussion
In this work we have introduced a simple parameterization which allows to iteratively reconstruct a fused image that automatically combines the strengths of absorption- and phase-contrast images in GI-BCT. We have shown that the proposed fused intensity-based iterative reconstruction algorithm (FIBIR) allows grating-interferometry CT to outperform conventional CT in a clinical dose regime below 26.1mGy [9]. In particular, we were able to show that GI-CT could enable dose reduction or enhanced spatial resolution of breast CT examinations. The degree of this improvement strongly depends on the visibility of the GI setup. At very low doses conventional CT is superior to GI-CT, which is attributable to the phase channel’s noise sensitivity since, at low counts, statistical phase wrapping [5] severely degrades the quality of the signal.
The function $g_\theta$ relating the attenuation and phase signals is the key concept of our approach. In this article we have assumed a simple piece-wise linear relation between the absorption- and phase-contrast channel, but the model can be extended to include more tissue types and therefore include additional, more precise a-priori knowledge.
Limitations of this simple $g_\theta$ relation are that the mapping only takes into account grey levels, as well as that it is surjective and not bijective. In other words, different $\mu$ values can be mapped to the same $\delta$ value, whereas the opposite is not true. This limitation could be lifted by using a more powerful neural network-based mapping which takes into account not only the grey levels, but also the spatial context by using neighborhood information. Care must be taken if a more complex mapping is used since the loss with respect to the input channels must remain sufficiently smooth to allow for efficient gradient-based optimization. Therefore, input-convex neural networks [40] could be a suitable candidate.
There are three critical parameters that are necessary to make our approach superior to classical CT: the photon count, the GI setups’s visibility and DPC sensitivity. The photon count is limited by the maximum tolerable dose, whereas the other two depend on the system geometry. A larger Talbot order allows for achieving higher sensitivity but at the same time the spectral acceptance of the gratings decreases, which in turn causes the visibility of the interferometer to decrease. Larger Talbot orders are thus expected to yield FIBIR gradients which are sharper (because of the higher DPC sensitivity) but noisier (because of the lower visibility), whereas shorter setups are expected to yield blurrier but less noisy gradients. Future work will have to investigate which trade-off works best in conjunction with the proposed FIBIR algorithm.
We have shown that our approach is compatible with both phase-stepping and single-shot acquisitions. It can therefore readily be implemented for more complex acquisitions such as on a continuously rotating gantry in a helical geometry.
Future work will be done to improve the forward model of the FIBIR algorithm to more realistically simulate the physics of our system and further improve image quality of real samples.
We have seen that the FIBIR algorithm, while having to deal with nosier input data (since it operates at half the counts), it introduces more signal at each iteration compared to classical CT. Since noise can be dealt with by using an appropriate prior whereas a lack of signal cannot be handled by regularization, we expect a regularizer to increase the improvement in FIBIR compared to classical CT even further. Likewise, we expect phase-stepping scans to benefit more from regularization than single-shot data since in the former more signal and more noise are added per iteration.
In conclusion, FIBIR could become an essential tool in enabling GI-BCT to transition into the clinics. However, several challenges must first be overcome. Improvements in grating fabrication will increase the visibility of the setup and will thus reduce the dose at which GI-BCT outperforms conventional CT. Furthermore, we are currently extending our FIBIR algorithm to better deal with mechanical vibrations that lead to temporal instabilities in the interference pattern. Lastly, at present we are translating the GI-BCT technology from a static setup to a continuously rotating gantry, which is essential for clinical compatibility.
Funding
Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung (CRSII5 183568, SNF Sinergia Grant); ETH-Research Commission (ETH-12 20-2); ETH Doc.Mobility (Fellowship); Promedica Stiftung Chur; Swisslos Lottery Fund of Kanton Aargau.
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
References
1. A. Momose, S. Kawamoto, I. Koyama, Y. Hamaishi, K. Takai, and Y. Suzuki, “Demonstration of x-ray Talbot interferometry,” Jpn. J. Appl. Phys. 42(Part 1, No. 1), 1–15 (2003). [CrossRef]
2. S. A. Zhou and A. Brahme, “Development of phase-contrast X-ray imaging techniques and potential medical applications,” Phys. Medica 24(3), 129–148 (2008). [CrossRef]
3. A. Momose, “X-ray phase imaging reaching clinical uses,” Phys. Medica 79, 93–102 (2020). [CrossRef]
4. W. Yashiro, Y. Terui, K. Kawabata, and A. Momose, “On the origin of visibility contrast in x-ray talbot interferometry,” Opt. Express 18(16), 16890–16901 (2010). [CrossRef]
5. R. Raupach and T. G. Flohr, “Analytical evaluation of the signal and noise propagation in x-ray differential phase-contrast computed tomography,” Phys. Med. Biol. 56(7), 2219–2244 (2011). [CrossRef]
6. Z. Wang, N. Hauser, G. Singer, M. Trippel, R. A. Kubik-Huch, C. W. Schneider, and M. Stampanoni, “Non-invasive classification of microcalcifications with phase-contrast X-ray mammography,” Nat. Commun. 5(1), 3797 (2014). [CrossRef]
7. C. J. D’Orsi, E. A. Sickles, and E. B. Mendelson, ACR BI-RADS Atlas, Breast Imaging Reporting and Data System (American College of Radiology, 2013).
8. P. Skaane, S. Sebuødegård, A. I. Bandos, D. Gur, B. H. Østerås, R. Gullien, and S. Hofvind, “Performance of breast cancer screening using digital breast tomosynthesis: results from the prospective population-based oslo tomosynthesis screening trial,” Breast Cancer Res. Treat. 169(3), 489–496 (2018). [CrossRef]
9. Y. Zhu, A. M. O’Connell, Y. Ma, A. Liu, H. Li, Y. Zhang, X. Zhang, and Z. Ye, “Dedicated breast CT: State of the art—Part ii. Clinical application and future outlook,” Eur. Radiol. 32, 2286–2300 (2022). [CrossRef]
10. M. Stampanoni, Z. Wang, T. Thüring, C. David, E. Roessl, M. Trippel, R. A. Kubik-Huch, G. Singer, M. K. Hohl, and N. Hauser, “The first analysis and clinical evaluation of native breast tissue using differential phase-contrast mammography,” Invest. Radiol. 46(12), 801–806 (2011). [CrossRef]
11. C. Arboleda, Z. Wang, and M. Stampanoni, “Tilted-grating approach for scanning-mode X-ray phase contrast imaging,” Opt. Express 22(13), 15447 (2014). [CrossRef]
12. C. Arboleda, V. Lutz-Bueno, Z. Wang, P. Villanueva-Perez, M. Guizar-Sicairos, M. Liebi, Z. Varga, and M. Stampanoni, “Assessing lesion malignancy by scanning small-angle x-ray scattering of breast tissue with microcalcifications,” Phys. Med. Biol. 64(15), 155010 (2019). [CrossRef]
13. M. Rawlik, A. Pereira, S. Spindler, et al., “Refraction beats attenuation in breast CT,” arXiv, arXiv:2301.00455 (2023). [CrossRef]
14. M. Engelhardt, J. Baumann, M. Schuster, C. Kottler, F. Pfeiffer, O. Bunk, and C. David, “High-resolution differential phase contrast imaging using a magnifying projection geometry with a microfocus x-ray source,” Appl. Phys. Lett. 90(22), 224101 (2007). [CrossRef]
15. C. Kottler, F. Pfeiffer, O. Bunk, C. Grünzweig, and C. David, “Grating interferometer based scanning setup for hard x-ray phase contrast imaging,” Rev. Sci. Instrum. 78(4), 043710 (2007). [CrossRef]
16. H. F. Talbot, “LXXVI. Facts relating to optical science. No. IV,” The London, Edinburgh, Dublin Philos. Mag. J. Sci. 9(56), 401–407 (1836). [CrossRef]
17. M. Bech, O. Bunk, T. Donath, R. Feidenhans’l, C. David, and F. Pfeiffer, “Quantitative x-ray dark-field computed tomography,” Phys. Med. Biol. 55(18), 5529–5539 (2010). [CrossRef]
18. S. van Gogh, S. Mukherjee, J. Xu, Z. Wang, M. Rawlik, Z. Varga, R. Alaifari, C.-B. Schönlieb, and M. Stampanoni, “Iterative phase contrast ct reconstruction with novel tomographic operator and data-driven prior,” PLoS One 17(9), e0272963 (2022). [CrossRef]
19. A. Momose, S. Kawamoto, I. Koyama, Y. Hamaishi, K. Takai, and Y. Suzuki, “Demonstration of x-ray Talbot interferometry,” Jpn. J. Appl. Phys. 42(Part 2, No. 7B), L866–L868 (2003). [CrossRef]
20. F. Pfeiffer, T. Weitkamp, O. Bunk, and C. David, “Phase retrieval and differential phase-contrast imaging with low-brilliance X-ray sources,” Nat. Phys. 2(4), 258–261 (2006). [CrossRef]
21. M. V. Teuffenbach, T. Koehler, A. Fehringer, M. Viermetz, B. Brendel, J. Herzen, R. Proksa, E. J. Rummeny, F. Pfeiffer, and P. B. Noël, “Grating-based phase-contrast and dark-field computed tomography: A single-shot method,” Sci. Rep. 7(1), 7476 (2017). [CrossRef]
22. B. Brendel, M. Von Teuffenbach, P. B. Noël, F. Pfeiffer, and T. Koehler, “Penalized maximum likelihood reconstruction for x-ray differential phase-contrast tomography,” Med. Phys. 43(1), 188–194 (2015). [CrossRef]
23. J. Xu, Z. Wang, S. Van Gogh, M. Rawlik, S. Spindler, and M. Stampanoni, “Intensity-based iterative reconstruction for helical grating interferometry breast CT with static grating configuration,” Opt. Express 30(8), 13847–13863 (2022). [CrossRef]
24. P. Zhu, K. Zhang, Z. Wang, Y. Liu, X. Liu, Z. Wu, S. A. McDonald, F. Marone, and M. Stampanoni, “Low-dose, simple, and fast grating-based X-ray phase-contrast imaging,” Proc. Natl. Acad. Sci. 107(31), 13576–13581 (2010). [CrossRef]
25. Z. Wu, K. Gao, Z. Wang, S. Wang, P. Zhu, Y. Ren, and Y. Tian, “Generalized reverse projection method for grating-based phase tomography,” J. Synchrotron Radiat. 28(3), 854–863 (2021). [CrossRef]
26. “DABAX library,” http://ftp.esrf.eu/pub/scisoft/xop2.3/DabaxFiles.
27. M. Berger, J. Hubbell, S. Seltzer, J. Chang, J. Coursey, R. Sukumar, D. Zucker, and K. Olsen, “NIST Standard Reference Database 8 (XGAM),” NIST, PML, Radiation Physics Division (2010).
28. D. R. White, R. V. Griffith, and I. J. Wilson, “ICRU Report 46,” J. ICRU Meas. os-24, 203–205 (1992). [CrossRef]
29. M. Willner, M. Bech, J. Herzen, I. Zanette, D. Hahn, J. Kenntner, J. Mohr, A. Rack, T. Weitkamp, and F. Pfeiffer, “Quantitative X-ray phase-contrast computed tomography at 82 keV,” Opt. Express 21(4), 4155 (2013). [CrossRef]
30. J. Nocedal, “Updating Quasi-Newton Matrices with Limited Storage,” Math. Comp. 35(151), 773–782 (1980). [CrossRef]
31. M. Abadi, A. Agarwal, P. Barham, et al., “TensorFlow: Large-Scale Machine Learning on Heterogeneous Distributed Systems,” software available from tensorflow.org (2015).
32. NVIDIA, P. VingelmannF. H. Fitzek, “Cuda, release: 10.2.89,” (2020).
33. S. van Gogh, Z. Wang, M. Rawlik, C. Etmann, S. Mukherjee, C. Schönlieb, F. Angst, A. Boss, and M. Stampanoni, “INSIDEnet: Interpretable nonexpansive data-efficient network for denoising in grating interferometry breast CT,” Medical Physics 49(6), 3729–3748 (2022). [CrossRef]
34. W. van Aarle, W. J. Palenstijn, J. De Beenhouwer, T. Altantzis, S. Bals, K. J. Batenburg, and J. Sijbers, “The ASTRA Toolbox: A platform for advanced algorithm development in electron tomography,” Ultramicroscopy 157, 35–47 (2015). [CrossRef]
35. A. Rose, Vision: Human and Electronic (Springer, 1973).
36. R. P. J. Nieuwenhuizen, K. A. Lidke, M. Bates, D. L. Puig, D. Grünwald, S. Stallinga, and B. Rieger, “Measuring image resolution in optical nanoscopy,” Nat. Methods 10(6), 557–562 (2013). [CrossRef]
37. B. Münch, P. Trtik, F. Marone, and M. Stampanoni, “Stripe and ring artifact removal with combined wavelet—fourier filtering,” Opt. Express 17(10), 8567–8591 (2009). [CrossRef]
38. D. Gürsoy, F. De Carlo, X. Xiao, and C. Jacobsen, “Tomopy: a framework for the analysis of synchrotron tomographic data,” J. Synchrotron Radiat. 21(5), 1188–1193 (2014). [CrossRef]
39. Z. Wang, C. A. Clavijo, E. Roessl, U. Van Stevendaal, T. Koehler, N. Hauser, and M. Stampanoni, “Image fusion scheme for differential phase contrast mammography,” J. Instrum. 8(07), C07011 (2013). [CrossRef]
40. B. Amos, L. Xu, and J. Z. Kolter, “Input convex neural networks,” arXiv, arXiv:1609.07152 (2016). [CrossRef]
