INSIDEnet: Interpretable NonexpanSIve Data-Efficient network for denoising in grating interferometry breast CT.

PURPOSE: Breast cancer is the most common malignancy in women. Unfortunately, current breast imaging techniques all suffer from certain limitations: they are either not fully three dimensional, have an insufficient resolution or low soft-tissue contrast. Grating interferometry breast computed tomography (GI-BCT) is a promising X-ray phase contrast modality that could overcome these limitations by offering high soft-tissue contrast and excellent three-dimensional resolution. To enable the transition of this technology to clinical practice, dedicated data-processing algorithms must be developed in order to effectively retrieve the signals of interest from the measured raw data.
METHODS: This article proposes a novel denoising algorithm that can cope with the high-noise amplitudes and heteroscedasticity which arise in GI-BCT when operated in a low-dose regime to effectively regularize the ill-conditioned GI-BCT inverse problem. We present a data-driven algorithm called INSIDEnet, which combines different ideas such as multiscale image processing, transform-domain filtering, transform learning, and explicit orthogonality to build an Interpretable NonexpanSIve Data-Efficient network (INSIDEnet).
RESULTS: We apply the method to simulated breast phantom datasets and to real data acquired on a GI-BCT prototype and show that the proposed algorithm outperforms traditional state-of-the-art filters and is competitive with deep neural networks. The strong inductive bias given by the proposed model's architecture allows to reliably train the algorithm with very limited data while providing high model interpretability, thus offering a great advantage over classical convolutional neural networks (CNNs).
CONCLUSIONS: The proposed INSIDEnet is highly data-efficient, interpretable, and outperforms state-of-the-art CNNs when trained on very limited training data. We expect the proposed method to become an important tool as part of a dedicated plug-and-play GI-BCT reconstruction framework, needed to translate this promising technology to the clinics.

INTRODUCTION

Breast cancer is the most prevalent malignancy in women with one out of eight developing the disease in her lifetime. To fight this public health burden, recent years have witnessed the introduction of many screening programs to early detect and consequently to better treat this disease. Unfortunately, both false positive and false negative rates remain high, leading to unnecessary psychological distress and missed tumors, respectively. The reason for this is attributable to the fact that current breast imaging techniques, most notably mammography, ultrasound, tomosynthesis, breast MRI, and absorption‐based breast CT, , all suffer from some limitations. In fact, none of these techniques simultaneously yields fully three‐dimensional (3D) data with sufficient soft‐tissue contrast and spatial resolution to detect crucial imaging biomarkers (small soft tissue lesions and their margins, microcalcifications, architectural distortions, and tiny soft tissue density differences), thereby making it difficult for radiologists to take confident decisions.

Consequently, there has been an ever‐increasing effort to utilize X‐ray phase contrast imaging which can potentially lead to higher soft‐tissue contrast compared to absorption‐based imaging, without sacrificing spatial resolution. In fact, with synchrotron sources and in laboratory applications phase contrast already delivers far superior soft‐tissue delineation. , X‐ray grating interferometry (GI) is a phase contrast technique holding most of the prerequisites for clinical compatibility. For this reason, our group has designed and is currently building a grating interferometry breast computed tomography (GI‐BCT) prototype, that is, translating the technology into a first‐of‐its‐kind compact medical device.

Obtaining high‐quality phase‐contrast images in clinically compatible settings remains a challenge. The theoretically achievable higher contrast in phase imaging compared to absorption cannot yet be fully exploited on compact X‐ray sources because of the high‐noise amplitudes and noise distributions. Phase contrast images are in fact characterized by intrinsic low‐frequency noise, especially in the low photon‐count case. Furthermore, those very same images contain a nonuniform noise distribution due to imperfect gratings. Finally, adhering to the severe constraints imposed by the clinical environment (such as radiation dose, scanning time, and patient comfort) demands for novel solutions to handle sparse sampling and photon starvation, thereby making it even more cumbersome to achieve high image quality.

It is important to mention here that this problem is specific to the phase contrast channel. In fact, in the absorption image the noise is predominantly present at the high frequencies, making it much easier to effectively denoise these data and for which a variety of established denoising methods , work well.

A powerful denoising algorithm is thus necessary to suppress the noise and let the higher intrinsic phase contrast emerge. In particular, a pipeline is needed to cope with the high‐noise amplitudes and heteroscedasticity, while simultaneously offering high algorithm interpretability, reliability, and robustness.

The two major groups of denoising algorithms, which exist today, namely traditional methods and deep learning methods, have both significant limitations. Algorithms in the first category rely on hand‐crafted priors , such as nonlocal similarity, low rankness, small gradient norms, nonnegative values, and sparsity in some transform domains, such as the wavelet domain. They form reliable pipelines which work extremely well on images with relatively little noise, but their performance drastically deteriorates when dealing with lower image quality, as they are unable to adapt to the data at hand. In contrast, deep learning methods have been shown to yield impressive results , , by implicitly learning a prior in a data‐driven fashion without the need of human intervention. A big downside of deep learning–based algorithms, however, is their high sensitivity to training data and their limited interpretability caused by the concatenation of (many) linear and nonlinear operations. It is well known that deep network–based denoisers might add structures which are not present and remove the ones which are present in the ground‐truth data. Such behavior is unacceptable in particular when it comes to denoising in medical applications. A third alternative is to combine the structure of handcrafted regularization with data‐driven learning, thereby leveraging the power of data while maintaining interpretability.

To date, to the best of our knowledge, a single article has been published on the use of a data‐driven method for denoising of GI, and in particular differential phase contrast (DPC), projection data. However, the authors used a black‐box model and restricted their analysis to radiography. In contrast, we focus on GI phase CT and propose a hybrid denoising algorithm, which we call “Interpretable NonexpanSIve Data‐Efficient network” (INSIDEnet), that attempts to leverage the strengths of both worlds: the interpretability of classical filters and the flexibility of data‐driven models. Importantly, the model has been parameterized to maximize interpretability and reliability, which are both imperative conditions for clinical applicability.

In this article, we demonstrate the performance of our hybrid algorithm on simulated breast phantoms and real data acquired on our GI‐BCT prototype. Our INSIDEnet achieves better results compared to traditional, non‐learning‐based models, without paying in robustness. The INSIDEnet shows that it is possible to achieve an excellent denoising performance with superior interpretability and robustness compared to deep neural networks.

MATERIALS AND METHODS

Grating interferometry breast CT

Conventional X‐ray imaging is based on photon absorption in the imaged tissues. Unfortunately, biological soft tissues have very similar attenuation coefficients as they are all mainly composed of carbon, oxygen, and water. Therefore, there is limited contrast between different body constituents. For phase, in the absence of noise, the theoretically achievable contrast is higher because the real part of the index of refraction (related to phase) is orders of magnitude larger than the imaginary part (related to absorption). It still has to be explored though if such superior contrast can also be achieved in a clinically compatible setting. Contrary to attenuation‐based imaging, which directly measures the intensity of the transmitted X‐ray beam, it is not possible to directly detect the phase in a polychromatic setting, which therefore has to be measured indirectly. Several techniques have been proposed in this regard with the most notable ones being propagation based, crystal interferometry, analyzer based, edge‐illumination, and GI. , , While all methods can be applied at synchrotron light sources, GI has received special attention as it satisfies the prerequisites for clinical applicability: it has high mechanical robustness, can be scaled up to large fields of view (FOV) and only requires moderate spatial coherence and monochromaticity.

GI encodes propagation‐induced phase changes in the beam wavefront—when passing through a specimen—into an intensity modulation measured by a detector placed downstream. Its simplest configuration, Talbot interferometry, consists of two gratings placed in a partially coherent beam. The latter is usually provided by a third/fourth‐generation synchrotron source or, with significantly less intensity, by a microfocus X‐ray tube. The first grating is usually a phase grating, that is, it does not absorb the beam but imposes a phase‐shift resulting in a controlled wavefront modulation at a specific distance downstream, usually where the second, absorbing analyzer grating is placed. When the source does not provide a sufficiently high spatial coherence, like in the case of a conventional X‐ray tube, a third grating can be introduced right after the source yielding to the so‐called Talbot–Lau configuration, as shown in Figure 1.

FIGURE 1

Left: Schematic drawing of GI setup in Talbot–Lau configuration. Right: phase stepping curves, with (orange) and without a sample (blue). The logarithm of the ratio of the average of the blue and orange curve gives the absorption signal, the reduction in its amplitude gives the dark‐field signal, whereas the relative shift of the curve with respect to its reference gives the differential phase signal

The intensity modulation of the resulting fringe pattern is characterized by its visibility, while retrieval of the absorption, phase, and dark‐field signals can be done with various methods with phase stepping and fringe scanning , being the most commonly used approaches. By combining the interferograms obtained with and without sample, one can retrieve the absorption signal, the differential phase signal, which is proportional to the larger‐than‐pixel‐size refraction, and the dark‐field signal, which is proportional to the incoherent refraction on a scale smaller than the pixel size. Finally, by rotating the sample or the X‐ray source and detector, GI naturally extends to GI‐CT.

High‐quality phase‐contrast tomograms of breast tissue have been demonstrated with GI on small FOV laboratory setups and with propagation‐based, large FOV imaging on a synchrotron. In the first case, good image quality is obtained with long scanning times. The second case benefits from high sensitivity, thanks to a large propagation distance, and near‐perfect coherent source. A clinically compatible device must on the other hand handle sparse sampling and photon starvation. As known from absorption‐based X‐ray imaging, low photon counts lead to lower signal‐to‐noise ratio (SNR) and thus lower image quality. Unfortunately, in GI this effect is even more detrimental. First, the characteristic imaging scheme based on the acquisition of an interferogram in GI leads to an amplification of the counting noise, with each of the three channels having a different noise propagation. Second, noise amplification is not uniform across the FOV as, in fact, grating fabrication defects and grating misalignment cause strong local noise amplitude variations. Finally, during reconstruction, correlations between pixels are introduced and, crucially, integration of DPC leads to an amplification of low‐frequency noise, making the noise pattern even more difficult to deal with.

We would like to stress that since a clinically compatible one‐shot acquisition method does not allow to explicitly retrieve the sinograms for the three signals, the proposed INSIDEnet is ultimately envisioned to act as a denoising prior (or proximal operator) in a plug‐and‐play framework by iteratively denoising image iterates within a gradient‐based optimization scheme. This approach will be presented elsewhere and will not be further discussed in this article, as we would like to focus here on the introduction of the denoising engine itself.

Simulated breast phantoms

Clean breast phantoms

We generated 30 in silico breast phantoms of voxels with a voxel size of 100 (real phantom size of 0.44 15.36 15.36 cm) containing three main tissue types, namely adipose, glandular, and skin.

To start with, we randomly generated 10 000 ellipsoids of different sizes, shapes, positions, and orientations, followed by two thresholds. We then multiplied this preliminary phantom with the central part of two randomly rotated binary masks (one for the whole breast and one for the skin) obtained from a 3D mask of a real breast that has been acquired with a breast CT scanner at the University Hospital Zürich. These data were part of a retrospective analysis of patient data approved by the local ethics committee. All patients gave their written informed consent. The edges of the ellipsoids were used to simulate duct‐like structures radiating out of the glandular tissue.

Realistic attenuation coefficients and phase shift coefficients for the absorption and the phase image, respectively, were assigned to the different regions representing adipose tissue, glandular tissue, and skin. Phase shift values were calculated starting from decrements in the real part of the index of refraction and using the known relation , where is the X‐ray wavelength corresponding to the design energy of our prototype. and values have been calculated using NIST XCOM and f1f2 Kissel.dat of the DABAX library, respectively, based on the tissue definition in ICRU 46. Finally, we added anatomical noise to model more realistic slight and inhomogeneities in the tissues as well as breast lesions of different sizes and contrasts to investigate the algorithms' performance in terms of lesion detectability. One slice of a phantom pair for absorption and phase is shown in Figure 2. It is important to mention here that these values represent an ideal setting. In reality, it will not be possible to reconstruct these and values precisely due to X‐ray polychromaticity, Compton scattering, and limited phase sensitivity.

FIGURE 2

Absorption and phase contrast breast phantoms. The units for absorption are attenuation coefficients, while for phase it is phase shift coefficients

We would like to highlight that in this simulation, apart from differing contrasts between tissues, phase, and absorption images contain exactly the same structures, that is, the same information. While research into X‐ray phase contrast imaging is evidently being carried out in the hope that extra diagnostic information becomes available, to date we do not know how such extra information would look like at a macroscopic scale in mammary tissue at clinically compatible doses. We are therefore not able to realistically simulate such structures. Since absorption and phase contrast data will be processed independently in this paper, the lack of extra information in phase contrast images will not affect the validity of the proposed method.

Noisy breast phantoms

Each of the simulated breast phantom pairs was then used to generate noisy counterparts. First, we used the ASTRA toolbox projector to obtain differential phase and transmission sinograms as follows :

where is the sample‐G2 distance, is the pitch of G2. The sample‐G2 distance was 73 cm, the pitch of G2 was 4.2 m, and the system's design energy was 46 keV. The source‐to‐sample distance was 103 cm, and the sample‐to‐detector distance was 74 cm. More detailed information about our prototype is provided in Section 2.3.

Flat‐field data (intensity map , visibility map , and phase map ) were obtained from our scanner to provide simulations as realistic as possible. Due to grating imperfections, , were all highly inhomogeneous, as shown in Figure 3. Combining these data, we then simulated phase‐stepping curves with the sample in place:

FIGURE 3

Flat‐field data of our GI‐BCT prototype: intensity map , phase map , and visibility map

Here, is the th phase step, uniformly spaced between 0 and and is the intensity value measured at the th phase step. We neglected visibility reduction in the sample, that is, , as no significant small‐angle scattering is to be expected in breast tissue apart from microcalcifications, which, however, were not simulated here.

Likewise, the background phase stepping curves were simulated as follows:

We then simulated detector quantum noise by sampling from the Poisson distribution with mean . We empirically determined the necessary photons to match the image quality of our real data. We thus simulated 40 000 photons leaving the source at every exposure, which resulted in the flat‐field data displayed in the first row of Figure 3.

We would like to mention here that this corresponds to a much higher dose than what is allowed in clinical practice. Such a higher dose allows to compensate for the yet insufficient grating quality which severely impacts the visibility. Once grating quality, and especially flat‐field visibility, will improve, the same image quality will be achievable with significantly less dose.

By combining the two‐phase stepping curves and , we then retrieved both the differential phase and the attenuation signal with simple Fourier analysis. Finally, absorption and phase images were obtained by reconstructing the retrieved signals with analytical reconstruction algorithms available in the ASTRA toolbox. Importantly, the Hilbert filter had to be applied for reconstructing the phase contrast image.

All denoising algorithms were deployed on two‐dimensional (2D) slices of the reconstructed data. While the algorithm could be easily extended to work on 3D data, we did not pursue this because of GPU memory constraints.

While denoising could also have been performed in the projection domain, we found that due to strong local noise amplitudes in the measurements, this yielded far inferior results than denoising the reconstructions in which these strong local noise amplitudes have previously been attenuated by backprojection. Moreover, since the algorithm is envisioned to be used as a proximal operator in the image space, we wanted to test its ability to denoise the latter rather than the sinogram data.

10 volumes were used for training, 10 for validation during training, and 10 for testing. All reported metrics in the simulations study have been obtained on the 10 testing volumes.

Real data

To demonstrate the performance of our algorithm on real data we applied our algorithm trained on simulated data to a tomogram acquired on our GI‐BCT prototype of a chunk of meat with roughly the same volume as a human breast.

We used a fifth Talbot‐order symmetric interferometer with 4.2 m pitch, with a single G0, single G1, and three G2 gratings. All gratings are bent. As X‐ray source we used a Comet MXR‐225HP/11 tube operated at 70 kVp and 10 mA and as a detector a CdTe photon counting Dectris prototype with 75 m pixel size operated at 10 Hz. We acquired five scans in continuous rotation, each consisting of 600 projections and then averaged the tomograms to obtain the input for our denoising algorithm. This averaging of five scans has been performed to compensate for the (yet) insufficient flat field visibility of our scanner.

INSIDEnet: Interpretable NonexpanSIve Data‐Efficient network

The proposed method combines different image‐processing paradigms, namely transform domain thresholding, transform learning, multiscale processing, explicit orthogonality (thus nonexpansiveness), and deep learning.

With the assumption that the signal can be expressed as a linear combination of few basis elements, transform domain thresholding transforms image patches into another domain such as the wavelet domain or the discrete cosine transform (DCT) domain. Owing to the high levels of cross‐correlation within patches, this representation of the data will be highly sparse. Filtering can be performed by thresholding or shrinking the coefficients. Finally, the filtered patches are transformed back to the image domain and the patches are aggregated to form the final denoised image.

Data‐adaptive bases have been shown to be superior compared to handcrafted ones. Transform learning in particular builds upon the idea that it might be beneficial to replace a fixed hand‐designed operator such as the wavelet transform or the DCT with an operator adapted to the data at hand. It thus aims to learn a transform matrix which moves the images to a space in which they have a highly sparse representation. This transform operator can then be used to regularize ill‐posed inverse problems like denoising.

We propose to combine these two ideas with deep learning and multiscale image processing to efficiently denoise highly corrupted images. Crucially, an explicit orthogonality constraint has been used to achieve high model robustness and interpretability. To exclude channel cross‐talk and be able to inspect the quality of the two signals of interest independently, we applied the INSIDEnet to absorption and phase data separately. It should be noted though that the proposed method can be applied to multiple channels if correlations between the two signals should be leveraged.

The fixed transform is replaced with learnable matrices, and the images are processed across four different scales. This ensures that our model can remove the noise across a broad frequency band. The entire pipeline is end‐to‐end trainable, thus allowing to not only learn the transforms but also the thresholds, thereby leaving very few tunable hyperparameters.

In the simulation study, the models have been trained in a supervised manner to map noisy images to their clean counterparts. For quantitative and qualitative evaluation, the trained models have been applied to simulated testing data and real data, respectively.

We will first introduce the stacked orthogonal transform learning denoising pipeline which is applied across each image scale (see Figure 4), before explaining how the multiple scales are generated and combined (see Figure 5).

FIGURE 4

Denoising pipeline applied at each image scale. The trainable parameters are displayed in light gray. The learning signal enters the pipeline after the aggregation step. In yellow: noisy data, in blue: transformed data, and in green: denoised data

FIGURE 5

INSIDEnet overview. Image decomposition: the noisy input image gets first denoised at full resolution, with the aim to remove the highest frequencies. Next, downsampling leads to a lower resolution image with still all lower frequency noise present. Applying the denoising pipeline again removes the noise corresponding to high frequencies at this scale. This process is repeated three times, effectively removing noise at both high and low frequencies. Image reconstruction: high‐frequency components are obtained by subtracting the low frequencies at each scale. These denoised high‐frequency components are then iteratively added to the lowest frequency components in a Laplacian pyramid fashion

Image preprocessing

Before entering the denoising pipeline, the data are scaled to be within [0,1]. The same scaling is then applied to their noisy counterparts. During training, the images are randomly shuffled to ensure unbiased learning.

Stacked orthogonal transform learning

Let be the noisy input image and the corresponding target image.

is divided into overlapping patches (stride 2 was empirically determined to be optimal), here represented by tensor , where is the patch size in horizontal and vertical directions. Next, the patches are flattened, thereby effectively reshaping , which now is in . This tensor gets then multiplied through an Einstein sum with the orthogonal transform matrix yielding the transformed image patches

We enforce to be orthogonal by employing the Cayley transform. , An orthogonal matrix can in fact be obtained by computing , where is a skew‐symmetric matrix which can in turn be obtained by an arbitrary matrix :

We thus let be our trainable matrix, with random initialization, which is explicitly transformed into an orthogonal matrix.

The filtering step itself is very simple and draws inspiration from the proximal operator of the norm, that is, a hard threshold on the coefficient magnitudes. To allow this threshold to be trainable, we approximate the hard threshold with a steep sigmoid function: where is the of the sigmoid and is not trainable. It was selected by visually plotting the resulting activation function and exploiting the fact that is approximately in the same range (0–1) as . , which defines the threshold, was initialized to at all scales to ensure that at the start of the training all coefficients are kept, thus mapping the input image to itself. Element‐wise multiplication of with yields the filtered representation of our image:

The data are then transformed back to the image domain through an Einstein sum

The steps above are repeated times. This effectively coincides with iterative denoising since every iteration of transform domain thresholding improves the quality of the image by removing some noise. Alternatively, it can be interpreted as a layer of a neural network. However, where in convolutional neural networks (CNNs) each layer is just an unconstrained forward operator followed by an arbitrary nonlinearity, here each layer has a clear mathematical meaning and is stabilized by projecting back the filtered coefficients to the image space. Our network is in a way also more general than a convolutional network in the sense that the transform operators are not necessarily convolutional. This can compensate for a possible loss of transform expressiveness imposed by the orthogonality restriction.

Finally, the patches are rearranged back to their original position in the image to obtain .

We would like to stress that, while many algorithms have been proposed to learn sparsifying transforms, this end‐to‐end approach allows (1) to learn the “regularization” parameters (i.e., the thresholds) as well, which would otherwise have to be tuned with an expensive grid search and (2) to jointly learn stacked transform matrices which would otherwise have to be optimized separately in a greedy fashion.

Image decomposition

Since the noise in our phase data corrupts the entire frequency spectrum, applying the steps above solely to the input image will only attenuate the high‐frequency noise, leaving all low‐frequency noise intact. This is because the lowest frequency value accessible to the pipeline above is pixels (the patch size). Therefore, in our particular case, it is important to process the images at multiple scales. After applying our proposed pipeline to the input image , we thus blur the denoised output image with a 2D Gaussian kernel of radius (or size) and standard deviation and then downsample it by a factor of 2. has been chosen so that of the downsampled area is covered by the kernel, was sufficiently large coverage for the kernel. The resulting image is again denoised and then downsampled. This process is iteratively repeated three times, thereby generating downscaled and denoised versions of the input image for (see Figure 5). We stopped at because successive layers did not improve denoising performance, that is, the mean squared error (MSE) did not further improve (see Figure 7).

FIGURE 7

MSE of the images in Figure 6. Different scales (rows in Figure 6) are plotted in different colors for clarity. The biggest improvement in MSE happens at the finest scales, with less and less improvement when going to coarser scales

FIGURE 6

Denoising progression within INSIDEnet. The outputs after every second filter across the four scales are shown. Each row corresponds to one scale. From left to right, subsequent filters are shown. Only every second filter is shown for better visualization. Starting at the top left and moving to the bottom right, we observe that the image quality steadily improves. This is also confirmed by the results in Figure 7. As expected, the filters in the upper row only remove the high‐frequency noise. The lower frequency noise is removed at the lower scales (subsequent rows)

Let be the denoising pipeline illustrated in the previous section. and (used below) are downsampling and upsampling operators, respectively, and is a zero‐centered Gaussian smoothing filter with . The image decomposition of our model is then summarized by Algorithm 1.

It is important to note here that the denoising at multiple scales happens sequentially. While it could also be applied in parallel, this led to less accurate results. This is also more intuitive, since it seems redundant to denoise the same frequencies at more than one scale. As shown in Figure 5, this sequential process resembles the encoder of the U‐net. In the U‐net, however, the filtering is not as interpretable since the filtering steps are not based on an established image filtering technique such as transform domain thresholding, but instead apply convolutions followed by a rather arbitrary activation function. The highly nonlinear and complicated decoder of the U‐net is replaced in our model by a simple Laplacian pyramid assembling.

U‐Nets are state‐of‐the‐art models for image‐to‐image problems, and our approach opens up the possibility of being competitive with U‐Nets using an interpretable algorithm.

Image reconstruction

Once all images have been denoised at different scales, they are used to reconstruct the final image. The approach has been inspired by Burger and Harmeling and consists of the steps displayed in Algorithm 2 and Figure 5.

At scale , we separate the low‐frequency (blue in Figure 5) and high‐frequency components (yellow in Figure 5). is obtained by first smoothing and downsampling the image, thus removing the high‐frequency details, followed by upsampling and smoothing. The high frequencies in contrast are simply obtained by subtracting the low frequencies from the image. Combining the best of two adjacent scales, that is, the low frequencies from the lower scale and the high frequencies from the upper scale, we obtain , a better version of . By starting with the coarsest two scales and iteratively applying these steps, the final image is assembled.

Loss function and optimization

We used the MSE on the full‐resolution image as a loss function and propagated back all gradients with respect to both the transform matrices and the thresholds:

We deployed the Adam optimization algorithm with an exponentially decaying learning rate (initial learning rate of 0.0001) and trained all models with a batch size of 1 (because of GPU memory constraints) until convergence, that is until the MSE on the validation set was not improving anymore.

Algorithm interpretability and nonexpansiveness

Our algorithm is essentially composed of linear filters and thresholding. By reshaping the learned transforms into filter kernels, it is possible to get insights about what types of features are being used to build sparse data representations. Likewise, the learned threshold values indicate how strongly the images get filtered.

While the inspection of filter kernels is also possible in CNNs, and the bias terms can give some hints about filtering strength, CNNs do not allow to easily inspect how the images are being denoised within the network. This is because of the immense number of filter channels that are applied at each layer, and which are not projected back into the image space after each filtering step. In contrast, our model allows us to easily visualize the results in the intermediate layers of the network as we iteratively move back and forth between image space and transform space. By looking at Figures 6 and  7, we observe that our model gradually improves image quality, both visually and quantitatively. From this plot, it also emerges that four scales and 10 filtering steps are good hyperparameters for this application.

Moreover, by design, our model starts training with an excellent parameter initialization. In fact, owing to its peculiar architecture and by setting the starting thresholds to , prior to training, a forward pass through the network will leave the image unchanged. This is in strong contrast with conventional models such as the U‐net in which a first forward pass will yield a very different output than the image target. This implies that smaller parameter adjustments need to be made during training in the INSIDEnet as compared to the U‐net.

What further sets our model apart from most state‐of‐the‐art denoising networks is that the linear transformation is explicitly constrained to be orthogonal, thereby increasing the stability of the network. In fact, it has been shown that orthogonality is sufficient for a 1‐Lipschitz and nonexpansiveness property which in turn makes networks more robust to adversarial attacks. Moreover, nonexpansiveness is an important property when it comes to proximal operators in plug‐and‐play frameworks, where the presented model is envisioned to be used.

All in all, our model resembles the U‐net architecture. However, each part in our architecture has a clear mathematical rationale which comes with desirable properties and the possibility to better look inside the denoising process, thus making our approach both more interpretable and more reliable compared to standard CNNs.

Computational aspects

The algorithm has been implemented in Tensorflow 2.1, and all computations have been carried out on an NVIDIA Titan RTX GPU with 24 GB of memory. Processing of a single slice takes an average of 0.15 s, depending on the hyperparameters and . This is three orders of magnitude faster than the BM3D filter and approximately as fast as a deep CNN.

RESULTS

Simulated data

To assess the effectiveness of the INSIDEnet, we performed a comparative study on the simulated breast phantoms (absorption and phase). We compared the INSIDEnet with a classical state‐of‐the‐art algorithm, namely the BM3D filter, as well as with a deep CNN. We performed various experiments by changing the hyperparameters and and found and to work well. We would like to mention that by using a larger , less image scales would be necessary as large patches allow to denoise lower frequencies compared to small patches. However, this would go at the expense of a much higher parameter number needed to transform large patches. Therefore, to keep a smaller parameter number for higher generalizability, it is recommended to use a small patch size with more image scales. All performance metrics have been computed on 440 testing images and are provided in Table 1. To be able to compute SNR and contrast‐to‐noise ratio (CNR) efficiently over all 440 slices, we calculated mean and standard deviations in regions of interest where the gray level values in the ground truth image were approximately constant.

TABLE 1

Denoising results

Model	nMAE	SNR	CNR	SSIM
Phase
Input	0.030 (0.002)	29.776 (1.998)	2.905 (0.193)	0.861 (0.023)
INSIDEnet	0.017 (0.003)	64.305 (5.944)	6.179 (0.470)	0.962 (0.012)
BM3D	0.017 (0.002)	53.420 (4.380)	5.208 (0.417)	0.965 (0.010)
U‐net	0.015 (0.002)	72.919 (7.095)	6.934 (0.480)	0.973 (0.007)
Absorption
Input	0.059 (0.004)	15.143 (1.111)	2.975 (0.207)	0.694 (0.014)
INSIDEnet	0.017 (0.002)	62.463 (6.473)	11.663 (0.995)	0.926 (0.011)
BM3D	0.018 (0.002)	63.103 (6.515)	11.742 (0.948)	0.924 (0.010)
U‐net	0.016 (0.002)	64.541 (6.378)	11.945 (0.958)	0.939 (0.008)

Note: Standard deviations of the metrics across all 440 slices are displayed in parentheses.

Abbreviations: CNR, contrast‐to‐noise ratio; Interpretable NonexpanSIve Data‐Efficient network; INSIDEnet, nMAE, normalised mean absolute error; SNR, signal‐to‐noise ratio; SSIM, structural similarity index.

We used the BM3D filter of the bm3d.py software package. The sigma value was obtained by computing the standard deviation in a uniform region of the noisy input image. Processing of a single pixels slice took 76 s on a CPU.

As a deep CNN, we implemented a U‐net with 269 176 parameters trained with am MSE loss in Tensorflow 2.1 which, as for the INSIDEnet, separately processes either phase or absorption images. The parameter number has been kept relatively small to be approximately in the same range as the 166 403 parameters of the INSIDEnet model and to avoid overfitting. We used an initial learning rate of 0.0001 with exponential decay, along with the Adam optimization algorithm (). Before entering the network, the input images were brought to zero mean and unit variance. As for our proposed algorithm, the training set consisted of 440 pairs of noisy and clean simulated breast phantom slices. We would like to point out that we purposely chose a rather small training set, as this will be our real‐world scenario once we will start acquiring real data. Therefore, we wanted to investigate the algorithm's capability to generalize well from little data. Early stopping has been used to arrest the training when the validation loss was not improving anymore. The processing time for a single slice during prediction was 0.06 s.

The denoising results of the three algorithms on phase contrast phantoms are shown in Figure 8. We display the results on a zoomed‐in part of the image for better visualization. In the top row, a noisy image is shown along with its clean counterpart. In the other rows, the performance of the three algorithms is shown, along with the difference image of the denoised and the input data. We see that the INSIDEnet is able to effectively remove noise across the frequency spectrum while keeping sharp edges. This is also supported by the difference image in which no signal, but a large range of frequencies are present. The third row shows that, as our proposed model, also the U‐net is able to satisfactorily denoise the data. The performance of the BM3D filter in contrast is much inferior: it is unsurprisingly unable to remove low‐frequency noise owing to its small patch‐based strategy. In fact, the difference image shows only high‐frequency noise. A quantitative comparison between the three models (see the upper part of Table 1) reveals that the proposed model is only slightly inferior compared to the U‐net which, in the absence of any architectural constraints, achieves the best results. In agreement with a visual inspection, the BM3D model performs significantly worse in terms of SNR and CNR. Interestingly, in terms of normalised mean absolute error (nMAE) and structural similarity index (SSIM), the BM3D algorithm is competitive with the two data‐driven pipelines. This shows the limitations of such metrics as visually the superior performance of the latter methods is evident.

FIGURE 8

Denoising result on the phase contrast phantom obtained with the INSIDEnet, the U‐net, and the BM3D filter. Top left: noisy input image, top right: clean image. In subsequent rows left: denoised image, on the right: difference image between noisy input and denoised image. Gray value units are phase shift coefficients . The black scale bar is 5 mm

By looking at the denoising performance of the three algorithms on absorption data in Figure 9, we see a different pattern than in phase. In particular, owing to the different noise spectrum, concentrated in the high frequencies, all three filters achieve a satisfying result. The metrics in the lower part of Table 1 confirm that the three algorithms achieve a very similar performance.

FIGURE 9

Denoising result on the absorption contrast phantom obtained with the INSIDEnet, the U‐net, and the BM3D filter. Top left: noisy input image, top right: clean image. In subsequent rows left: denoised image, on the right: difference image between noisy input and denoised image. Gray value units are attenuation coefficients . The black scale bar is 5 mm

As expected, a comparison between denoising results on phase and absorption data reveals that it is easier to denoise the latter data. In fact, all image quality metrics improve more significantly in absorption than in phase.

Data efficiency

To assess the data‐driven models' generalization performance, we performed an experiment in which we trained both architectures on very limited data, that is, a single image and tested them on 440 slices as in the previous experiment. The results in Table 2 show that under such conditions the INSIDEnet quantitatively outperforms the U‐net. A visual inspection of the results in Figure 10 confirms these findings. In fact, a close look at the images reveals that the U‐net blurs the phase image and introduces dark artifacts around some edges in the absorption image.

TABLE 2

Denoising results with a model trained on a single image pair

	nMAE	SNR	CNR	SSIM
Phase
INSIDEnet	0.018 (0.003)	60.347 (5.480)	5.832 (0.413)	0.956 (0.012)
U‐net	0.021 (0.002)	52.369 (5.241)	4.928 (0.333)	0.961 (0.009)
Absorption
INSIDEnet	0.017 (0.002)	61.459 (6.313)	11.515 (0.973)	0.924 (0.011)
U‐net	0.020 (0.002)	53.149 (8.013)	9.265 (0.874)	0.924 (0.010)

Note: Standard deviations of the metrics across all 440 slices are displayed in parentheses.

Abbreviations: CNR, contrast‐to‐noise ratio; Interpretable NonexpanSIve Data‐Efficient network; INSIDEnet, nMAE, normalised mean absolute error; SNR, signal‐to‐noise ratio; SSIM, structural similarity index.

FIGURE 10

Illustration of higher data efficiency of the INSIDEnet compared to the U‐net. The models have been trained on a single training image. The black scale bar is 5 mm

These results suggest that (1) the INSIDEnet architecture imposes a strong inductive bias on the denoising problem, thereby enabling it to fit its parameters with very limited data; and (2) that the INSIDEnet's efficient initialization strategy indeed helps to more easily fit the model. These two aspects clearly set our model apart from classical CNNs. This has a great practical significance as training data is always scarce in a clinical setting and especially in the development phase of a clinical prototype. Therefore, having a model which can effectively be trained with small amounts of clean images is crucial.

Lesion detectability task

To assess our algorithm's performance in terms of lesion detectability, we added identical lesions with different contrasts at random locations to our simulated testing data. In Figure 11, five lesions with varying contrast (1.02, 1.04, 1.06, 1.08, and 1.10) are shown for the clean phase contrast phantom data, for the corresponding noisy data as well as for all algorithms considered in this paper. We see that none of the algorithms is able to recover the lesion below a contrast of 1.06. The BM3D filter leads to relatively good lesion delineation, despite its inability to remove the low‐frequency noise. The INSIDEnet and the U‐net achieve a comparable performance when trained on 440 images. A close look reveals that the INSIDEnet seems more robust compared to the U‐net as it leads to slightly better lesion delineation than the U‐net, when trained on a single image.

FIGURE 11

Lesion detectability task on phase contrast data. First row: phantom data, second row: noisy data, third row: denoised with INSIDEnet, fourth row: denoised with INSIDEnet trained on a single image, fifth row: denoised with U‐net, sixth row: denoised with U‐net trained on a single image, seventh row: denoised with BM3D. From left to right the lesion contrast increases (1.02, 1.04, 1.06, 1.08, 1.10). The displayed area is 7 mm 10 mm

By looking at Figure 12 for detectability on absorption, we see a slightly different result. The INSIDEnet and the BM3D filter enable lesion delineation down to the contrast of 1.06, with the latter achieving the best overall performance. Somewhat surprisingly, the U‐net is unable to retrieve the lesion at contrast 1.06. As for the phase, we see that the INSIDEnet is more robust compared to the U‐net when trained on a single image. In fact, the former show very little loss in performance when trained on a single image.

FIGURE 12

Lesion detectability task on absorption contrast data. First row: phantom data, second row: noisy data, third row: denoised with INSIDEnet, fourth row: denoised with INSIDEnet trained on a single image, fifth row: denoised with U‐net, sixth row: denoised with the U‐net trained on a single image, seventh row: denoised with BM3D. From left to right, the lesion contrast increases (1.02, 1.04, 1.06, 1.08, 1.10). The displayed area is 7 mm 10 mm

Figure 13 shows the meat sample scanned on our GI‐BCT prototype in phase and absorption contrast. It is easy to see that indeed in absorption the noise is concentrated in the high frequencies, whereas in phase also lower frequencies are affected, thus making denoising more challenging. However, as it emerges from Table 3, phase contrast data actually have both a higher SNR as well as a higher CNR than in absorption contrast, thus highlighting why there is such a high interest in bringing phase contrast to clinical practice.

FIGURE 13

Phase contrast (left) and absorption (right) contrast of meat scan acquired on our GI‐BCT prototype. The red boxes have been used to compute SNR and CNR values in Table 3. The white scale bars are 15 mm

TABLE 3

Denoising results on real data

	Phase		Absorption
Model	SNR	CNR	SNR	CNR
Input	25.607	3.618	13.996	1.869
INSIDEnet	57.967	8.332	82.777	13.540
BM3D	38.161	5.391	79.956	12.321
U‐net	63.949	8.785	86.270	14.278

Abbreviations: CNR, contrast‐to‐noise ratio; Interpretable NonexpanSIve Data‐Efficient network; INSIDEnet, nMAE, normalised mean absolute error; SNR, signal‐to‐noise ratio; SSIM, structural similarity index.

However, on this particular sample there does not appear to be a qualitative advantage of phase contrast data over absorption as no extra information appears to be present compared to absorption.

Figure 14 shows the denoising results of the two data‐driven models as well as BM3D. As on the simulations, on phase contrast both data‐driven models achieve a satisfactory performance, whereas BM3D fails to effectively denoise the tomograms. A close look suggests though that the U‐net slightly blurs the image as it was the case when we trained the model on a single slice. This might suggest that the U‐net tends to blur images when it has to generalize to unseen data. Quantitatively, the U‐net achieves the highest performance. BM3D performs significantly worse on both metrics. On absorption, all methods achieve a satisfying performance. The U‐net achieves the highest SNR and CNR, closely followed by both the INSIDEnet and the BM3D filter. The line profiles in Figure 15 indicate that all models are able to maintain sharp edges, except for the U‐net on phase contrast data.

FIGURE 14

Denoising results on a scan of a chunk of meat. The red lines (10 mm) show the region at which the line profiles in Figure 15 have been computed

FIGURE 15

Line profiles taken at the locations displayed in Figure 14

The fact that the models could be trained on simulations and applied to real measurements suggests that both image and noise statistics of the simulations match well to real measurements. This is an important finding as it allows us to efficiently train models in the absence of ground truth.

DISCUSSION

This paper is part of a larger effort to translate GI‐BCT to the clinics and provides physicians and patients with a technology which can offer higher tissue contrast and thus increased chances of spotting potential malignancies. To achieve this, data‐processing algorithms are needed which can handle the low raw data quality as it is currently acquired in GI‐BCT. The goal of this paper was to investigate denoising of reconstructed tomograms. We used two established methods (BM3D and U‐net), proposed a new dedicated method (INSIDEnet), and tested their ability to reliably deal with the complex noise pattern in GI's phase contrast channel.

Where the BM3D filter was unable to satisfactorily denoise phase contrast images, we could show that it is possible to significantly increase phase contrast image quality with the two data‐driven methods. The U‐net yielded the best quantitative results when trained on relatively large amounts of data, whereas the proposed model has proven to be very data efficient. In fact, a single image sufficed to obtain a satisfactory denoising performance. The power of combining data‐ and model‐based methods has thus been highlighted by the INSIDEnet in outperforming the U‐net when trained on very limited data.

The INSIDEnet has been designed to find a good trade‐off between performance and interpretability/robustness, crucial in the medical field. It is possible to inspect all of its trainable weights which are interpretable as they (1) indicate how strong the sparsity enforcing regularization is or (2) show the filter kernels used to achieve sparse data representation. Furthermore, as shown in Figure 6, one can visualize the progress the model makes in its deeper layers, thus offering superior interpretability compared to CNNs. The steady image quality improvement that emerges in this regard supports our claim of high model reliability.

We started developing this algorithm by leveraging classical ideas in image denoising. Curiously, the final architecture resembles the U‐net in many ways. First of all, the multiscale processing is very similar to the “U” structure of the U‐net. Second, our model also consists of linear matrix multiplications, followed by nonlinear activations. In our case, however, the activation function is not a simple Rectified linear unit (ReLU) function but is rather mathematically motivated by the proximal operator of the norm. Finally, the concatenation of the transpose of the current transform and the next transform could be parameterized with a single matrix without losing expressive power. However, we believe that the inductive bias realized by going back to image space after each thresholding operation, by leveraging explicit orthogonality, instead of applying solely feed‐forward matrix multiplications, is the key feature that gives our model higher robustness compared to the U‐net.

A crucial part of our algorithm is thus given by the explicit orthogonality constraint enforced with the Cayley transform. , In fact, with no such constraint, model training and performance were much more unstable. While such explicit constraints could also be implemented in a convolutional layer, this would require significantly higher computational costs, which are avoided by enforcing kernel orthogonality on transform matrices operating on image patches.

We would like to further point out that, besides providing more than satisfactory results, our algorithm is fast, which is critical in our project as we must process large image volumes. Finally, while the algorithm has been developed to cope with phase contrast CT, it can be applied to a variety of image data.

Finally, we could show that real data could be efficiently denoised by our model, even if trained on simulated data, thus confirming the validity of our simulation study. Unfortunately, both simulations and real data have shown that with the current setup we are unable to generate phase contrast data which are superior to absorption, neither prior to denoising nor after. Future efforts will thus be oriented towards improving grating quality to reduce the noise propagation in the phase contrast images.

CONCLUSION

In conclusion, we have shown that it is possible to efficiently denoise both simulated and real noisy phase contrast images with a data‐efficient, fast, and interpretable data‐driven algorithm. We expect the INSIDEnet to become an important tool as part of a dedicated plug‐and‐play iterative reconstruction framework, which is currently under development. We hope that in future, together with improved grating quality, this will allow us to effectively deal with the low‐frequency noise of the phase contrast channel and let its higher intrinsic contrast emerge.

CONFLICT OF INTEREST

The authors have no conflicts to disclose.