Ultra Low-Parameter Denoising: Trainable Bilateral Filter Layers in Computed
Tomography
Fabian Wagner1, Mareike Thies1, Mingxuan Gu1, Yixing Huang1, Sabrina
Pechmann2, Mayank Patwari1, Stefan Ploner1, Oliver Aust3, Stefan Uderhardt4,5,
Georg Schett4,5, Silke Christiansen2,6, Andreas Maier1
1Pattern Recognition Lab, Friedrich-Alexander University Erlangen-Nürnberg,
Germany
2Fraunhofer Institute for Ceramic Technologies and Systems IKTS, Forchheim,
Germany
3Leibniz Institute for Analytical Sciences ISAS, Dortmund, Germany
4Department of Internal Medicine 3 - Rheumatology and Immunology, Friedrich-
Alexander University Erlangen-Nürnberg, Germany
5University Hospital Erlangen, Germany
6Institute for Nanotechnology and Correlative Microscopy e.V. INAM, Forchheim,
Germany

Author to whom correspondence should be addressed. email<EMAIL_ADDRESS>

###### Abstract

Background: Computed tomography is widely used as an imaging tool to visualize
three-dimensional structures with expressive bone-soft tissue contrast.
However, CT resolution and radiation dose are tightly entangled, highlighting
the importance of low-dose CT combined with sophisticated denoising
algorithms.
Purpose: Most data-driven denoising techniques are based on deep neural
networks and, therefore, contain hundreds of thousands of trainable
parameters, making them incomprehensible and prone to prediction failures.
Developing understandable and robust denoising algorithms achieving state-of-
the-art performance helps to minimize radiation dose while maintaining data
integrity.
Methods: This work presents an open-source CT denoising framework based on the
idea of bilateral filtering. We propose a bilateral filter that can be
incorporated into a deep learning pipeline and optimized in a purely data-
driven way by calculating the gradient flow toward its hyperparameters and its
input. Denoising in pure image-to-image pipelines and across different domains
such as raw detector data and reconstructed volume, using a differentiable
backprojection layer, is demonstrated.
Results:Although only using three spatial parameters and one range parameter
per filter layer, the proposed denoising pipelines can compete with deep
state-of-the-art denoising architectures with several hundred thousand
parameters. Competitive denoising performance is achieved on x-ray microscope
bone data (0.7053 and 33.10) and the 2016 Low Dose CT Grand Challenge dataset
(0.9674 and 43.07) in terms of SSIM and PSNR.
Conclusions: Due to the extremely low number of trainable parameters with
well-defined effect, prediction reliance and data integrity is guaranteed at
any time in the proposed pipelines, in contrast to most other deep learning-
based denoising architectures.

##  I. Introduction

Ionizing radiation used in Computed Tomography (CT) is known to be harmful to
any living tissue. Therefore, the amount of deposited energy in each
investigated sample, the so-called dose, is to be minimized following the As
Low As Reasonably Achievable (ALARA) principle 1. However, radiation dose and
image resolution are inherently entangled as a certain number of x-rays must
be absorbed by the scanned patient to achieve the desired contrast 2. Further
reducing the patient dose results in degraded image quality due to increased
Poisson noise in CT projections.
Since the emergence of the first CT scanners, denoising algorithms were
developed and applied to restore image quality while keeping the radiation
dose moderate 3. Non-linear filters 4, 5, 6, 7 and iterative reconstruction
techniques 8, 9, 10 have been successfully applied. Whereas such approaches
usually require hand-tuned hyperparameters, cannot abstract complex features,
or are known to be computationally expensive, purely data-driven, end-to-end
trainable, approaches have been proposed in recent years 11, 12, 13, 14, 15,
16, 17 fueled by the emergence of deep learning and, in particular,
convolutional neural networks. Most of these models achieve competitive
denoising performance but are built on deep neural networks with multiple
layers and contain hundreds of thousands of trainable parameters. Although
deep architectures help networks extract complex features from data, such
models are often regarded as black boxes as it is impossible to fully
comprehend their data processing and control failing network predictions.
Besides, adversarial examples in the form of small perturbations of the
network input can lead to undesired drastic changes in the prediction due to
the high dimension of extracted features 18. Such uncertainties often prohibit
deep learning applications in the medical imaging field where the reliability
of the data must be maintained 19, 20. Additionally, training large numbers of
parameters in neural networks usually requires broad medical datasets with
paired ground truth data, which can be hard to obtain.
The bilateral filter has been successfully applied on CT data 6 as it performs
combined filtering in spatial and range domain with Gaussian kernels smoothing
the noise fluctuations while preserving edges 21. However,
hyperparameters—fundamentally determining the filter performance—usually have
to be hand-picked. Multiple works have been proposed aiming for automatically
finding optimal filter parameters 22, 23, 24, 25, 26 using risk estimators or
grasshopper optimization. All these techniques are not suitable for
integration into an end-to-end, data-driven, optimization pipeline and are
based on different statistical assumptions on the noise. A different work
presents a bilateral filter based on convolutional filtering the permutohedral
lattice, a higher-dimensional representation of the image data, that can be
optimized 27. However, their algorithm requires so-called splatting and
slicing operations, which can only approximate the data and, therefore,
introduce uncertainties. Additionally, the requirement of a higher-dimensional
grid increases computational demands. Patwari et al. 12 proposed the JBFnet, a
deep neural network architecture inspired by joint bilateral filtering that
can be trained in a purely data-driven way. They showed competitive
performance to much deeper networks, although reducing the number of trainable
parameters by choosing a shallow convolutional architecture. Besides,
denoising approaches indirectly optimizing hyperparameters of (joint)
bilateral filters using reinforcement learning have been proposed 28, 29.
However, their training is more sophisticated as it is again based on deep
architectures for choosing the correct parameter updates. Additionally,
Patwari et al. 29 make use of a residual network to generate an image prior as
well as a reward network.
We aim to extend the aforementioned approaches by proposing a trainable
bilateral filter layer using only the inherent four filter parameters
(including three spatial filter dimensions) during training and inference. By
analytically deriving the gradient flow toward the parameters as well as to
the layer input, we can directly optimize all filter parameters via
backpropagation and incorporate the C++/CUDA accelerated layer into any
trainable pipeline using the PyTorch framework 30. Additionally, we show
simultaneous optimization of filter parameters in the projection and image
domain, using a differentiable backprojection layer 31. Due to the very low
number of trainable parameters, we can optimize our pipeline with only very
little training data and in a self-supervised manner using Noise2Void training
32, while still achieving competitive performance compared to state-of-the-art
deep neural networks. We explain the competitive denoising performance with
the theoretical findings of Maier et al. 33, 34, who proved that incorporating
prior knowledge into artificial neural networks lowers the upper error bound
of the model prediction. The experiments are performed on the 2016 Low Dose CT
Grand Challenge dataset 35 ($25\,\%$ dose) for benchmark purposes as well as
on a low-dose x-ray microscope (XRM) bone dataset ($10\,\%$ dose).

##  II. Materials and Methods

###  II.A. Trainable bilateral filter

The bilateral filter had great success in CT denoising due to its ability to
smooth image content while preserving edges. This is achieved by a composed
filter kernel with spatial and range contributions. With the noisy input
reconstruction $\mathbf{X}$ and the denoised prediction $\mathbf{\hat{Y}}$,
the discrete filter operation can be written as 21

$\displaystyle\hat{Y}_{k}=\frac{1}{w_{k}}\underbrace{\sum_{n\in\mathcal{N}}G_{\sigma_{s}}(\mathbf{p_{k}}-\mathbf{p_{n}})G_{\sigma_{r}}(X_{k}-X_{n})X_{n}}_{\begin{subarray}{c}=:\,\alpha_{k}\end{subarray}}$
(1)

with the definition of the normalization factor $w_{k}$

$\displaystyle
w_{k}:=\sum_{n\in\mathcal{N}}G_{\sigma_{s}}(\mathbf{p_{k}}-\mathbf{p_{n}})G_{\sigma_{r}}(X_{k}-X_{n}),$
(2)

voxel index $k\in\mathbb{N}$, and the Gaussian function

$\displaystyle
G_{\sigma_{r}}(c):=\exp\left(-\frac{c^{2}}{2\sigma_{r}^{2}}\right).$ (3)

Each predicted output voxel $\hat{Y}_{k}$ is calculated from the neighborhood
$\mathcal{N}$ of voxel $X_{k}$, indexed by $n\in\mathcal{N}$. The spatial
kernel $G_{\sigma_{s}}$ performs image smoothing and is dependent on the
distance between the positions $\mathbf{p_{k}}\in\mathbb{N}^{d}$ and
$\mathbf{p_{n}}\in\mathbb{N}^{d}$. In the $d$-dimensional case
$G_{\sigma_{s}}$ is composed of multiple Gaussian kernels, e.g.,

$\displaystyle
G_{\sigma_{s}}(\mathbf{c})=\prod_{\beta\in\\{x,y,z\\}}\exp\left(-\frac{c_{\beta}^{2}}{2\sigma_{\beta}^{2}}\right)$
(4)

assuming $d=3$ spatial filter dimensions and the spatial distance vector
$\mathbf{c}=(c_{x},c_{y},c_{z})$. The additional range kernel $G_{\sigma_{r}}$
is responsible for preserving edges during filtering as similar voxel values
are stronger weighted by incorporating the intensity distance $(X_{k}-X_{n})$.

Figure 1: Working scheme of the proposed trainable bilateral filter layer.
Black arrows mark the forward pass, whereas red arrows illustrate the gradient
flow toward the input $X_{i}$ and the filter parameters $\sigma_{\gamma}$
during backpropagation.

Kernel widths $\sigma_{s}$ and $\sigma_{r}$ are hyperparameters of the
bilateral filter and are usually hand-picked by the user dependent on image
content and voxel range. However, tuning the filter layers by hand is
cumbersome, and finding an optimal parameter set can not be guaranteed. To
optimize the parameters in a data-driven fashion and incorporate them in a
fully differentiable pipeline, the gradient, given by the derivative of the
loss function $L$ with respect to each parameter $\sigma_{\gamma}$ must be
calculated

$\displaystyle\frac{\partial L}{\partial\sigma_{\gamma}}=\frac{\partial
L}{\partial\mathbf{\hat{Y}}}\frac{\partial\mathbf{\hat{Y}}}{\partial\sigma_{\gamma}}=\sum_{k}\frac{\partial
L}{\partial\hat{Y}_{k}}\frac{\partial\hat{Y}_{k}}{\partial\sigma_{\gamma}}.$
(5)

Automatically deriving the gradient with a framework like PyTorch is in
principle possible, but comes with a huge computational and memory overhead.
It would require calculating gradients for each computation in the filter
forward pass, which is computationally infeasible with conventional hardware
in a reasonable run time. Instead, we directly calculate the analytical
solution of $\frac{\partial\hat{Y}_{k}}{\partial\sigma_{\gamma}}$, which is
described in detail in the Appendix. Additionally, the gradient with respect
to every voxel of the noisy input reconstruction $X_{i}$ needs to be derived
to propagate a loss into previous filter layers during backpropagation to
allow stacking multiple bilateral filters or incorporating the layer into a
deep architecture

$\displaystyle\frac{\partial L}{\partial X_{i}}=\frac{\partial
L}{\partial\mathbf{\hat{Y}}}\frac{\partial\mathbf{\hat{Y}}}{\partial
X_{i}}=\sum_{k}\frac{\partial
L}{\partial\hat{Y}_{k}}\frac{\partial\hat{Y}_{k}}{\partial X_{i}}.$ (6)

With each predicted $\hat{Y}_{k}$ being dependent on intensity differences of
two input voxels $X_{k}$ and $X_{n}$, the analytical calculation of the
gradient flow toward the filter input

$\displaystyle\frac{\partial\hat{Y}_{k}}{\partial
X_{i}}=-w_{k}^{-2}\alpha_{k}\frac{\partial w_{k}}{\partial
X_{i}}+w_{k}^{-1}\frac{\partial\alpha_{k}}{\partial X_{i}}$ (7)

is more complex than the derivative with respect to the kernel widths
$\sigma_{\gamma}$. Eventually, it requires distinguishing the two cases $k=i$
and $k\neq i$ for finding an applicable solution. They refer to the off-center
and center elements of the bilateral filter kernel respectively. The detailed
calculation is again conducted in the Appendix, accompanied by implementation
notes. Fig. 1 illustrates the data flow in a bilateral filter layer. The
forward and backward pass of the layer is implemented in C++ and CUDA to
leverage performance and integrated into the PyTorch framework using the
pybind11 module 36. Our open-source filter layer is publicly available under
https://github.com/faebstn96/trainable-bilateral-filter-source.

###  II.B. Denoising pipelines

Figure 2: Illustration of the employed denoising pipelines containing the
proposed bilateral filter layers. Multiple filters are trained in a supervised
(first row) and unsupervised (Noise2Void 32, second row) fashion.
Additionally, bilateral filter layers are optimized in projection and image
domain simultaneously using a differentiable backprojection operator (dual-
domain BFs, third row). The $\sigma_{\gamma}$ parameters represent the only
trainable parameters of the networks and are all optimized independently. ND
refers to the normal-dose reconstructions whereas LD denotes the low-dose data
incorporating noise.

Multiple denoising pipelines containing trainable bilateral filter layers,
illustrated in Fig. 2, are employed to investigate the performance and
limitations of the layers. First, a simple post-processing approach is
studied, directly applying three subsequent bilateral filter layers on the
reconstructed volume. The filters are optimized for the mean squared error
(MSE) loss calculated between prediction $\mathbf{\hat{Y}}$ and normal-dose
volume $\mathbf{Y}$.
Additionally, self-supervised training is performed using the Noise2Void (N2V)
training scheme. The noisy volume is perturbed by randomly replacing a defined
ratio of voxels ($1\,\%$) with voxel intensities from their respective $5^{3}$
voxel neighborhoods and fed through the denoising pipeline. Subsequently, an
MSE loss is calculated between prediction and noisy non-perturbed volume at
the voxel positions of the perturbations. If the noise contribution of
neighboring voxels is approximately uncorrelated, the resulting gradient will
converge to the denoised version of the volume, as proved by Krull et al. 32.
The advantage of the N2V training scheme is that no paired high-dose data is
required. Accordingly, no additional ground truth knowledge can be employed
during training.
Finally, an FDK-based reconstruction pipeline with a single bilateral filter
layer in the projection and two filters in the image domain is trained end-to-
end on XRM data. The employed pipeline is illustrated in the last row of Fig.
2 and makes use of the fully differentiable reconstruction pipeline of Thies
et al. 37, which is described in the following section. All filter parameters
are optimized for the MSE loss between $\mathbf{\hat{Y}}$ and the normal-dose
reconstruction $\mathbf{Y}$ in the image domain.

###  II.C. XRM reconstruction pipeline

The FDK-based 38 reconstruction pipeline by Thies et al. 37 connects acquired
CT projection data with its 3D interpretable reconstruction. The
backprojection is incorporated as a known operator33, using the differentiable
layer provided by the Pyro-NN framework by Syben et al. 31. The reconstruction
pipeline is adapted to XRM projection data from a Zeiss Xradia 620 Versa
microscope, which is a high-resolution cone-beam CT for small samples, and
allows to propagate a loss calculated in the image domain to any location
within the reconstruction pipeline enabled by its differentiable
implementation. Accordingly, bilateral filter layers can be trained at
multiple locations in the pipeline in a purely data-driven manner using the
PyTorch framework. Particularly, we perform experiments demonstrating the
superior performance of bilateral filter layers applied in both projection and
image domain simultaneously over denoising in image domain only.

###  II.D. Experimental setup

We compare the performance of the trainable bilateral filter layer-based
pipelines to four deep state-of-the-art denoising architectures on the ten
abdomen scans from the 2016 Low Dose CT Grand Challenge dataset 35 with
$1\,\text{mm}$ slice thickness reconstructed in slices of size $512\times
512$. Here, patient L291 is selected for validation, patients L310, L333 and
L506 for testing and all $3411$ slices from the remaining seven scans for
training.
In a second experiment, we investigate the denoising performance on a dataset
of high-resolution ex vivo XRM scans of mouse tibia bones. In the context of
preclinical osteoporosis research, lacunae, bone structures with a diameter of
3-20 µm, were investigated 39. This is instructive for understanding the bone
metabolism 40, 41 as they contain osteocytes, cells being heavily involved in
the bone remodeling process 42, 43. However, the severe amount of radiation
dose required for resolving micrometer-sized lacunae in a mouse so far
prohibits in vivo imaging 44, 45 and, therefore, requires sophisticated dose
reduction techniques to push the XRM in vivo limit to the micrometer scale. We
create a low-dose XRM dataset from ex-vivo mouse tibia bone samples to
investigate in-vivo XRM acquisitions. The sample preparation was carried out
as approved by the local animal ethics committee in compliance with all
ethical regulations (license TS-10/2017). $1400$ high-resolution projections
are acquired for every sample in a short-scan setting with $28\,\text{s}$
exposure time per projection image. Low-dose and normal-dose projections are
reconstructed in a $2048^{3}$ grid of $1.4\,\text{\textmu m}$ isotropic voxel
size using the pipeline of Thies et al. 37. Data acquired with $10\,\%$ dose
is simulated by incorporating noise following Yu et al. 46. Two separate
stacks of $30$ slices each are reconstructed from every bone scan as the total
amount of reconstructed slices would by far exceed the capacity of the
denoising models and the variety of learnable local features. The final XRM
bone dataset contains five scans for training, one for validation, and four
for testing, which is in total equivalent to $9600$ $512\times 512$ image
patches, exceeding the size of the 2016 Low Dose CT Grand Challenge dataset.
The quadratic autoencoder architecture (QAE) proposed by Fan et al. 13, as
well as the JBFnet by Patwari et al. 12, both achieve remarkable results on
the CT denoising task, outperforming many architectures proposed during the
2016 Low Dose CT Grand Challenge. The RED-CNN network, published by Chen et
al. 11 is a deep, convolutional architecture with more than $1.8\,\text{Mio}$
parameters achieving state-of-the-art results likewise. Another denoising
approach based on a generative adversarial network is WGAN 14, using a
combination of Wasserstein and perceptual loss. Encouraged by their
performance, these four methods are selected for comparison with multiple
pipelines containing the proposed trainable bilateral filter layers on the
aforementioned datasets.
In order to allow a fair comparison between the different denoising
approaches, the networks are implemented and trained based on officially
published code repositories and the description from their publications. The
Adam optimizer with learning rates $4\cdot 10^{-4}$ (decaying), $1\cdot
10^{-4}$ (constant), $1\cdot 10^{-4}$ (decaying), and $1\cdot 10^{-5}$
(constant) is used for QAE, JBFnet, RED-CNN, and WGAN respectively. Both QAE
and RED-CNN are trained using the MSE loss, whereas JBFnet makes use of a
custom loss described in their publication 12 as well as pre-training of its
prior module. Generator and discriminator of WGAN are trained alternately
using Wasserstein loss for stable convergence and a perceptual loss as
presented in 14. All models are trained until the validation loss plateaus or
increases.
In all bilateral filter-based models, $\sigma_{x,y,z}$ and $\sigma_{r}$ are
initialized with $0.5$ and $0.01$ and optimized with two separate Adam
optimizers with learning rates $0.01$ and $0.005$ respectively, as the spatial
and range parameters operate on two independent scales. The tiny amount of
required training data is demonstrated by training the bilateral filter
pipelines only on a single stack of $21$ neighboring slices with a size of
$512\times 512$ voxels from one training scan of the Grand Challenge dataset.
For the XRM bone dataset a stack of $15$ neighboring patches of size
$512\times 512$ voxels is used. Note that the bilateral filter layers cannot
overfit the data assuming a comparable amount of noise within all scans due to
the low number of trainable parameters with well-defined influence.
Optimization is performed until convergence of the training loss for up to
$5000$ iterations which took up to $20$ minutes (on an NVIDIA Quadro RTX
4000).

Figure 3: Example input and target slices from the investigated Grand
Challenge dataset ($25\,\%$ dose) and the XRM bone dataset ($10\,\%$ dose).
Red squares highlight the ROIs for the visual comparison. The display windows
are $[-150,250]\,\text{HU}$ and $[0.25,0.7]\,\text{arb.}\,\text{unit}$.

##  III. Results

###  III.A. Denoising results

Figure 4: The first and third rows show denoising predictions from the 2016 Low Dose CT Grand Challenge dataset in the ROIs illustrated in Fig. 3 for all employed methods. Our proposed pipelines are highlighted in bold letters. The display window is $[-150,250]\,\text{HU}$. The second and fourth rows present difference images between prediction and target where prediction artifacts like blurred edges become particularly visible. Difference images are plotted in the window $[-150,150]\,\text{HU}$. Table 1: Quantitative denoising results of the compared pipelines on the test patients of the 2016 Low Dose CT Grand Challenge dataset. For each model, the number of trainable parameters and the average inference time per $512\times 512$ image slice (Quadro RTX 4000) are presented. The names of our proposed pipelines and the best-performing models are highlighted.  | Grand Challenge Dataset |  | $[\text{s}]$
---|---|---|---
|  SSIM (mean $\pm$ std) |  PSNR (mean $\pm$ std) | # params | runtime
Low-dose CT | $0.8991\pm 0.033$ | $38.64\pm 1.62$ | - | -
WGAN 14 | $0.9341\pm 0.085$ | $40.47\pm 0.86$ | $7,800,000$ | $0.1$
RED-CNN 11 | $0.9680\pm 0.010$ | $43.63\pm 1.19$ | $1,800,000$ | $4.2$
JBFnet 12 | $0.9666\pm 0.012$ | $42.77\pm 1.13$ | $118,000$ | $5.6$
QAE 13 | $0.9650\pm 0.010$ | $43.35\pm 1.18$ | $137,000$ | $2.0$
1 BF layer | $0.9656\pm 0.013$ | $42.73\pm 1.38$ | $4$ | $0.2$
3 BF layers | $0.9674\pm 0.012$ | $43.07\pm 1.42$ | $12$ | $0.5$
3 BFs (N2V 32) | $0.9577\pm 0.015$ | $41.15\pm 1.17$ | $12$ | $0.5$
Table 2: Quantitative denoising results of the compared pipelines on the test patients of the XRM bone dataset. For each model, the number of trainable parameters is shown in the right column.  | XRM Bone Dataset |
---|---|---
|  SSIM (mean $\pm$ std) |  PSNR (mean $\pm$ std) | # parameters
Low-dose CT | $0.1652\pm 0.0084$ | $19.95\pm 0.27$ | -
WGAN 14 | $0.5170\pm 0.0135$ | $29.45\pm 0.23$ | $7,800,000$
RED-CNN 11 | $0.7122\pm 0.0087$ | $33.09\pm 0.20$ | $1,800,000$
JBFnet 12 | $0.6564\pm 0.0139$ | $30.98\pm 0.31$ | $118,000$
QAE 13 | $0.7126\pm 0.0088$ | $33.09\pm 0.20$ | $137,000$
1 BF layer | $0.6599\pm 0.0167$ | $31.40\pm 0.30$ | $4$
3 BF layers | $0.6698\pm 0.0158$ | $32.19\pm 0.26$ | $12$
3 BF layers (N2V 32) | $0.4590\pm 0.0280$ | $27.61\pm 0.57$ | $12$
Dual-domain BFs | $0.7053\pm 0.0120$ | $33.10\pm 0.22$ | $12$

We present qualitative denoising results on selected slices from both test
datasets, visualized in Fig. 3. Multiple magnified regions of interest (ROIs)
allow comparing the denoising performance on small image features of the Grand
Challenge dataset for better visualization and are shown in Fig. 4. The
additionally provided difference images (prediction - target) particularly
highlight artifacts in the model predictions. Closer studying local features
reveals over-smoothed results in low-contrast regions for the QAE and the RED-
CNN compared to the normal-dose target, highlighted by orange circles and
arrows. However, high-frequency details like edges are well preserved as there
are few structures visible in the difference images. The predictions of
bilateral filter-based pipelines, JBFnet, and WGAN visually appear more
natural and closer to the target images, compared to the other approaches.
However, blurred edges are visible in the difference images for JBFnet, WGAN,
and the self-supervised bilateral filter pipeline (3BFs N2V). The predictions
of the supervised-trained bilateral filter pipeline (3BFs) contain the noise
pattern closest to the target ROI, while slightly stronger removing edges
compared to QAE and RED-CNN. Quantitatively, RED-CNN, 3BFs, and QAE outperform
the other models on the Grand Challenge dataset in terms of structural
similarity index measure (SSIM) and peak signal-to-noise ratio (PSNR),
presented in Table 1.
The model predictions on the XRM bone dataset are depicted in Fig. 5,
including the entire reconstruction pipeline, simultaneously denoising in
sinogram and reconstruction domain (Dual-domain BFs). Visually, Dual-domain
BFs, RED-CNN, and QAE outperform the other models in terms of noise removal
and edge sensitivity. Predictions of the self-supervised bilateral filter
pipeline preserve edges but still contain a substantial amount of noise. In
contrast, the difference images of WGAN and JBFnet reveal more removed high-
frequency details at edges, compared to all other methods. Quantitatively,
denoising with three trained bilateral filters in the reconstruction domain
performs slightly worse than QAE and RED-CNN in terms of SSIM and PSNR, as
listed in Table 2. However, moving one of the bilateral filter layers into the
sinogram domain improves the denoising performance to match the SSIM and PSNR
of the best-performing deep models. All trainable bilateral filter-based
pipelines use several orders of magnitude fewer parameters compared to the
deep reference architectures. Simultaneously, quantitatively and qualitatively
competitive denoising performance is achieved, outperforming multiple deep
reference methods.

Figure 5: The first and third rows show denoising predictions from the XRM
bone dataset in the ROIs illustrated in Fig. 3 for all employed methods. Our
proposed pipelines are highlighted in bold letters. The display window is
$[0.25,0.7]\,\text{arb.}\,\text{unit}$. The second and fourth rows present
difference images between prediction and target where prediction artifacts
like blurred edges become particularly visible. Difference images are plotted
in the window $[-0.1,0.1]\,\text{arb.}\,\text{unit}$.

###  III.B. Depth of the bilateral filter pipeline

Figure 6: The denoising performance of multiple stacked bilateral filter
layers is compared based on the average SSIM and PSNR on the test patients of
the 2016 Low Dose CT Grand Challenge dataset. Denoising with two or more
trained layers is advantageous compared to only using a single layer.

The benefit of iteratively removing noise with subsequent bilateral filters is
investigated by training bilateral filter pipelines of different depths.
Pipelines containing up to six filter layers are trained until convergence and
tested on the Grand Challenge dataset. Mean SSIM and PSNR suggest that the
denoising performance indeed benefits from stacking three filter layers, as
shown in Fig. 6. Adding more filters slightly reduces the performance as more
parameters have to properly converge to relatively small values to preserve
edges.

##  IV. Discussion

The proposed pipelines quantitatively and qualitatively show comparable
denoising to state-of-the-art deep architectures. Although information from a
single CT projection is spread through the whole volume, the appearance of
noise in the reconstructed volume is a local phenomenon. Smart filtering
within a finite neighborhood can, therefore, perform surprisingly well, as is
the case for the bilateral filter. However, especially the choice of optimal
range parameters $\sigma_{r}$ is very crucial for the denoising performance.
If $\sigma_{r}$ is chosen too large, regions of constant attenuation are well
restored, but edges are blurred concurrently, leading to degraded predictions.
Optimizing the hyperparameters in a purely data-driven way can, therefore,
leverage the applicability of the bilateral filter, as well as guarantee a
near-optimal choice of parameters, as we demonstrated by achieving
quantitative scores comparable to state-of-the-art methods.
We empirically found that the denoising performances of optimized pipelines
are independent of the parameter initialization. For pipelines containing more
than a single filter layer, the filter parameters first converge to very
similar values between the layers until the loss is almost minimized. During
further training, individual filter parameters start to converge to distinct
values, slowly reducing the training loss further. During this fine-tuning
phase, stacked bilateral filters learn to focus on different features, e.g.,
preserving edges or smoothing planes.
Compared to the ground truth targets, the predictions of RED-CNN and QAE
visually appear over-smoothed in low-contrast regions, which can remove
features being beneficial for the physician. Whereas supervised training of
the bilateral filter-based pipelines performs comparable to deep reference
models on both investigated datasets, the self-supervised approach using the
Noise2Void technique does not fully converge to the optimal parameter set. In
general, improved image quality can be recognized, however, the presented
difference images identify blurred edges due to insufficiently learned filter
parameters. If noise is not fully distinguishable from content, the range
parameters $\sigma_{r}$ do not properly converge to the optimum values. In
their work, Krull et al. 32 already identified this limitation of their
Noise2Void algorithm for images with inherently high irregularities.

Figure 7: Example of failing network predictions of the trained deep QAE model
on strongly absorbing objects in the scan of patient L006 and L193 from the
Low Dose CT Image and Projection Dataset 47. The display window is
$[-150,250]\,\text{HU}$.

Approximation and estimation error are commonly combined to estimate the risk
bound of neural networks. The two contributions refer first to the distance
between target function and closest neural network function and second to the
distance between the estimated and the ideal network function 48. Our
quantitative and visual results show that deep neural networks and trainable
bilateral filter pipelines can both well predict denoised reconstructions
after extensive training and, therefore, show comparable approximation risk.
However, we think trainable bilateral filter layers have a considerably
smaller estimation error due to the restricted influence of the few trained
parameters. The decreased error bound can be explained by treating the filter
algorithm similar to a known operator 33, 34. In contrast, deep neural
network-based approaches can learn network functions that only work well
within the training data distribution but fail to produce reliable results for
other samples from the CT data distribution, which results in a larger
estimation error.
In a clinical environment, data integrity at any point in the medical data
processing pipeline is crucial for a reliable diagnosis from the physician.
During training, deep models are optimized on a finite number of training
samples to extract features generalizing over the entire distribution of data.
However, in a clinical environment, data can be acquired which is not properly
represented by the finite training data distribution. Image processing
algorithms must be able to handle such samples to avoid failing downstream
tasks, including the diagnosis by the physician. Fig. 7 exhibits failing
network predictions of one deep reference method on out-of-distribution CT
slices around medical implants in abdomen scans of the Low Dose CT Image and
Projection Dataset 47. The reference method was before trained on the 2016 Low
Dose CT Grand Challenge training dataset until convergence and achieves state-
of-the-art performance on the test data, as presented in Table 1. A different
predicted artifact of a deep reference method is highlighted in Fig. 4 with
red arrows. Here, a shadow-like structure right next to the bone region is
introduced which is also visible in the difference image. By algorithmic
design, trainable bilateral filter layers only filter in a finite neighborhood
with well-defined Gaussian kernels and can not produce such artifacts in their
predictions. The inherent lack of complex feature abstraction in the proposed
filter layer enforces proximity between prediction and measured ground truth
and is advantageous for maintaining reliable data compared to the
nontransparent data processing performed in deep state-of-the-art denoising
models. Additionally, the low amount of required training data makes trainable
bilateral filter layers convenient to use and easy to integrate into complex
data processing pipelines.

##  V. Conclusions

We present a trainable, fully differentiable bilateral filter layer that can
be directly incorporated in any deep architecture using the PyTorch framework
with GPU acceleration. We show that trainable bilateral filter layers with
only four parameters can achieve state-of-the-art performance of deep neural
networks with hundred of thousands of trainable parameters on the low-dose CT
denoising task, while performing fully comprehensible data processing and
producing physically reliable predictions at any time. The low number of
parameters allows training with little data, self-supervised training using
the Noise2Void scheme, and incorporating the filter layer at multiple
locations in a CT reconstruction pipeline to further leverage denoising
performance. In summary, trainable bilateral filters allow extensive dose
reduction while maintaining high image quality without sacrificing data
integrity.

##  VI. Acknowledgments

The research leading to these results has received funding from the European
Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation program (ERC Grant No. 810316).

##  VII. Conflict of Interest

The authors have no relevant conflicts of interest to disclose.

## Appendix I: Analytical derivative of the bilateral filter

We aim to make the spatial and range hyper parameters $\sigma_{s}$ and
$\sigma_{r}$ of the bilateral filter trainable in a gradient descent
algorithm. In order to update the parameters during optimization the
derivative of the loss function $L$ with respect to each parameter
$\sigma_{\gamma}$ is required

$\displaystyle\frac{\partial L}{\partial\sigma_{\gamma}}=\frac{\partial
L}{\partial\mathbf{\hat{Y}}}\frac{\partial\mathbf{\hat{Y}}}{\partial\sigma_{\gamma}}=\sum_{k}\frac{\partial
L}{\partial\hat{Y}_{k}}\frac{\partial\hat{Y}_{k}}{\partial\sigma_{\gamma}}.$
(8)

Following the definition of the bilateral filter in (1), the derivative
$\frac{\partial\hat{Y}_{k}}{\partial\sigma_{\gamma}}$ yields

$\displaystyle\frac{\partial\hat{Y}_{k}}{\partial\sigma_{\gamma}}=-w_{k}^{-2}\alpha_{k}\frac{\partial
w_{k}}{\partial\sigma_{\gamma}}+w_{k}^{-1}\frac{\partial\alpha_{k}}{\partial\sigma_{\gamma}}$
(9)

with

$\displaystyle\frac{\partial w_{k}}{\partial\sigma_{\gamma}}$
$\displaystyle=\sum_{n\in\mathcal{N}}\frac{\partial}{\partial\sigma_{\gamma}}G_{\sigma_{s}}(\mathbf{p_{k}}-\mathbf{p_{n}})G_{\sigma_{r}}(X_{k}-X_{n}),$
(10) $\displaystyle\frac{\partial\alpha_{k}}{\partial\sigma_{\gamma}}$
$\displaystyle=\sum_{n\in\mathcal{N}}X_{n}\frac{\partial}{\partial\sigma_{\gamma}}G_{\sigma_{s}}(\mathbf{p_{k}}-\mathbf{p_{n}})G_{\sigma_{r}}(X_{k}-X_{n}),$
(11)

and the derivative of the Gaussian kernel

$\displaystyle\frac{\partial}{\partial\sigma_{\gamma}}G_{\sigma}(c)=G_{\sigma}(c)\frac{c^{2}}{\sigma_{\gamma}^{3}}.$
(12)

Additionally, the gradient with respect to every voxel of the noisy input
volume $X_{i}$ needs to be derived, as it should be possible to propagate a
loss into previous filter layers during backpropagation

$\displaystyle\frac{\partial L}{\partial X_{i}}=\frac{\partial
L}{\partial\mathbf{\hat{Y}}}\frac{\partial\mathbf{\hat{Y}}}{\partial
X_{i}}=\sum_{k}\frac{\partial
L}{\partial\hat{Y}_{k}}\frac{\partial\hat{Y}_{k}}{\partial X_{i}}.$ (13)

The derivative of the output voxel $\hat{Y}_{k}$ with respect to the input
voxel $X_{i}$ yields

$\displaystyle\frac{\partial\hat{Y}_{k}}{\partial
X_{i}}=-w_{k}^{-2}\alpha_{k}\frac{\partial w_{k}}{\partial
X_{i}}+w_{k}^{-1}\frac{\partial\alpha_{k}}{\partial X_{i}}$ (14)

by again using the definition of the bilateral filter from (1) and applying
the product rule analog to (9).
The filtered output voxel $\hat{Y}_{k}$ is dependent on both input voxels
$X_{k}$ and $X_{n}$. To carry out the derivative with respect to $X_{i}$ the
two cases $k\neq i$ and $k=i$ are distinguished.

Case 1: ($k\neq i$)

$\displaystyle\begin{split}\left.\frac{\partial w_{k}}{\partial
X_{i}}\right|_{k\neq
i}&=\,\sum_{n\in\mathcal{N}}G_{\sigma_{s}}(\mathbf{p_{k}}-\mathbf{p_{n}})\frac{\partial}{\partial
X_{i}}G_{\sigma_{r}}(X_{k}-X_{n})\\\
&=\,G_{\sigma_{s}}(\mathbf{p_{k}}-\mathbf{p_{i}})G_{\sigma_{r}}(X_{k}-X_{i})\frac{X_{k}-X_{i}}{\sigma_{r}^{2}}\end{split}$
(15) $\displaystyle\begin{split}\left.\frac{\partial\alpha_{k}}{\partial
X_{i}}\right|_{k\neq
i}=&\,\sum_{n\in\mathcal{N}}G_{\sigma_{s}}(\mathbf{p_{k}}-\mathbf{p_{n}})\frac{\partial}{\partial
X_{i}}G_{\sigma_{r}}(X_{k}-X_{n})X_{n}\\\
=&\,G_{\sigma_{s}}(\mathbf{p_{k}}-\mathbf{p_{i}})\left[\left(\frac{\partial}{\partial
X_{i}}G_{\sigma_{r}}(X_{k}-X_{i})\right)X_{i}+G_{\sigma_{r}}(X_{k}-X_{i})\left(\frac{\partial}{\partial
X_{i}}X_{i}\right)\right]\\\
=&\,G_{\sigma_{s}}(\mathbf{p_{k}}-\mathbf{p_{i}})\cdot
G_{\sigma_{r}}(X_{k}-X_{i})\left[\frac{X_{k}-X_{i}}{\sigma_{r}^{2}}X_{i}+1\right]\end{split}$
(16)

In both expressions only the $i$th term of the sum ($n=i$) contributes.

Case 2: ($k=i$)

$\displaystyle\begin{split}\left.\frac{\partial w_{k}}{\partial
X_{i}}\right|_{k=i}=&\,\frac{\partial}{\partial
X_{i}}\sum_{n\in\mathcal{N}}G_{\sigma_{s}}(\mathbf{p_{i}}-\mathbf{p_{n}})G_{\sigma_{r}}(X_{i}-X_{n})\\\
=&\,\frac{\partial}{\partial
X_{i}}\left[\underbrace{G_{\sigma_{s}}(\mathbf{p_{i}}-\mathbf{p_{i}})G_{\sigma_{r}}(X_{i}-X_{i})}_{\begin{subarray}{c}=1\end{subarray}}+\sum_{n\in\mathcal{N},n\neq
i}G_{\sigma_{s}}(\mathbf{p_{i}}-\mathbf{p_{n}})G_{\sigma_{r}}(X_{i}-X_{n})\right]\\\
=&\,\sum_{n\in\mathcal{N},n\neq
i}G_{\sigma_{s}}(\mathbf{p_{i}}-\mathbf{p_{n}})\frac{\partial}{\partial
X_{i}}G_{\sigma_{r}}(X_{i}-X_{n})\\\
=&\,\sum_{n\in\mathcal{N}}G_{\sigma_{s}}(\mathbf{p_{i}}-\mathbf{p_{n}})G_{\sigma_{r}}(X_{i}-X_{n})\frac{X_{n}-X_{i}}{\sigma_{r}^{2}}\end{split}$
(17) $\displaystyle\begin{split}\left.\frac{\partial\alpha_{k}}{\partial
X_{i}}\right|_{k=i}=&\,\frac{\partial}{\partial
X_{i}}\sum_{n\in\mathcal{N}}G_{\sigma_{s}}(\mathbf{p_{i}}-\mathbf{p_{n}})G_{\sigma_{r}}(X_{i}-X_{n})X_{n}\\\
=&\,\frac{\partial}{\partial
X_{i}}\left[\underbrace{G_{\sigma_{s}}(\mathbf{p_{i}}-\mathbf{p_{i}})G_{\sigma_{r}}(X_{i}-X_{i})}_{\begin{subarray}{c}=1\end{subarray}}X_{i}+\sum_{n\in\mathcal{N},n\neq
i}G_{\sigma_{s}}(\mathbf{p_{i}}-\mathbf{p_{n}})G_{\sigma_{r}}(X_{i}-X_{n})X_{n}\right]\\\
=&\,1+\sum_{n\in\mathcal{N},n\neq
i}\left[G_{\sigma_{s}}(\mathbf{p_{i}}-\mathbf{p_{n}})\cdot
G_{\sigma_{r}}(X_{i}-X_{n})X_{n}\frac{X_{n}-X_{i}}{\sigma_{r}^{2}}\right]\\\
=&\,1+\sum_{n\in\mathcal{N}}\left[G_{\sigma_{s}}(\mathbf{p_{i}}-\mathbf{p_{n}})\cdot
G_{\sigma_{r}}(X_{i}-X_{n})X_{n}\frac{X_{n}-X_{i}}{\sigma_{r}^{2}}\right]\end{split}$
(18)

Note that in the last steps the sum can be carried out over the entire
neighborhood $N$ as the contribution of the term $n=i$ is zero respectively.

## Appendix II: Implementation

In practice, (13) describes a convolution of the kernel
$\frac{\partial\hat{Y}_{k}}{\partial X_{i}}$ with the backpropagated loss
$\frac{\partial L}{\partial\hat{Y}_{k}}$ in a neighborhood denoted by $k$
which is given by the finite kernel size—analog to a conventional
convolutional layer. However, due to the dependency of the kernel on two
voxels of the input volume $X_{k}$ and $X_{i}$, the analytical derivative of
the convolutional kernel is more elaborate. When calculating the sum in (13),
both (15) and (16) are used for the contributions $k\neq i$ and only the term
$k=i$ is derived in (17) and (18). Note that for faster computation the sums
in (17) and (18) can be pre-calculated in the forward pass of the filtering.
The finite spatial kernel is calculated by taking into account voxels covering
five times the spatial standard deviation in each spatial dimension, but at
least $5\times 5\times 5$ voxels. This could affect the performance for large
spatial kernels with $\sigma_{s}\gg 1$, but was never a relevant factor in any
of our experiments.
To validate the correct implementation of the filter derivative, we compared
the analytical gradient, provided by the backward pass of the filter layer,
with a numerical approximation of the gradient, computed via small finite
differences. Our public repository contains a gradient check script using the
torch.autograd.gradcheck function from the PyTorch framework.

## References

##

  * 1 J. M. Boone, W. R. Hendee, M. F. McNitt-Gray, and S. E. Seltzer,  Radiation exposure from CT scans: how to close our knowledge gaps, monitor and safeguard exposure—proceedings and recommendations of the Radiation Dose Summit, sponsored by NIBIB, February 24–25, 2011,  Radiology 265, 544–554 (2012).
  * 2 H. H. Barrett, S. Gordon, and R. Hershel,  Statistical limitations in transaxial tomography,  Computers in biology and medicine 6, 307–323 (1976).
  * 3 A. Maier and R. Fahrig,  GPU denoising for computed tomography,  Graphics Processing Unit-Based High Performance Computing in Radiation Therapy 113 (2015).
  * 4 Z. S. Kelm, D. Blezek, B. Bartholmai, and B. J. Erickson,  Optimizing non-local means for denoising low dose CT,  in 2009 IEEE International Symposium on Biomedical Imaging: From Nano to Macro, pages 662–665, IEEE, 2009.
  * 5 A. Maier, L. Wigström, H. G. Hofmann, J. Hornegger, L. Zhu, N. Strobel, and R. Fahrig,  Three-dimensional anisotropic adaptive filtering of projection data for noise reduction in cone beam CT,  Med. Phys. 38, 5896–5909 (2011).
  * 6 A. Manduca, L. Yu, J. D. Trzasko, N. Khaylova, J. M. Kofler, C. M. McCollough, and J. G. Fletcher,  Projection space denoising with bilateral filtering and CT noise modeling for dose reduction in CT,  Med. Phys. 36, 4911–4919 (2009).
  * 7 M. Manhart, R. Fahrig, J. Hornegger, A. Doerfler, and A. Maier,  Guided noise reduction for spectral CT with energy-selective photon counting detectors,  in Proceedings of the Third CT Meeting, pages 91–94, 2014.
  * 8 X. Han, J. Bian, D. R. Eaker, T. L. Kline, E. Y. Sidky, E. L. Ritman, and X. Pan,  Algorithm-enabled low-dose micro-CT imaging,  IEEE Transactions on Medical Imaging 30, 606–620 (2010).
  * 9 M. Beister, D. Kolditz, and W. A. Kalender,  Iterative reconstruction methods in X-ray CT,  Physica Medica 28, 94–108 (2012).
  * 10 P. Gilbert,  Iterative methods for the three-dimensional reconstruction of an object from projections,  Journal of theoretical biology 36, 105–117 (1972).
  * 11 H. Chen, Y. Zhang, M. K. Kalra, F. Lin, Y. Chen, P. Liao, J. Zhou, and G. Wang,  Low-Dose CT With a Residual Encoder-Decoder Convolutional Neural Network,  IEEE Transactions on Medical Imaging 36, 2524–2535 (2017).
  * 12 M. Patwari, R. Gutjahr, R. Raupach, and A. Maier,  JBFnet - Low Dose CT Denoising by Trainable Joint Bilateral Filtering,  in International Conference on Medical Image Computing and Computer-Assisted Intervention—MICCAI 2020, pages 506–515, Springer, 2020\.
  * 13 F. Fan, H. Shan, M. K. Kalra, R. Singh, G. Qian, M. Getzin, Y. Teng, J. Hahn, and G. Wang,  Quadratic autoencoder (Q-AE) for low-dose CT denoising,  IEEE Transactions on Medical Imaging 39, 2035–2050 (2019).
  * 14 Q. Yang, P. Yan, Y. Zhang, H. Yu, Y. Shi, X. Mou, M. K. Kalra, Y. Zhang, L. Sun, and G. Wang,  Low-Dose CT Image Denoising Using a Generative Adversarial Network With Wasserstein Distance and Perceptual Loss,  IEEE Transactions on Medical Imaging 37, 1348–1357 (2018).
  * 15 E. Kang, W. Chang, J. Yoo, and J. C. Ye,  Deep convolutional framelet denosing for low-dose CT via wavelet residual network,  IEEE Transactions on Medical Imaging 37, 1358–1369 (2018).
  * 16 E. Kang, J. Min, and J. C. Ye,  A deep convolutional neural network using directional wavelets for low-dose X-ray CT reconstruction,  Med. Phys. 44, e360–e375 (2017).
  * 17 M. D. Ketcha, M. Marrama, A. Souza, A. Uneri, P. Wu, X. Zhang, P. A. Helm, and J. H. Siewerdsen,  Sinogram + image domain neural network approach for metal artifact reduction in low-dose cone-beam computed tomography,  Journal of Medical Imaging 8, 1–16 (2021).
  * 18 X. Yuan, P. He, Q. Zhu, and X. Li,  Adversarial Examples: Attacks and Defenses for Deep Learning,  IEEE Transactions on Neural Networks and Learning Systems 30, 2805–2824 (2019).
  * 19 V. Antun, F. Renna, C. Poon, B. Adcock, and A. C. Hansen,  On instabilities of deep learning in image reconstruction and the potential costs of AI,  Proceedings of the National Academy of Sciences 117, 30088–30095 (2020).
  * 20 Y. Huang, T. Würfl, K. Breininger, L. Liu, G. Lauritsch, and A. Maier,  Some investigations on robustness of deep learning in limited angle tomography,  in International Conference on Medical Image Computing and Computer-Assisted Intervention—MICCAI 2018, pages 145–153, Springer, 2018\.
  * 21 C. Tomasi and R. Manduchi,  Bilateral filtering for gray and color images,  in Sixth international conference on computer vision, pages 839–846, IEEE, 1998.
  * 22 Y. Chen and Y. Shu,  Optimization of bilateral filter parameters via chi-square unbiased risk estimate,  IEEE Signal Processing Letters 21, 97–100 (2013).
  * 23 V. Anoop and P. R. Bipin,  Medical image enhancement by a bilateral filter using optimization technique,  Journal of medical systems 43, 1–12 (2019).
  * 24 H. Kishan and C. S. Seelamantula,  Optimal parameter selection for bilateral filters using Poisson Unbiased Risk Estimate,  in 2012 19th IEEE International Conference on Image Processing, pages 121–124, IEEE, 2012.
  * 25 T. Dai, Y. Zhang, L. Dong, L. Li, X. Liu, and S. Xia,  Content-Aware Bilateral Filtering,  in 2018 IEEE Fourth International Conference on Multimedia Big Data (BigMM), pages 1–6, IEEE, 2018.
  * 26 H. Peng and R. Rao,  Bilateral kernel parameter optimization by risk minimization,  in 2010 IEEE International Conference on Image Processing, pages 3293–3296, IEEE, 2010.
  * 27 V. Jampani, M. Kiefel, and P. V. Gehler,  Learning sparse high dimensional filters: Image filtering, dense crfs and bilateral neural networks,  in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 4452–4461, 2016.
  * 28 W. Kang and M. Patwari,  Low Dose Helical CBCT denoising by using domain filtering with deep reinforcement learning,  arXiv preprint arXiv:2104.00889 (2021).
  * 29 M. Patwari, R. Gutjahr, R. Raupach, and A. Maier,  Low Dose CT Denoising via Joint Bilateral Filtering and Intelligent Parameter Optimization,  arXiv preprint arXiv:2007.04768 (2020).
  * 30 A. Paszke et al.,  PyTorch: An Imperative Style, High-Performance Deep Learning Library,  in Proc. NeurIPS, pages 8024–8035, Curran Associates, Inc., 2019\.
  * 31 C. Syben, M. Michen, B. Stimpel, S. Seitz, S. Ploner, and A. K. Maier,  PYRO-NN: Python reconstruction operators in neural networks,  Med. Phys. 46, 5110–5115 (2019).
  * 32 A. Krull, T.-O. Buchholz, and F. Jug,  Noise2void-learning denoising from single noisy images,  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 2129–2137, 2019.
  * 33 A. Maier, C. Syben, B. Stimpel, T. Würfl, M. Hoffmann, F. Schebesch, W. Fu, L. Mill, L. Kling, and S. Christiansen,  Learning with known operators reduces maximum error bounds,  Nature machine intelligence 1, 373–380 (2019).
  * 34 A. Maier, F. Schebesch, C. Syben, T. Würfl, S. Steidl, J.-H. Choi, and R. Fahrig,  Precision learning: towards use of known operators in neural networks,  in 2018 24th International Conference on Pattern Recognition (ICPR), pages 183–188, IEEE, 2018.
  * 35 C. H. McCollough et al.,  Low-dose CT for the detection and classification of metastatic liver lesions: results of the 2016 low dose CT grand challenge,  Med. Phys. 44, e339–e352 (2017).
  * 36 W. Jakob,  pybind11 – Seamless operability between C++11 and Python, 2021.
  * 37 M. Thies,  Differentiable reconstruction for x-ray microscopy data,  Code Ocean (2021).
  * 38 L. A. Feldkamp, L. C. Davis, and J. W. Kress,  Practical cone-beam algorithm,  Josa a 1, 612–619 (1984).
  * 39 K. M. Hannah, C. D. Thomas, J. G. Clement, F. De Carlo, and A. G. Peele,  Bimodal distribution of osteocyte lacunar size in the human femoral cortex as revealed by micro-CT,  Bone 47, 866–871 (2010).
  * 40 A. Grüneboom et al.,  A network of trans-cortical capillaries as mainstay for blood circulation in long bones,  Nature Metabolism 1, 236–250 (2019).
  * 41 A. Grüneboom, L. Kling, S. Christiansen, L. Mill, A. Maier, K. Engelke, H. H. Quick, G. Schett, and M. Gunzer,  Next-generation imaging of the skeletal system and its blood supply,  Nature Reviews Rheumatology 15, 533–549 (2019).
  * 42 R. Gruber, P. Pietschmann, and M. Peterlik,  Introduction to bone development, remodelling and repair,  in Radiology of osteoporosis, pages 1–23, Springer, 2008.
  * 43 P. R. Buenzli and N. A. Sims,  Quantifying the osteocyte network in the human skeleton,  Bone 75, 144–150 (2015).
  * 44 L. Mill, L. Kling, A. Grüneboom, G. Schett, S. Christiansen, and A. Maier,  Towards In-Vivo X-Ray Nanoscopy,  in Bildverarbeitung für die Medizin 2019, pages 251–256, Springer, 2019.
  * 45 Y. Huang, L. Mill, R. Stoll, L. Kling, O. Aust, F. Wagner, A. Grüneboom, G. Schett, S. Christiansen, and A. Maier,  Semi-permeable Filters for Interior Region of Interest Dose Reduction in X-ray Microscopy.,  in Bildverarbeitung für die Medizin, pages 61–66, 2021.
  * 46 L. Yu, M. Shiung, D. Jondal, and C. H. McCollough,  Development and validation of a practical lower-dose-simulation tool for optimizing computed tomography scan protocols,  Journal of computer assisted tomography 36, 477–487 (2012).
  * 47 T. R. Moen, B. Chen, D. R. Holmes III, X. Duan, Z. Yu, L. Yu, S. Leng, J. G. Fletcher, and C. H. McCollough,  Low-dose CT image and projection dataset,  Med. Phys. 48, 902–911 (2021).
  * 48 A. R. Barron,  Approximation and estimation bounds for artificial neural networks,  Machine learning 14, 115–133 (1994).