I.

INTRODUCTION

MALIGNANCIES of the breast still represent the most prevalent cancer in women [1]. Unfortunately, none of the currently used breast imaging techniques (mammography, breast ultrasound, breast MRI and absorption-based breast CT and tomosynthesis [2], [3]) is able to provide fully three-dimensional images with sufficiently high isotropic resolution and soft-tissue contrast necessary to identify critical breast cancer imaging biomarkers [4]. Therefore, better imaging modalities are needed to improve early detection and increase survival rates. X-ray phase-contrast CT could potentially offer a solution by combining the high three-dimensional resolution, which comes with CT, with superior soft-tissue contrast.

When X-ray waves interact with matter, their amplitude and phase are modified according to the refractive index of the material they interact with. The refractive index of a material is given by n =1 − δ + iβ. The real part δ dictates the change in the beam’s phase Φ as

From this, the refraction angle α can be computed using

The imaginary part β is directly linked to the attenuation coefficient via μ = 4πβ/λ, which can then be used to compute the beam’s attenuation by using the Beer-Lambert law.

It is widely known that soft tissues are characterized by similar β’s [5], which makes it difficult to distinguish different tissue types in conventional absorption-based CT. Conversely, larger differences in δ’s [5] can theoretically yield higher soft tissue contrast in reconstructed phase-contrast CT volumes.

The X-ray’s phase must be computed indirectly and many approaches to achieve this have been proposed over the years. Grating interferometry [6], [7], [8] arguably has the highest chance of making the transition to clinical practice. In fact, it has non-restrictive requirements in terms of temporal and spatial coherence of the X-ray beam, it can be operated at large fields-of-view (FOV) and it has a comparably high mechanical robustness [6].

Grating interferometry detects the X-ray’s refraction angle α by exploiting a peculiar interference pattern called Talbot carpet [9]. When an X-ray beam is refracted, this results in a lateral shift in the interference pattern. Therefore, by measuring this shift, the wavefront’s change in phase can be easily obtained by integrating (2). To obtain this Talbot carpet, three gratings are positioned between the source and the detector [7]. The first grating (source grating or G0) is composed of a highly absorbing material such as gold and is placed immediately in front of the X-ray tube to improve beam coherence. The second grating (phase grating or G1), which is not designed to absorb photons, imposes a significant phase-shift to the X-ray beam and creates the interference pattern. To measure the lateral shift of this interference pattern induced by the sample, which is in the μm range, a highly resolving detector would be required. Unfortunately, to date no such detectors exist. Therefore, to circumvent this problem, a highly absorbing third grating (analyzer grating or G2) is placed in front of the detector. By moving one of the gratings with respect to the others in x-direction, it is then possible to obtain an interferogram called phase stepping curve [10], from which it is in turn possible to compute the lateral shift of the interference pattern.

The interferogram is modeled as

where I₀, V₀, and Φ₀ are the flat-field intensity, visibility and phase maps, respectively, and k is the k-th phase step. The transmission T, the dark-field D and the differential phase sinograms are given by:

with λ being the wavelength, g₂ the pitch of the G2 grating and d₂ the distance between the origin and G2.

The interference pattern’s shift φ, which is directly linked to the beam’s refraction and thus phase, can then be retrieved with Fourier analysis. The same holds for the absorption signal, which is related to the average intensity of the curve, and for the dark-field signal, which is related to the curve’s amplitude. In this work though we will focus exclusively on phase. When combined with a CT acquisition protocol, grating interferometry (GI) naturally extends to GI-CT.

In an attempt to bring GI-CT to clinical practice, our group has embarked in a long term effort to build a first-of-its-kind Grating Interferometry Breast CT (GI-BCT) prototype. As such, deposited radiation dose, scanning times and patient comfort must be compatible with clinical standards.

With these constraints in place, high quality images with grating interferometry CT represent a challenging objective. In fact, to date, the successful use of grating interferometry phase-contrast CT has been limited to synchrotron beamlines [11] and laboratory setups [12] where high image quality is achieved by a high X-ray flux or by long scanning times.

There are two main reasons why it is so challenging to obtain high-quality images with GI-CT. First, an intrinsic noise amplification takes place during signal retrieval [13]. Second, the differential nature of the phase-contrast forward operator causes the inverse problem to be more ill-conditioned as compared to conventional CT.

The main figure of merit which determines the quality of a grating interferometer is its visibility, i.e. the amplitude of the interference pattern in a flat-field scan. The higher the visibility, the more precisely one can compute the interference pattern’s lateral shift and, consequently, the X-ray beam’s refraction. Visibilities of 30% have been reported for polychromatic 46keV setups [14]. Unfortunately, high visibilities are challenging to achieve for systems which are shorter and/or have higher sensitivity, which requires the grating structures to be smaller.

Given the highly noisy sinogram data, it is important to reconstruct the tomograms with a stable inversion algorithm. A pseudo-inverse such as the filtered backprojection (FBP) algorithm could be applied in conjunction with the Hilbert filter [15] to solve this task. However, it is widely accepted in the CT community that iterative reconstruction algorithms are better suited to deal with highly ill-conditioned problems.

Typically, in iterative reconstruction, we need to define a variational loss function comprising a data-fidelity term and a regularization functional which incorporates prior knowledge about the expected reconstruction. Minimization of this loss with an optimizer of choice then allows to reconstruct an image which is simultaneously consistent with the measurements as well as with the prior knowledge.

The most widely used prior in image reconstruction is the total variation (TV) prior [16], which promotes homogeneous regions separated by sharp edges in the solution by assuming piece-wise constant signals. While TV is still considered to be a powerful baseline algorithm to regularize ill-conditioned inverse problems, recent years have witnessed the rise of data-driven algorithms which outperform traditional methods and thus constitute the current state-of-the-art in the field.

Two main data-driven approaches that draw inspiration from classical variational optimization schemes have been proposed in the literature. The first approach comprises end-to-end methods which unroll iterative schemes, thereby transforming each iteration of the iterative reconstruction algorithm into a distinct layer of a neural network [17], [18]. Since the imaging physics is embedded into the network, these models are generally believed to be more robust to noise and adversarial perturbations compared to pure black-box neural networks. In the second approach, the idea is to learn the regularizer a-priori on a representative training set, and to then use the trained regularizer in conjunction with the data-fidelity term in a classical variational optimization framework [19], [20], [21], [22], [23].

The first type of approach tends to yield superior results and deliver faster reconstructions [18]. However, it has the important disadvantage that no convergence guarantees can be derived and that it needs explicit supervision. Algorithms in the second category tend to be slower and potentially yield slightly inferior results. However, they have several advantages. The algorithms can be trained independently of the forward operator, they are more data-efficient (and possibly unsupervised), and they are amenable to stability and convergence analysis [21]. Since it is of utmost importance in the medical field to reliably reconstruct the tomograms, we propose a novel algorithm that fits into the second category.

Many different ideas have been proposed to learn a regularizer in a data-driven manner. Some notable ones are adversarial regularization [20] and its convex counterpart [21], score matching networks [24], regularization by denoising (RED) [19], network Tikhonov (NETT) [25] and data-driven Plug-and-Play denoisers [22], [23].

To address some of the challenges in GI-BCT and to be able to reconstruct phase-contrast tomograms from highly corrupted measurement data, we propose an iterative reconstruction algorithm which leverages the power of deep learning to regularize the highly ill-conditioned tomographic inversion problem. In particular, we propose an algorithm that alternates between data updates governed by the L-BFGS algorithm [26] and regularization steps performed with a denoising deep neural network in a Plug-and-Play fashion [27]. We apply the proposed approach to both simulated data and real measurements and show that it achieves excellent results.

II.

METHODS

The optimization problem we aim to solve is

where δ is the image containing the real part of the index of refraction, φ is the retrieved differential phase-contrast sinogram, and A is a linear forward operator modeling equation (6).

Unfortunately, (7) is a highly ill-conditioned problem. Since an unregularized optimization of (7) will converge to a highly unstable solution, a powerful regularization strategy is necessary to stabilize the reconstruction.

As mentioned in the introduction, data-driven regularizers have emerged as the new state-of-the-art in the field and now routinely outperform classical regularization strategies such as total variation (TV) [16]. Among the many proposed strategies, we found data-driven Plug-and-Play regularizers to work best. Therefore, we here propose to learn a regularization network which is able to remove both noise and artefacts from the image iterates as they converge to the final reconstruction.

Given corrupted and clean images, we wish to train a network f_θ which maps the corrupted image δ to its clean counterpart δ*. The objective function to achieve this is

where f_θ is an deep neural network with sufficient expressive power and {δ_i, } are training data pairs.

We parameterized f_θ with a 7 million parameter, bias-free U-net [28]. We removed the biases because it has been shown that this leads to 1) more interpretable denoising performance, and 2) that biasless networks generalize better to different noise amplitudes [29]. The higher interpretability comes from the fact that we can regard the denoising process as being locally linear [29]. Therefore, computing the Jacobian shows how the pixel neighborhood is used for denoising a particular pixel. The higher robustness is important in our case since 1) the noise encountered during image reconstruction might slightly differ from iteration to iteration, and 2) it makes it easier to train the network as we do not have to find the perfect amount of noise to train on.

To solve (7) we use the L-BFGS optimization scheme proposed by [26] and apply a denoising step after every k-th iteration. The proposed method is summarized in Algorithm 1.

Algorithm 1:

L-BFGS optimization with data-driven Plug-and-Play denoiser

By removing noise and artefacts, the denoising step can be interpreted as a projection of the current image iterate to the data manifold to which clean reconstructions belong. Iteratively alternating between data-fidelity optimization and projection to the data manifold thereby allows to get close to the measured data, while staying close to the data manifold of clean reconstructions.

III.

RESULTS

A.

Simulated data

To quantitatively assess the performance of our method we applied it to simulated breast phantoms developed in-house. We simulated 600 projections and added Poisson noise that matched real measurements. The results in Figure 1 show that the proposed method achieves excellent results and clearly outperforms analytical reconstruction, both qualitatively as well as in terms of structural similarity index (SSIM) and mean squared error (MSE).

Fig. 1.

Reconstruction results on simulated breast phantoms. First column: full slice, second column: the part inside the red region enlarged.

B.

Real data

To demonstrate the effectivenss of our method on real data, we scanned a fixed mastectomy obtained at the University Hospital Zürich on our GI-BCT prototype. We acquired 600 projections under continuous circular rotation. A total of 10 scans have been averaged to compensate for the low visibility (thus high noise amplitudes) we are currently working with.

We applied exactly the same algorithm as for the simulated data. Importantly, we used the neural network that has been trained on simulated data. Figure 2 shows that the proposed method once again achieves excellent results and clearly outperforms analytical reconstruction. We want to emphasize that the applicability of a denoiser, that was trained on simulated data, to real-world data, is a crucial advantage of our method since it is very difficult to obtain high-quality real data in a medical setting.

Fig. 2.

Reconstruction results on mastectomy data. First column: full slice, second column: the part inside the red region enlarged.

IV.

CONCLUSION

In this article we have proposed a novel iterative reconstruction algorithm with a data-driven denoising prior and have shown that it is able to produce excellent results. Importantly, the regularizer can be trained on simulated data and later be applied to real measurements. Ongoing work focuses on the theoretical convergence guarantees of the proposed algorithm.