research papers\(\def\hfill{\hskip 5em}\def\hfil{\hskip 3em}\def\eqno#1{\hfil {#1}}\)

Journal logoJOURNAL OF
SYNCHROTRON
RADIATION
ISSN: 1600-5775

Ring artifact reduction via multiscale nonlocal collaborative filtering of spatially correlated noise

crossmark logo

aTampere University, Finland, and bSLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025, USA
*Correspondence e-mail: ymir.makinen@tuni.fi

Edited by M. Yamamoto, RIKEN SPring-8 Center, Japan (Received 13 July 2020; accepted 17 February 2021; online 16 April 2021)

X-ray micro-tomography systems often suffer severe ring artifacts in reconstructed images. These artifacts are caused by defects in the detector, calibration errors, and fluctuations producing streak noise in the raw sinogram data. In this work, these streaks are modeled in the sinogram domain as additive stationary correlated noise upon logarithmic transformation. Based on this model, a streak removal procedure is proposed where the Block-Matching and 3-D (BM3D) filtering algorithm is applied across multiple scales, achieving state-of-the-art performance in both real and simulated data. Specifically, the proposed fully automatic procedure allows for attenuation of streak noise and the corresponding ring artifacts without creating major distortions common to other streak removal algorithms.

1. Introduction

Ring artifacts are ubiquitous in computed tomography (Jha et al., 2013[Jha, A. K., Purandare, N. C., Shah, S., Agrawal, A., Puranik, A. D. & Rangarajan, V. (2013). Indian J. Nucl. Med. 28, 232-233.]; Artul, 2013[Artul, S. (2013). BMJ Case Rep. 2013, bcr-2013-201379.]; Boas & Fleischmann, 2012[Boas, F. E. & Fleischmann, D. (2012). Imaging Med. 4, 229-240.]); they originate from angular streak noise in measured raw sinogram data used to reconstruct a tomographic volume (Croton et al., 2019[Croton, L. C., Ruben, G., Morgan, K. S., Paganin, D. M. & Kitchen, M. J. (2019). Opt. Express, 27, 14231-14245.]) and appear as darker or lighter circles or arcs centered on the axis of rotation for data acquisition. Streak noise can be caused by mis-calibration of detector linear response, beam fluctuations, beam hardening, or dusty or damaged scintillator screens (Haibel, 2008[Haibel, A. (2008). Advanced Tomographic Methods in Materials Research and Engineering, Monographs on the Physics and Chemistry of Materials, pp. 141-160.]; Vidal et al., 2005[Vidal, F. P., Létang, J. M., Peix, G. & Cloetens, P. (2005). Nucl. Instrum. Methods Phys. Res. B, 234, 333-348.]; Anas et al., 2010[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.]).

Minimization of ring artifacts by using adequate scanning protocols (Pelt & Parkinson, 2018[Pelt, D. M. & Parkinson, D. Y. (2018). Meas. Sci. Technol. 29, 034002.]), high quality scintillator screens and detectors is possible. It is, however, difficult to completely avoid such artifacts and therefore achieve highest quality reconstruction solely by experimental measures. Several algorithms have been proposed to reduce ring artifacts in tomographic imaging, including wavelet-FFT filters (Münch et al., 2009[Münch, B., Trtik, P., Marone, F. & Stampanoni, M. (2009). Opt. Express, 17, 8567-8591.]), combinations of polynomial smoothing filters and careful calibration of the detector response function (Vo et al., 2018[Vo, N. T., Atwood, R. C. & Drakopoulos, M. (2018). Opt. Express, 26, 28396-28412.]; Croton et al., 2019[Croton, L. C., Ruben, G., Morgan, K. S., Paganin, D. M. & Kitchen, M. J. (2019). Opt. Express, 27, 14231-14245.]), or iterative algorithms (Paleo & Mirone, 2015[Paleo, P. & Mirone, A. (2015). J. Synchrotron Rad. 22, 1268-1278.]) that combine regularized reconstruction with denoising.

In this work, we model the streak noise as a spatially correlated noise in the sinogram domain, and propose a denoising procedure aiming to remove the streak noise before reconstruction. The denoising procedure is based on collaborative filtering, which employs both non-local self-similarity and transform-domain shrinkage to denoise a noisy signal through jointly transformed grouped blocks. In particular, we use the image denoising algorithm BM3D (Dabov et al., 2007[Dabov, K., Foi, A., Katkovnik, V. & Egiazarian, K. (2007). IEEE Trans. Image Process. 16, 2080-2095.], 2008[Dabov, K., Foi, A., Katkovnik, V. & Egiazarian, K. O. (2008). Proc. SPIE, 6812, 681207.]), leveraging the recent inclusion of exact transform-domain noise variances (Mäkinen et al., 2020[Mäkinen, Y., Azzari, L. & Foi, A. (2020). IEEE Trans. Image Process. 29, 8339-8354.]), which allow for accurate modeling of long noise correlation within the jointly transformed blocks.

Noting that some streaks may be too wide to be adequately captured by a group of standard-sized BM3D blocks, we further propose multiscale streak removal with BM3D. The proposed procedure is fully automatic and includes self-calibration of the filtering strength. We demonstrate the superior performance of the proposed approach on real data from the table-top Prisma XRM microCT at Sigray, and from the synchrotron-based microCT at the Advanced Photon Source (APS) in Argonne, available through Tomobank (De Carlo et al., 2018[De Carlo, F., Gürsoy, D., Ching, D. J., Batenburg, K. J., Ludwig, W., Mancini, L., Marone, F., Mokso, R., Pelt, D. M., Sijbers, J. & Rivers, M. (2018). Meas. Sci. Technol. 29, 034004.]).

2. Transform-domain collaborative filtering of correlated noise

In this section, we interpret sinogram streaks as spatially correlated noise, formalizing streak removal as filtering of correlated noise where the streaks follow a basic stationary model. As a powerful tool to deal with this model, we adopt a recent BM3D designed for dealing with long-range correlation such as that which characterizes the streaks. This constitutes the denoising module at the core of a multiscale and nonstationary filtering architecture that will be presented in Section 3[link] for the more general case of real-world streak noise.

2.1. Correlated noise model

We consider a noisy input [{z\!:\!X\!\rightarrow\!{\bb{R}}}] to be a combination of underlying data y and additive stationary spatially correlated noise η to be filtered,

[z(x) = y(x)+\eta(x),\qquad x\in X, \eqno(1)]

where [x\!\in\!X\!\subset\!{\bb{Z}}^{2}] is the coordinate in the finite two-dimensional image domain X (representing angles and displacements when z is a sinogram) and

[\eta = \nu\,{\bigcirc\kern-9pt\ast}\,\,g, \qquad \nu\left(\cdot\right)\sim{\cal N}\left(0,1 \right), \eqno(2)]

with ν being the zero-mean independent and identically distributed (i.i.d.) Gaussian noise with unit variance, and `[{\bigcirc\kern-9pt{*}}\,]' denoting 2-D convolution with the kernel g. The kernel g defines the spatial correlation of the noise as well as the noise strength, with [{\|g\|_{2} = {\rm{std}}(\eta)}]. An equivalent way of representing correlated noise is by its power spectral density (PSD) Ψ,

[\Psi = {\rm{E}}\left\{\left|{\cal{F}}\left[\eta\right]\right|^{2}\right\} = { \rm var}\left\{{\cal{F}}\left[\eta\right]\right\} = \left|X\right|\left|{ {\cal{F}}\left[g\right]}\right|^{2}, \eqno(3)]

with F being the 2-D Fourier transform and |X| denoting the cardinality (i.e. number of elements) of X. Equivalently, a kernel g satisfying (2[link])–(3[link]) can be defined from Ψ as

[g = |X|^{-1/2} \, {\cal{F}}^{\,-1}\big[{{\rm{std}}\left({\cal{F}}\left[\eta \right]\right)}\big] = |X|^{-1/2} \, {\cal{F}}^{\,-1}\left[{\sqrt{\Psi} }\right]. \eqno(4)]

2.1.1. Basic model for constant stationary streaks

Although the streaks are originally multiplicative in nature, sinogram data are considered upon a logarithmic transformation and therefore the streaks can be modeled by the additive noise η in (1[link]). The sinogram streak noise is fairly constant in the angular dimension, presenting very long-range correlation in the noise along this dimension. Treating angle as the vertical dimension and displacement as horizontal, we consider the basic case of horizontally white and vertically constant streak noise. Such noise can be modeled through (2[link]), when setting g as a vertically constant line. This simple kernel as well as the corresponding PSD Ψ are demonstrated in Figure 1[link]; an example of real streak noise viably approximated through this model is shown in Figure 2[link].

[Figure 1]
Figure 1
Example noise [\eta = \nu\,{\bigcirc\kern-8.5pt\ast}\,\,g], [\nu\left(\cdot\right)\sim{\cal N}\left(0,1\right)], the corresponding correlation kernel g, and PSD Ψ. For the kernel and the PSD, black pixels of the image correspond to value 0 in the data.
[Figure 2]
Figure 2
Example of streak noise in a sinogram (a fragment of Fly) that can be seen as well approximated by the model in Figure 1[link].

In practice, the above simple model cannot be used but for small segments of the sinograms, as streak noise can often feature horizontal correlation, vertical variations, or nonstationarities that are not described by the model. In Section 3[link], a complete processing pipeline is further proposed to allow modeling more complex cases of streak noise through (1[link])–(4[link]), enabling their attenuation through the collaborative filter.

2.2. Transform-domain collaborative filtering and BM3D

The rationale of transform-domain filtering is to work with a representation of the signal where most of the signal is compacted to only a few coefficients, whereas the remaining coefficients mostly comprise noise. Hence, by attenuating the coefficients with a non-linear shrinkage operator, it is possible to attenuate noise while keeping most of the signal intact. Nonlocal collaborative filters utilize this property in the context of collective transform coefficients of groups of similar blocks extracted from the image. In all of the following sections, we consider the recently proposed variant (Mäkinen et al., 2020[Mäkinen, Y., Azzari, L. & Foi, A. (2020). IEEE Trans. Image Process. 29, 8339-8354.]) of BM3D for correlated noise denoising where the input is z and the goal of denoising is to estimate y based on the statistics of η or equivalently knowledge of Ψ or g.

In BM3D, all operations are made with regard to a reference block moving through the image. For each position of the reference block, the following steps are executed:

(i) Collect similar blocks into a group through block-matching.

(ii) Obtain the 3-D transform spectrum by collectively transforming the obtained blocks.

(iii) Perform shrinkage.

(iv) Transform the shrunk spectra back to block estimates and aggregate them to the original locations from which they were collected.

The 3-D transform spectrum of the grouped noisy blocks is obtained through first applying a 2-D transform [{\cal T}^{\,2{\rm{D}}}] locally to each block, then a 1-D transform [{\cal T}^{\,{\rm{NL}}}] through the `stack' of grouped blocks. Denoting by [\left\{z_{x_{1}},\ldots,z_{x_{M}}\right\}] a group of M blocks of N pixels extracted from z at coordinates x1,…, xM, we obtain the [{\cal T}^{\,2{\rm{D}}}] spectrum coefficients as sixt = [\left\langle z_{x_{t}},b^{2{\rm{D}}}_{i}\right\rangle], for i = 1,…, N, t = 1,…, M, where [b^{2{\rm{D}}}_{i}] is the ith basis function of [{\cal T}^{\,2{\rm{D}}}]. The 3-D spectrum coefficients are calculated through the direct tensor product of the [{\cal T}^{\,2{\rm{D}}}] and [{\cal T}^{\,{\rm{NL}}}] transforms, as

[\eqalignno{ s_{i,j}^{{\,x_{1},\ldots,\,x_{M}}} & = \left\langle\left[z_{x_{1}}, \ldots, z_{x_ {M}}\right], b^{2{\rm{D}}}_{i} \!\otimes b^{{\rm{NL}}}_{j}\right\rangle \cr& = \left\langle\left[s_{i}^{\,x_{1}}, \ldots, s_{i}^{\,x_{M}}\right], b^{\rm{NL}}_{j}\right\rangle = \sum_{t = 1}^{M}b^{{\rm{NL}}}_{j}(t)\,{s_{i}^{\,x_{t}}}, &(5)}]

where [b^{{\rm{NL}}}_{j}\!(t)] is the tth element of the jth basis function [b^{{\rm{NL}}}_{j}] of [{\cal T}^{\,{\rm{NL}}}], ⊗ denotes the tensor product, and [\left[\cdot\,,\ldots,\cdot\right]] denotes the stacking of the blocks along the third dimension. The indexing and notation of the transform spectrum coefficients are illustrated in Figure 3[link].

[Figure 3]
Figure 3
Notation and indexing of patch coordinates xl, patches zxl, and coefficients si xl and [s_{i,j}^{\,x_{1},\ldots,x_{M}}] in the the corresponding [{\cal T}^{\,2\rm {D}}] spectra and [{{\cal T}^{\,3\rm {D}}}] spectrum, reproduced with permission from Mäkinen et al. (2020[Mäkinen, Y., Azzari, L. & Foi, A. (2020). IEEE Trans. Image Process. 29, 8339-8354.]). The illustration is for a group of three blocks of size 2×2 at coordinates x1 = (4,3), x2 = (7,5), x3 = (8,6) within a 10×10-pixel image.

In each step of the algorithm, the variances of [s_{i,j}^{\,{x_{1},\ldots,\,x_{M}}}] play a key role; we denote them by [v_{i,j}^{\,{x_{1},\ldots,\,x_{M}}}]. For their calculation, we refer the reader to Mäkinen et al. (2020[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Mäkinen, Y., Azzari, L. & Foi, A. (2020). IEEE Trans. Image Process. 29, 8339-8354.]. Here, we provide a summary of the macroscopic operations of the algorithm.

2.2.1. Block-matching

For each reference block, BM3D defines a local neighborhood from which similar blocks are collected. Each block in the neighborhood is ranked by

[{L_{x_{{\rm{R}}}}}{\!}\left(x_{j}\right) = \left\|z_{x_{{\rm{R}}}}\!-z_{x_{j }}\right\|_{2}^{2}\,\,-\,\,2\gamma\sum_{i = 1}^{N}v_{i,2}^{\,x_{R},\,x_{j}}, \eqno(6)]

where [z_{x_{{\rm{R}}}}] is the reference block, zxj is a potential match, vi,2 xRxj is the ith coefficient of the block-pair transform-domain variance corresponding to block difference, and [{\gamma\in{\bb{R}}}]. The common aim of block-matching is to find blocks which are the most similar to the reference block in terms of the underlying noise-free content. When only a noisy image is available, the similarity is evaluated between noisy blocks and the term scaled by γ in (6[link]) compensates for bias in the ranking caused by noise correlation. Specifically, with [{\gamma\! = \!0}] the matches would be mainly guided by the strong vertical correlation of the streak noise and thus be located along the streaks, largely ignoring any similarity of the underlying signal; setting [{\gamma\,\gt\,0}] mitigates this bias by promoting matching of blocks in which the noise is not correlated with that of the reference block. In particular, we employ [{\gamma\! = \!3}] as proposed by Mäkinen et al. (2020[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Mäkinen, Y., Azzari, L. & Foi, A. (2020). IEEE Trans. Image Process. 29, 8339-8354.] for the general case of correlated noise, facilitating further the matching of blocks which differ from the reference block mainly due to the variance of the block difference.

The common design of BM3D includes two distinct stages of denoising with different shrinkage operators, meaning that the full image is processed twice. In the second stage, the block-matching is commonly executed on the image estimate produced by the first denoising stage. As this image can be presumed noise-free, the second block-matching is executed without any compensation for noise correlation.

2.2.2. Shrinkage of the 3-D spectra

The core of BM3D is shrinkage performed on the 3-D transform spectrum of the grouped noisy blocks. For a given transform-domain coefficient of the group, a generic shrinkage can be expressed as

[s_{i,\,j}^{\,{x_{1},\ldots,x_{M}}} \,\,\longmapsto\,\,\alpha_{i,\,j}\, s_{i,\,j}^{\,{x_{1},\ldots,x_{M}}}, \eqno(7)]

where αi, j is a shrinkage attenuation factor which depends on [s_{i,\,j}^{\,{x_{1},\ldots,\,x_{M}}}], the noise statistics, and possible other priors.

BM3D utilizes two shrinkage operations: in the first denoising stage, the denoising process performs shrinkage by hard-thresholding; the second stage employs a Wiener filter, utilizing the hard-thresholding image estimate as a pilot signal.

In hard-thresholding, the shrinkage is performed by setting spectrum coefficients smaller than a threshold to zero, as they are mostly composed of noise,

[\alpha_{i,j}^{\,\rm{HT}} = \left\{ \! \matrix{ 1,\hfill & {\rm{if}}\ \left|\,s_{i,\,j}^{\,{x_{1},\ldots,x_{M}}}\right| \geq \left({{v_{i,\,j}^{\,{x_{1},\ldots,x_{M}}}}} \right)^{1/2}\lambda,\hfill \cr 0,\hfill & {\rm{otherwise,}}\hfill } \right. \eqno(8)]

where [{\lambda\,\geq\,0}] is a fixed constant.

In Wiener filtering, the attenuation coefficients of the transfer function are computed from the previous estimate, used as pilot signal, and from the variance of the noise spectrum coefficients as

[\alpha_{i,\,j}^{\rm{wie}} = {{ \big\|\left\langle\left[\hat{y}^{\rm{HT}}_{x_{1}}, \ldots, \hat{y}^{\rm{HT}}_{x_{M}}\right], b^{2{\rm{D}}}_{i} \!\otimes b^{{\rm{NL}}}_{j}\right\rangle\big\|^{2} }\over{ \Big\|\left\langle\left[\hat{y}^{\rm{HT}}_{x_{1}},\ldots,\hat{y}^{\rm{HT}}_{x_{M}}\right], b^{2{\rm{D}}}_{i} \!\otimes b^{{\rm{NL}}}_{j}\right\rangle\Big\|^{2} + \mu^{2}v_{i,\,j}^{{x_{1},\ldots,x_{M}}} }}, \eqno(9)]

where [\hat{y}^{\rm{HT}}] is the estimate of y obtained from the hard-thresholding stage, and μ2 is a scaling factor included due to aggregation to influence the bias-variance ratio we wish to minimize through the Wiener filter.

2.2.3. Aggregation

After calculating the attenuation factors of the group, they can be applied to the 3-D transform spectra to obtain estimates for the grouped blocks,

[\hat{y}_{x_{j}} = Q^{2{\rm{D}}\!}\Big\{\left\langle\alpha_{i,j}^{{x_{1},\ldots,x_{M} }}s_{i,j}^{\,{x_{1},\ldots,x_{M}}}\!, q_{j}^{{\rm{NL}}}\right\rangle\Big\},\eqno(10)]

where Q2D is the inverse transform of [{\cal T}^{\,2{\rm{D}}}], and [q_{j}^{{\rm{NL}}}] is the jth transform basis function of the inverse of [{\cal T}^{\,{\rm{NL}}}].

Hence, an estimate is produced for all blocks included in the group. As a new group is built for every position xR of the reference block, there is a large amount of block estimates providing a highly redundant covering of the image. Let XR be the set of coordinates of all reference blocks and denote by [{\{\,\hat{y}_{x_{j}}^{x_{{\rm{R}}}}\}_{j = 1}^{M_{{x}_{{\rm{R}}}}}}] the set of estimates (10[link]) for the group of blocks matched to the reference block at position xR. Then, the block estimates of an image are [\cup{}_{x_{{\rm{R}}} \in X_{{\rm{R}}}\!}{\{\,\hat{y}_{x_{j}}^{x_{ {\rm{R}}}}\}_{j = 1}^{M_{{{x_{{\rm{R}}}}}}}}] and they can all be distinct.

We aggregate all block estimates at their respective positions into the image through an adaptive weighted average,

[\hat{y} = {{ \sum_{x_{{\rm{R}}} \in X_{{\rm{R}}}} \sum_{j = 1}^{M_{{{x_{{\rm{R}}}}}}} \omega_{x_{j}}^{\,x_{{\rm{R}}}} \, W_{x_{j}}\,\hat{y}^{\,x_{{\rm{R}}}}_{x_{j}} }\over{ \sum_{x_{{\rm{R}}} \in X_{{\rm{R}}}} \sum_{j = 1}^{M_{{{x_{{\rm{R}}}}}}} \omega_{x_{j}}^{\,x_{{\rm{R}}}} \, W_{x_{j}} }}, \eqno(11)]

where [\omega_{x_{j}}^{\,x_{{\rm{R}}}}] is a block-specific weight and Wxj is a windowing function over blocks at position xj. The weights [\omega_{x_{j}}^{\,x_{{\rm{R}}}}], inversely proportional to the residuals of transform-domain noise variances, promote estimates with less residual noise to improve the quality of the final estimate.

The steps for denoising a group of blocks are demonstrated in Figure 4[link].

[Figure 4]
Figure 4
Denoising a part of Fly (a portion of Figure 2[link]) with vertical streak noise with Ψ as in Figure 1[link]. Left to right: (a) positions x1,…, x8 of one group of blocks with reference block in red; (b) contents of the eight matched blocks [{z_{x_{1}},\ldots,z_{x_{8}}}]; (c) the resulting [{{\cal T}^{\,3\rm {D}}}] spectrum coefficients [s_{i,j}^{\,{x_{1},\ldots,x_{8}}}]; (d) corresponding 3-D noise standard deviations [({v_{i,j}^{\,{x_{1},\ldots,x_{8}}}})^{1/2}]; (e) hard-thresholded coefficients [\alpha_{i,j}^{\rm{HT}}s_{i,j}^{\,{x_{1},\ldots,x_{8}}}]; (f) the denoised group of blocks [{\hat{y}^{{\rm{HT}}}_{x_{1}},\ldots,\hat{y}^{{\rm{HT}}}_{x_{8}}}], and (g) the denoising result of hard-thresholding [\hat{y}^{{\rm{HT}}}]. For the spectrum coefficients and the standard deviations (c, d, e), 50%-gray pixel color in the figure corresponds to value 0 in the data.

3. Processing pipeline

In this section, we consider the necessary steps for modeling the streak noise through (1[link]) for real sinogram data, hence allowing the effective application of BM3D for streak removal.

3.1. Bright-fielding and log-transformation

The optical attenuation through the sample is determined experimentally via bright-field corrections requiring two additional inputs: the bright-field and the dark-field (Seibert et al., 1998[Seibert, J. A., Boone, J. M. & Lindfors, K. K. (1998). Proc. SPIE, 3336, 348-354.]). The bright-field is an acquisition obtained by the imaging procedure with no sample, and the dark-field is obtained with no beam; both are 2-D arrays the size of effective pixels of the detector. Furthermore, the Beer–Lambert law relates the X-ray transform through the sample to the optical attenuation by a logarithmic transformation (Swinehart, 1962[Swinehart, D. F. (1962). J. Chem. Educ. 39, 333.]).

Hence, the raw projections Praw are first normalized as

[P_{{\rm{norm}}} = {{P_{{\rm raw}}-I_{{\rm{D}}}} \over {\,I_{{\rm B}}-I_{ {\rm{D}}}\,}}, \eqno(12)]

where ID is the dark-field and IB is the bright-field,1 and then log-transformed as

[P_{{\rm log}} = {\ln}\left(P_{{\rm{norm}}}\right). \eqno(13)]

Bright-fielding (12[link]) provides a partial, but not thorough, correction of the streak noise (Davidson et al., 2003[Davidson, D., Fröjdh, C., O'Shea, V., Nilsson, H.-E. & Rahman, M. (2003). Nucl. Instrum. Methods Phys. Res. A, 509, 146-150.]); the denoising pipeline of the following sections is designed to attenuate the remaining streak noise.

3.1.1. Noise in the projections

Apart from possible completely defective detectors2 we treat the variation in detector response as normally distributed. We further model the streak noise as locally stationary, meaning that the streak variance is presumed constant within sufficient area (i.e. the block-matching search neighborhood) for the application of BM3D. As the data are obtained through a photon-counting detector, the statistics of the measured raw data can be further modeled through a Poisson distribution. Considering both the approximately normally distributed streak noise and the Poissonian component, noise in projections normalized by (12[link]) can be modeled as

[P_{{\rm{norm}}} = A\left(1+\eta_{{}_{{\rm{P}}}}\right)+\pi = \left(A+{{\pi}\over{1+\eta_{{ }_{{\rm{P}}}}}}\right) \left(1+\eta_{{}_{{\rm{P}}}}\right), \eqno(14)]

where A are the noise-free projections, [\eta_{{}_{{\rm{P}}}}] is the normally distributed streak noise component, and π is (approximately) white Poissonian noise with zero mean; all components of (14[link]) are considered as 3-D arrays and multiplications are elementwise.

We note that the natural logarithm of (13[link]) acts as a variance-stabilizing transformation for the multiplicative noise component [(1+\eta_{{}_{{\rm{P}}}})]. Hence, we have

[P_{{\rm log}} = \ln\left(P_{{\rm{norm}}}\right) \approx {\ln}\left(A+{{\pi} \over {1\!+\!\eta_{{}_{{\rm{P}}}}}}\right)+\eta_{{}_{ {\rm{P}}}}, \eqno(15)]

where the approximation comes from [{{\ln}(1\!+\!\eta_{{}_{{\rm{P}}}}) \approx \eta_{{}_{{\rm{P}}}}}]. The additive noise component [\eta_{{}_{{\rm{P}}}}] in (15[link]) corresponds to the streak noise to be denoised. As here we only aim to attenuate the streak noise, through denoising we estimate [{{\ln}[A+{{\pi}/({1+\eta_{{}_{{\rm{P}}}}}})]}]; the embedded noise term [{{\pi}/{(1+\eta_{{}_{{\rm{P}}}})}}] although not i.i.d. is nevertheless (approximately) white and does not present streaks.

Individual sinograms, each of which is defined as a cross section of the stack of projections Plog, are denoted as

[Z = Y+\eta_{{}_{{\rm Z}}}, \eqno(16)]

where Y denotes the underlying streak-free sinogram, and [\eta_{{}_{{\rm Z}}}] is the corresponding cross section of [\eta_{{}_{{\rm{P}}}}]. The sinograms Z are used as the input for the processes in the following Section 3[link].2[link].

3.2. Multiscale filtering architecture

In the following, we assume that sinograms Z are oriented such that streaks are oriented vertically, i.e. the angular component is vertical and the displacement is horizontal.

The streak noise [\eta_{{}_{{\rm Z}}}] is characterized by very long-range correlation. In particular, because vertically there are no high-frequency streak noise components, the streaks can be filtered entirely at a coarse vertical scale, with consequent benefits in terms of efficacy and computational efficiency. Furthermore, BM3D operates using blocks of fixed size within a limited neighborhood which may be too small to fully denoise the wider streaks. Thus, we also want to denoise across multiple horizontal scales to effectively attenuate streaks of varying sizes.

Our multiscale implementation is based on a simple and efficient pixel binning to go towards coarser scales by replacing adjacent pixels by their sum. To go back towards finer scales we leverage the iterative debinning approach from Azzari & Foi (2016[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Azzari, L. & Foi, A. (2016). IEEE Signal Process. Lett. 23, 1086-1090.], which is based on spline upsampling. The multiscale denoising process is illustrated in Figure 5[link] and proceeds as follows.

[Figure 5]
Figure 5
Flowchart of the multiscale denoising process, starting from the noisy sinogram Z and resulting in the estimate [\hat{Y}] of the streak-free sinogram, both of size n×m. First, Z is vertically rescaled into a sinogram Z0 of size [{n\!\times\!m_{{\rm v}}}], [{m_{{\rm v}}\!\le m}], through the binning operator [{\cal B}_{{\rm v}}]. Then, by repeated horizontal binning [{\cal B}_{{\rm h}}], Z0 is progressively downscaled into a series of sinograms [{Z_{k}\! = \!{\cal B}_{\rm h}\left(Z_{k-1}\right)}], [{k\! = \!1,\ldots,K}], each of size [\lceil 2^{-k}n\rceil\!\times\!m_{{\rm v}}]. The coarsest scale noisy input Z *K = ZK is denoised with BM3D to produce [\hat{Y}_{K}]. Then, recursively for [{k\! = \!{K\!-\!1},\ldots,0}], the noisy input [{Z^{\,*}_{k}\! = \!Z_{k}-{\cal B}_{\rm h}^{\,-1}\!\left(Z_{k+1}\!-\!\hat{Y}_{k+ 1}\right)}] is denoised by BM3D to produce [\hat{Y}_{k}]; this definition of Z *k means that the coarse-scale horizontal components of Zk are replaced by [\hat{Y}_{k+1}]. The PSD for each scale is estimated as described in Section 3.3[link].2[link]. The resulting estimate [\hat{Y}_{0}] of the horizontal debinning (size [{n\!\times\!m_{{\rm v}}}]) similarly replaces the coarse-scale vertical components of Z to obtain the full-size estimate [\hat{Y}].

We begin with a single vertical binning of the full noisy sinogram Z of height m to a sinogram Z0 of height [m_{{\rm v}} \leq m] through a binning operator [{\cal B}_{{\rm v}}]. On Z0, we perform all consequent horizontal operations and denoising.

After vertical binning, the sinogram Z0 is progressively halved in size K times through a horizontal binning operator [{\cal B}_{{\rm h}}]: [{Z_{k}] = [{\cal B}_{{\rm h}}\left(Z_{k-1}\right)}] = [{\cal B}^{\,k}_{{\rm h}}\left(Z_{0}\right)] = [{\cal B}_{{\rm h}}^{\,k}\left[{\cal B}_{{\rm v}}\left(Z\right)\right]], [{k\! = \!1,\ldots,K}]. Denoting by n the width of Z and Z0, Zk has width ⌈2kn⌉; with every binning, the streak width also gets halved. The multiscale denoising is operated in a coarse-to-fine fashion, where progressively for each [{k = K,{K\!-\!1},\ldots,2,1,0}] we obtain an estimate [\hat{Y}_{k}] of [{{\cal B}_{{\rm h}}^{\,k}\left[{\cal B}_{{\rm v}}\left(Y\right) \right]}]. We start by taking as noisy input Z *K of BM3D the smallest binned sinogram ZK; in this way, we obtain from Z *K = ZK the coarsest estimate [\hat{Y}_{k}], which is taken as initialization for the following recursive steps executed for each scale [k\! = \!{K\!-\!1},\ldots,0]:

(1) Replace the horizontal coarser-scale components of Zk by those of the estimate [\hat{Y}_{k+1}],

[\eqalign{ Z^{\,*}_{k} & = Z_{k}-{\cal B}^{\,-1}_{{\rm h}} \left[{\cal B}_{{\rm h}}\left(Z_{k}\right)\right] + {\cal B}^{\,-1}_{{\rm h}} \left(\hat{Y}_{k+1}\right) \cr& = Z_{k}-{\cal B}^{\,-1}_{{\rm h}}\left(Z_{k+1}-\hat{Y}_{k+1} \right).}]

(2) Denoise Z *k with BM3D to produce the estimate [\hat{Y}_{k}].

The result [\hat{Y}_{0}] of the last denoising step is the fully denoised estimate the size of Z0. To produce the full-size estimate [\hat{Y}] of Y, we replace the vertical coarse-scale components of Z with those of [\hat{Y}_{0}], similar to the Step (1) above,

[\hat{Y} = Z-{\cal B}^{\,-1}_{{\rm v}} \left[{\cal B}_{{\rm v}}\left(Z _{0}\right)\right] + {\cal B}^{\,-1}_{{\rm v}}\left(\hat{Y}_{0}\right) = Z-{\cal B}^{\,-1}_{{\rm v}}\left(Z_{0}-\hat{Y}_{0}\right).]

Figure 6[link] illustrates the sinograms over the various stages of the multiscale denoising process.

[Figure 6]
Figure 6
Multiscale denoising of Fly. Left: the noisy sinogram Z. Center and right: three scales of the multiscale denoising process, each showing Zk, Z *k, and [\hat{Y}_{k}]. The full-size estimate [\hat{Y}] is displayed in Fig. 12[link].

3.3. Multiscale noise model

For BM3D denoising, we regard Z *k of each scale k as z of the model (1[link]), as

[Z^{\,*}_{k} = {\cal B}_{{\rm h}}^{\,k}\left[{\cal B}_{\rm {v}}(Y)\right]+\eta_{k}^{\,*}, \eqno(17)]

where

[\eta_{k}^{\,*} = \left\{ \matrix{ \eta_{K}\,,\hfill & k = K\,,\hfill \cr \eta_{k}-{\cal B}_{{\rm h}}^{\,-1}\left[{\cal B}_{{\rm h}}\left(\eta_{k}\right)\right]\,,\hfill & k\,\lt\,K,\hfill }\right. \eqno(18)]

and [{\eta_{k} = {\cal B}_{{\rm h}}^{\,k}\left[{\cal B}_{{\rm v}}\left(\eta _{{}_{{\rm Z}}}\right)\right]}].

This definition for [\eta_{k}^{\,*}], [{k\!<\!K}], follows from considering the coarser-scale estimate [\hat{Y}_{k}] as perfectly denoised. Similar to (3[link]), [{\eta_{k}^{\,*}}] is treated as correlated noise with PSD

[\Psi_{k}^{\,*} = {\rm var}\left({\cal{F}}\left[\eta^{\,*}_{k}\right]\right) = |X_{k}|\left|{\cal{F}}\left[g_{k}^{\,*}\right]\right|^{\,2},\qquad K\geq k \geq 0, \eqno(19)]

where gk * is a correlation kernel and [|X_{k}| = \lceil 2^{-k}n\rceil m_{{\rm v}}]. As per (4[link]), the kernel gk * can be defined as

[g_{k}^{\,*} = |X_{k}|^{-1/2}\,{\cal{F}}^{\,-1}\left[{{\rm{std}}\left({\cal{F}} \left[\eta^{\,*}_{{\rm k}}\right]\right)}\right]. \eqno(20)]

3.3.1. Multiscale PSD of white streak noise

Let [{\eta_{0}\! = \!{\cal B}_{{\rm v}}\left(\eta_{{}_{{\rm Z}}}\right)}] be horizontally white and vertically constant streak noise like in Figure 1[link]. Under this assumption for η0, we have that also all [{\eta_{k} = {\cal B}_{{\rm h}}^{\,k}\left(\eta_{0}\right)}] for [{0 \leq k\leq K}] are horizontally white and vertically constant, with variance [{{\rm var}(\eta_{k}) = 2^{k}{\rm var}(\eta_{0})}]. The doubling of the variance with every horizontal binning follows from the noise whiteness, which means that each pixel of the coarser scale sinogram is a sum of two pixels with independent noise of equal variance. Therefore, disregarding the specific support size of their actual finite realizations, we can identify these stationary random fields as

[\eta_{k} = 2^{k/2\,}{\rm{std}}(\eta_{0})\,\eta_{{\rm W}}, \eqno(21)]

where ηW is a white streak noise like in Figure 1[link] with [{{\rm var}(\eta_{{\rm W}}) = 1}]. We can hence rewrite (18[link]) as

[\eta_{k}^{\,*} = \left\{\matrix{2^{K/2\,}{\rm{std}}(\eta_{0})\, \eta_{{\rm W}}\,,\hfill & k = K,\hfill \cr 2^{k/2\,}{\rm{std}}(\eta_{0})\left(\eta_{{\rm W}}-{\cal B}_{{\rm h}}^{\,-1}\left({\cal B}_{{\rm h}}\left(\eta_{{\rm W}}\right) \right)\right)\,,\hfill & k\,\lt\,K.\hfill}\right. \eqno(22)]

This together with (20[link]) means that we can characterize [\eta^{\,*}_{k}] through a kernel gk * obtained by scaling either of two basic two kernels gc and gB,

[\eqalign{ g_{{\rm{c}}} & = |X|^{-1/2}{\cal F}^{\,-1} \left[{{\rm{std}}\left({\cal{F}}\left[\eta_{W}\right]\right)} \right], \cr g_{{\cal B}} & = |X|^{-1/2}{\cal{F}}^{\,-1} \left[{{\rm{std}}\left({\cal{F}}\left[\eta_{{\rm W}}-{\cal B}_{ {\rm h}}^{\,-1}\left({\cal B}_{{\rm h}}\left(\eta_{{\rm W}}\right) \right)\right]\right)}\right],} \eqno(23)]

by a factor

[{\varsigma_{k}\! = \!2^{k/2}{\rm{std}}(\eta_{0}) = {\rm{std}}(\eta_{k})}, \eqno(24)]

as

[g_{k}^{\,*} = \left\{\matrix{\varsigma_{k}g_{{\rm{c}}}\,,\hfill & k = K,\hfill \cr \varsigma_{k}g_{{\cal B}}\,,\hfill & k\,\lt\,K.\hfill}\right. \eqno(25)]

We note that although the equalities (23[link]) formally depend on the realization size |X|, in practice this term only renormalizes the kernel with regard to the Fourier transforms; hence (23[link]) can be computed for an arbitrary support.

The kernel gc is single-pixel wide and vertically constant like in Figure 1[link], with [{\left\|g_{{\rm{c}}}\right\|_{2} = {\rm{std}}(\eta_{{\rm W}}) = 1}]. Example noise [\eta_{k}^{\,*}], [{k\!<\!K}], and the corresponding kernel [{g_{k}^{\,*} = \varsigma_{k}g_{{\cal B}}}] and PSD [\Psi_{k}^{\,*}] (19[link]) are shown in Figure 7[link].

[Figure 7]
Figure 7
Example of noise [\eta_{k}^{\,*}] of Z *k, [{k\!<\!K}], the correlation kernel [{g_{k}^{\,*}\! = \!\varsigma_{k}g_{{\cal B}}}] where gB is produced by B and [{\cal B}^{-1\!}], and the corresponding PSD [\Psi_{k}^{\,*}]. For the kernel, 50%-gray pixel color in the figure corresponds to value 0 in the data; for the PSD, black is 0. Note the missing low frequencies at the center of the PSD, and the higher-frequency nature of [\eta_{k}^{\,*}] compared with the white streak noise η in Figure 1[link].

Estimation of ς0. To estimate [{{\rm{std}}(\eta_{0})\! = \!\varsigma_{0}}], we first convolve Z0 with a 2-D kernel [{g_{{\rm{d}}}\! = \!\phi\otimes\psi}] where ϕ is a 1-D column Gaussian function of length mv/2 and standard deviation mv/12 and ψ is a horizontal high-pass Daubechies wavelet `db3' of length 6, hence convolution with gd realizes low-pass filtering in the vertical and high-pass filtering in the horizontal. Thus, compared with Z0, [{Z_{0\,}{\bigcirc\kern-9pt\ast}\,\,g_{{\rm{d}}}}] offers a lower signal-to-noise ratio (SNR), which facilitates the estimation of noise statistics; an example of Z0 and the corresponding [{Z_{0\,}{\bigcirc\kern-9pt\ast}\,\,g_{{\rm{d}}}}] are shown in Figure 8[link] (top). One can compute an estimate of the standard deviation of [{\eta_{0\,}{\bigcirc\kern-9pt\ast}\,\,g_{{\rm{d}}}}] via its median absolute deviation (Hampel, 1974[Hampel, F. R. (1974). J. Am. Stat. Assoc. 69, 383-393.]),

[\widehat{{\rm{std}}\left(\eta_{0\,}{\bigcirc\kern-9pt\ast}\,\,g_{{\rm{d}}}\right)} = 1.4826\, \mathop{\rm{smed}}_{{X_{0}}} \left(\left|Z_{0\,}{\bigcirc\kern-9pt\ast}\,\,g_{\rm{d}} - \mathop{\rm{smed}}_{{X_{0}}} \left(Z_{0\,}{\bigcirc\kern-9pt\ast}\,\,g_{\rm{d}}\right)\right|\right), \eqno(26)]

where smed denotes the sample median and the factor 1.4826 calibrates the estimate with respect to a normal distribution of the noise. As [{{\rm{std}}\left(\eta_{0\,}{\bigcirc\kern-9pt\ast}\,\,g_{{\rm{d}}}\right)}] = [\|\varsigma_{0\,}g_{ {\rm{c}}}{\bigcirc\kern-9pt\ast}\,\,{g_{{\rm{d}}}}\|_{2}], an estimate [\hat{\varsigma}_{0}] of ς0 can be obtained through

[\hat{\varsigma}_{0} = \|g_{{\rm{c}}}{\bigcirc\kern-9pt\ast}\,\,{g_{{\rm{d}}}}\|_{2}^{-1} \, \widehat{{\rm{std}}\left(\eta_{0\,}{\bigcirc\kern-9pt\ast}\,\,g_{{\rm{d}}}\right)}. \eqno(27)]

[Figure 8]
Figure 8
A fragment of Z0 of Fly and the corresponding [{Z_{0\,}{\bigcirc\kern-8.5pt\ast}\,\,{g_{{\rm{d}}}}}] (top), and a fragment of Z *0 with the corresponding [{Z^{\,*}_{0\,}{\bigcirc\kern-8.5pt\ast}\,\,{g_{{\rm{d}}}}}] (bottom). For [{Z_{0\,}{\bigcirc\kern-8.5pt\ast}\,\,{g_{{\rm{d}}}}}] and [{Z^{\,*}_{0\,}{\bigcirc\kern-8.5pt\ast}\,\,{g_{{\rm{d}}}}}], 50%-gray pixel color in the figure corresponds to value 0 in the data. Note how most of the signal of the sinogram fragments is not present in the convolved arrays, facilitating the estimation of the streak noise statistics. Although [{Z_{0\,}{\bigcirc\kern-8.5pt\ast}\,\,{g_{{\rm{d}}}}}] and [{Z^{\,*}_{0\,}{\bigcirc\kern-8.5pt\ast}\,\,{g_{{\rm{d}}}}}] look very similar, careful visual inspection reveals slight differences similar to those between [{\eta\! = \!\eta_{k}}] in Figure 1[link] and [\eta^{\,*}_{k}] in Figure 7[link].
3.3.2. Adapting the model to non-white streak noise

The above model (21)[link]–(25)[link] assumes that the streak noise [\eta_{{}_{{\rm Z}}}] is horizontally white and stationary; however, real streak noise is never exactly white across the displacement, and may thus have significant differences in noise power between scales.

To adapt to these deviations from the model, we relax the definition (24)[link] of ςk and allow the scaling parameter [{\varsigma_{k}\!\geq\!0}] to vary with each scale k, while assuming the kernels as in (25)[link] for simplicity.3 In this way, to adaptively model the PSD (19)[link], we require only the estimation of ςk on each Z *k.

Estimation of ςk, [{k\!\geq\!0}]. For k = K, ςK can be estimated from ZK = Z *K by trivial substitutions of 0 with K in (26[link]) and (27[link]). Although also for [{k\!<\!K}] one could estimate ςk similarly from Zk, a more accurate estimate can be obtained using Z *k as this leverages the denoising of the coarser scales and thus [{Z^{\,*}_{k}{\bigcirc\kern-9pt\ast}\,\,g_{{\rm{d}}}}] offers an even lower SNR than [{Z_{k}{\bigcirc\kern-9pt\ast}\,\,g_{{\rm{d}}}}]. An example of Z *k and the corresponding [{Z^{\,*}_{k}{\bigcirc\kern-9pt\ast}\,\,g_{{\rm{d}}}}] are shown in Figure 8[link] (bottom). Similar to (26[link]) and for any [{K\!\geq\!k\!\geq\!0}], the standard deviation of the noise in [{Z^{\,*}_{k}{\bigcirc\kern-9pt\ast}\,\,{g_{{\rm{d}}}}}] can be estimated as

[\widehat{{\rm{std}}\left(\eta_{k}^{\,*\,}{\bigcirc\kern-9pt\ast}\,\,g_{{\rm{d}}}\right)} = 1.4826\,\mathop{\rm{smed}}_{{X_{k}}} \left(\left|{{Z^{\,*}_{k}{\bigcirc\kern-9pt\ast}\,\,{g_{{\rm{d}}}}}} - \mathop{\rm{smed}}_{{X_{k}}} \left(Z^{\,*}_{k}{\bigcirc\kern-9pt\ast}\,\,{g_{{\rm{d}}}}\right)\right|\right).]

Noting that [{{\rm{std}}\left(\eta_{k}^{\,*}{\bigcirc\kern-9pt\ast}\,\,g_{{\rm{d}}}\right)} = \|g^{\,*}_{k}{\bigcirc\kern-9pt\ast}\,\,{g_{{\rm{d}}}}\|_{2}], we then estimate ςk as

[\hat{\varsigma}_{k} = \left\{\matrix{\|g_{{\rm{c\,}}}{\bigcirc\kern-9pt\ast}\,\,{ g_{{\rm{d}}}}\|_{2}^{-1}\,\widehat{{\rm{std}}\left(\eta_{K}^{\,*}{\bigcirc\kern-9pt\ast}\,\,g_{ {\rm{d}}}\right)},\hfill & k = K,_{\vphantom{\big|}}\hfill \cr \|g_{{\cal B\,}}{\bigcirc\kern-9pt\ast}\,\,{g_{{\rm{d}}}}\|_{2}^{-1}\,\widehat{{\rm{std}}\left(\eta_{k}^{\,*}{\bigcirc\kern-9pt\ast}\,\,g_{{\rm{d}}}\right)},\hfill & k\,\lt\,K.\hfill}\right. \eqno(28)]

3.3.3. Horizontal nonstationarity of ηZ

Variance of the streak noise may differ across the sinogram due to changes in photon flux or noise in the bright-field. Thus, it may not be possible to denoise Z *k assuming a constant ςk for all spatial positions without either oversmoothing or leaving noise artifacts in some areas. To adapt to horizontal nonstationarity, we further relax the model allowing ςk to vary within each scale k. In particular, before noise estimation and denoising, we split Z *k and [{Z^{\,*}_{k}{\bigcirc\kern-9pt\ast}\,\,{g_{{\rm{d}}}}}] into overlapping, full-height segments. We apply BM3D separately on each segment of Z *k, using a PSD scaled by [\hat{\varsigma}_{k}] estimated on the corresponding segment of [{Z^{\,*}_{k}{\bigcirc\kern-9pt\ast}\,\,{g_{{\rm{d}}}}}], i.e. we consider each segment as a separate noisy image z with a corresponding Ψ. After denoising, the segment estimates produced by BM3D are recombined with a windowing function to form the full estimate [\hat{Y}_{k}].

3.4. Attenuation of extreme streaks

We note that the projections often include several streaks caused by defects in the scintillator. These streaks can be far stronger than what are reasonably produced by the distribution of [\eta_{{}_{{\rm{P}}}}] and therefore require a specific pre-processing. To this end, after applying the bright-field and before the multiscale denoising process, we run a simple procedure on Plog which aims to detect and attenuate only the most extreme streaks. First, we calculate the median across the angular dimension of the 3-D stack of projections as

[\tilde{P} = \mathop{\rm{smed}}_{\rm{angle}} \left(P_{\rm{log}}\right),]

resulting in a 2-D map in which the streaks present as pixels extremely brighter or darker than their surroundings. To detect extreme outliers, for each coordinate x representing a single pixel of the detector and hence of [\tilde{P}], we fit a bivariate cubic polynomial ℘x to a window [\tilde{P}_{x}] of [\tilde{P}] centered at x. Then, consistent with Gaussian modeling of [\eta_{{}_{{\rm{P}}}}], we mark the center pixel [\tilde{P}(x)] defective if [{|\tilde{P}(x)\!-\!\wp_{x}(x)|\,\,\gt\,\,4\,{\rm{sstd}}\{\tilde{P}_{x}\!-\!\wp_{x}\}}], where sstd denotes the sample standard deviation; each marked pixel in [\tilde{P}] corresponds to a full column of the sinograms.

Each pixel of a defective column is replaced with the median of non-defective pixels within a 2-D window considering the displacement dimensions around it. We note that columns corrected in this way are unlikely to be completely free of streak noise; instead, the aim is to introduce less extreme pixel values that can be further denoised by the following applications of BM3D. In order not to overload the notation, we denote the output of this step as Plog identical to its input.

The complete streak noise attenuation procedure is illustrated in Figure 9[link]. The procedure is fully automatic, requiring as an input only the raw projections and the bright- and dark-fields.

[Figure 9]
Figure 9
The complete denoising process, requiring as inputs the noisy projections Praw and the bright- and dark-fields IB, ID, and producing as the output an estimate of the streak-free stack of projections composed of sinogram estimates [\hat{Y}].

4. Experiments

We test our pipeline on synthetic data as well as two real acquisitions displaying ring artifacts. As a comparison, we show results for two leading streak-removal procedures from the Tomopy library (Gürsoy et al., 2014[Gürsoy, D., De Carlo, F., Xiao, X. & Jacobsen, C. (2014). J. Synchrotron Rad. 21, 1188-1193.]): Münch et al. (2009[Münch, B., Trtik, P., Marone, F. & Stampanoni, M. (2009). Opt. Express, 17, 8567-8591.]) and Vo et al. (2018[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Vo, N. T., Atwood, R. C. & Drakopoulos, M. (2018). Opt. Express, 26, 28396-28412.]. In particular, for the latter we combine `Algorithm 3', `Algorithm 5', and `Algorithm 6', which is demonstrated by Vo et al. (2018[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Vo, N. T., Atwood, R. C. & Drakopoulos, M. (2018). Opt. Express, 26, 28396-28412.] to attenuate a variety of different streaks.

For the synthetic experiments, we use a sinogram (627 × 180 pixels) of the Shepp–Logan phantom obtained through MATLAB Radon transform upon a sign change and an exponential transformation. We regard this sinogram as the noise-free projections A and generate noise according to (14[link]) with g as a one-pixel wide image-height vertical kernel like the one in Figure 1[link]. To obtain streak noise of different strengths, the streak noise component [(1\!+\!\eta_{{}_{{\rm{P}}}})] is generated with [{\rm{std}}(\eta_{{}_{{\rm{P}}}})] = 0.005, 0.01, 0.02, 0.05. Next, to generate noisy measurements with different SNR levels for the Poisson component, we separately scale A to the ranges [1280, 2560] (higher SNR) and [640, 1280] (lower SNR) and generate a Poisson variate with mean and variance [{A(1\!+\!\eta_{{}_{{\rm{P}}}})}], thus defining the Poissonian noise π as the difference between this Poisson variate and [{A(1\!+\!\eta_{{}_{{\rm{P}}}})}]. Furthermore, we include experiments with [{\pi = 0}] (infinite SNR), thus resulting in a total of 12 combinations of Poisson and streak noise strengths. We do not simulate extreme streaks. As the underlying data consist of only a single sinogram, we have [{Z = P_{{\rm log}}}] and we consider [{\ln}\left(A+\pi/(1\!+\!\eta_{{}_{{\rm{P}}}})\right)] as the streak-free sinogram Y. The results of the phantom experiments are collected in Table 1[link], and illustrated in Figures 10[link] and 11[link].

Table 1
Average signal-to-noise ratio (SNR) after attenuation of streaks in the Shepp–Logan phantom subject to mixed streak and Poissonian noise as in (14)[link], with different combinations of [{{\rm{std}}(\eta_{{}_{{\rm{P}}}})}] and peak values of A, with [{\rm{peak}} = \infty] being the limiting case for which [{\pi\! = \!0}]

As all of the algorithms aim to remove streak noise only, the SNR values are calculated with [Y = \ln[A+\pi/(1\!+\!\eta_{{}_{{\rm{P}}}})]] as [{{\rm SNR}(\hat{Y})] = [10\,{\rm{log}}_{10}({\rm{svar}}_{X}\{Y^{2}\}/][{\rm{smean}}_{X}((\hat{Y}\!-\!Y)^{2}))], where svar and smean denote sample variance and sample mean, respectively. Each value of the table is the average SNR over ten different noise realizations.

    SNR
Peak std(ηP) Noisy Münch et al. (2009[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Münch, B., Trtik, P., Marone, F. & Stampanoni, M. (2009). Opt. Express, 17, 8567-8591.] Vo et al. (2018[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Vo, N. T., Atwood, R. C. & Drakopoulos, M. (2018). Opt. Express, 26, 28396-28412.] Proposed
∞ (π = 0) 0.005 32.61 11.80 28.97 44.05
0.01 26.59 11.78 28.52 39.19
0.02 20.58 11.72 27.48 34.29
0.05 12.77 11.49 24.62 27.24
2560 0.005 32.66 11.85 28.32 38.41
0.01 26.64 11.82 27.89 35.90
0.02 20.63 11.77 26.90 32.63
0.05 12.82 11.54 24.21 26.67
1280 0.005 32.71 11.90 27.76 36.51
0.01 26.69 11.87 27.36 34.31
0.02 20.68 11.82 26.45 31.55
0.05 12.86 11.59 23.92 26.21
[Figure 10]
Figure 10
Left: comparison of sinograms after different denoising procedures on the Shepp–Logan phantom with noise as in (14)[link] with [{\rm{std}}(\eta_{{}_{{\rm{P}}}})] = 0.02 and signal peak 2560. Top to bottom: [Y = \ln[A+\pi/(1+\eta_{{}_{{\rm{P}}}})]]; Z; Münch et al. (2009[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Münch, B., Trtik, P., Marone, F. & Stampanoni, M. (2009). Opt. Express, 17, 8567-8591.]; Vo et al. (2018[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Vo, N. T., Atwood, R. C. & Drakopoulos, M. (2018). Opt. Express, 26, 28396-28412.]; proposed procedure based on BM3D denoising. Right: corresponding estimation errors. Note how both of the comparison methods create strong artifacts around the areas with the highest contrast, as pointed by the arrows. These artifacts are not present in either the noisy image Z or the BM3D-based estimate.
[Figure 11]
Figure 11
Corresponding tomograms of Figure 10[link]. Top: [{Y\! = \!{\ln}[A+\pi/(1+\eta_{{}_{{\rm{P}}}})]}] and Z. Bottom, left-to-right: Münch et al. (2009[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Münch, B., Trtik, P., Marone, F. & Stampanoni, M. (2009). Opt. Express, 17, 8567-8591.], Vo et al. (2018[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Vo, N. T., Atwood, R. C. & Drakopoulos, M. (2018). Opt. Express, 26, 28396-28412.], and proposed procedure based on BM3D denoising. Note the strong circular components on both Münch et al. (2009[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Münch, B., Trtik, P., Marone, F. & Stampanoni, M. (2009). Opt. Express, 17, 8567-8591.] and Vo et al. (2018[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Vo, N. T., Atwood, R. C. & Drakopoulos, M. (2018). Opt. Express, 26, 28396-28412.] which are method artifacts present only in the results by these two algorithms.

The Fly dataset consists of 180 projections with 50 s exposure (detector pixel size 27 µm, demagnified to 15.7 µm by cone-beam geometry) collected using a Sigray Prisma X-ray micro-tomography instrument at 34 kV; each sinogram is 512 pixels wide. The Fly contains both extreme streaks caused by defective detectors and approximately normally distributed streaks, although they are generally more intense towards the edges of the projections due to weak photon flux and Poisson noise affecting the bright-field. Thus, the Fly benefits greatly from both the extreme streak removal procedure of Section 3[link].4[link] and relaxing the stationarity assumption by performing PSD estimation and denoising in multiple parts for each sinogram as described in Section 3.3.3[link]. The denoising results of two different sinograms are shown in Figure 12[link]; the corresponding tomograms of the second sinogram of Figure 12[link] are shown in Figure 13[link].

[Figure 12]
Figure 12
Comparison of two sinograms of Fly after different denoising procedures. Top to bottom: noisy sinogram Z; Münch et al. (2009[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Münch, B., Trtik, P., Marone, F. & Stampanoni, M. (2009). Opt. Express, 17, 8567-8591.]; Vo et al. (2018[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Vo, N. T., Atwood, R. C. & Drakopoulos, M. (2018). Opt. Express, 26, 28396-28412.]; proposed procedure based on BM3D denoising. Although Vo et al. (2018[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Vo, N. T., Atwood, R. C. & Drakopoulos, M. (2018). Opt. Express, 26, 28396-28412.] is very effective at removing streaks, it also considerably affects the sinogram features; note, for example, the considerably weaker bold diagonal line (indicated by the first arrow) compared with the other algorithms. Both Münch et al. (2009[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Münch, B., Trtik, P., Marone, F. & Stampanoni, M. (2009). Opt. Express, 17, 8567-8591.] and Vo et al. (2018[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Vo, N. T., Atwood, R. C. & Drakopoulos, M. (2018). Opt. Express, 26, 28396-28412.] also distort larger areas of the sinograms, as pointed out by the second arrow; these problems are absent from the BM3D-based result. Although not visually obvious here, the differences cause severe artifacts in the tomograms, as can be seen in Figure 13[link].
[Figure 13]
Figure 13
Comparison of resulting tomograms after different denoising procedures on the second sinogram of Fly shown in Figure 12[link]. Top to bottom: noisy reconstructed object; Münch et al. (2009[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Münch, B., Trtik, P., Marone, F. & Stampanoni, M. (2009). Opt. Express, 17, 8567-8591.]; Vo et al. (2018[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Vo, N. T., Atwood, R. C. & Drakopoulos, M. (2018). Opt. Express, 26, 28396-28412.]; proposed procedure based on BM3D denoising. Although all methods achieve good results in removing the streaks, both Münch et al. (2009[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Münch, B., Trtik, P., Marone, F. & Stampanoni, M. (2009). Opt. Express, 17, 8567-8591.] and Vo et al. (2018[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Vo, N. T., Atwood, R. C. & Drakopoulos, M. (2018). Opt. Express, 26, 28396-28412.] introduce strong shadows absent from the proposed estimate, as indicated by the arrows.

We also test the algorithm on a soft tissue sample 00076 displaying severe ring artifacts freely available in TomoBank (De Carlo et al., 2018[De Carlo, F., Gürsoy, D., Ching, D. J., Batenburg, K. J., Ludwig, W., Mancini, L., Marone, F., Mokso, R., Pelt, D. M., Sijbers, J. & Rivers, M. (2018). Meas. Sci. Technol. 29, 034004.]). The data contain 2000[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.] projections with 2.2 µm pixels, 100 ms exposure time obtained at the Advanced Photon Source, 2-BM beamline; other experimental parameters are an X-ray energy of 60–70 keV, 10 µm LuAG scintillator, and sample-to-detector distance as 90 mm. The sinograms are 2560 pixels in width. Included are ten samples for bright- and dark-fields, which are averaged to obtain a single bright-field and dark-field. The denoising results of a single sinogram are shown in Figure 14[link], and the corresponding tomograms in Figure 15[link].

[Figure 14]
Figure 14
Comparison of sinograms after different denoising procedures on 00076. Top to bottom: noisy sinogram Z; Münch et al. (2009[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Münch, B., Trtik, P., Marone, F. & Stampanoni, M. (2009). Opt. Express, 17, 8567-8591.]; Vo et al. (2018[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Vo, N. T., Atwood, R. C. & Drakopoulos, M. (2018). Opt. Express, 26, 28396-28412.]; proposed procedure based on BM3D denoising. Note how in the zoom-in the proposed method manages to remove streak noise without creating additional artifacts. Münch et al. (2009[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Münch, B., Trtik, P., Marone, F. & Stampanoni, M. (2009). Opt. Express, 17, 8567-8591.] creates a horizontal streak-like artifacts as seen in the middle of the zoom-in, not present in the noisy sinogram; Vo et al. (2018[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Vo, N. T., Atwood, R. C. & Drakopoulos, M. (2018). Opt. Express, 26, 28396-28412.] does not fully denoise the sinogram.
[Figure 15]
Figure 15
Comparison of resulting tomograms after different denoising procedures on 00076. Top to bottom: noisy reconstructed object; Münch et al. (2009[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Münch, B., Trtik, P., Marone, F. & Stampanoni, M. (2009). Opt. Express, 17, 8567-8591.]; Vo et al. (2018[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Vo, N. T., Atwood, R. C. & Drakopoulos, M. (2018). Opt. Express, 26, 28396-28412.]; proposed procedure based on BM3D denoising. Münch et al. (2009[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Münch, B., Trtik, P., Marone, F. & Stampanoni, M. (2009). Opt. Express, 17, 8567-8591.] manage to remove almost all noise in both low and high frequencies, but create artifacts where the original did not have any, as seen from the rightmost zoom. The proposed denoising procedure removes most of the noise [including the wide streaks still present in Vo et al. (2018[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Vo, N. T., Atwood, R. C. & Drakopoulos, M. (2018). Opt. Express, 26, 28396-28412.]], and does not introduce further artifacts. Note also the central pixel, magnified in the middle, which is very dark for both Münch et al. (2009[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Münch, B., Trtik, P., Marone, F. & Stampanoni, M. (2009). Opt. Express, 17, 8567-8591.] and Vo et al. (2018[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Vo, N. T., Atwood, R. C. & Drakopoulos, M. (2018). Opt. Express, 26, 28396-28412.], whereas the proposed procedure does not leave any visible artifact.

4.1. Implementation details

We use the BM3D implementation for Python (available from the PyPI package bm3d) with the `refilter' profile and input PSD estimated as described in Section 3.3[link].2[link].

For the multiscale denoising procedure, we performed vertical binning with [m_{{\rm v}}] = [\lceil m/\lceil m/64\rceil\rceil] ≃ 64 pixels; this value, being slightly larger than the height of the BM3D search neighborhood of Section 2.2[link].1[link] (39 × 39 pixels), allows our method to attenuate streaks which change slowly across the angle – a larger value of mv might be used to deal with streaks featuring faster variation. The number of horizontal scales K for each sinogram was set as [{K\! = \!\lfloor{\rm log}_{2}(n/40)\rfloor}], which gives K = 3 for Fly and the Shepp–Logan phantom, and K = 6 for 00076; these values offer a compromise between denoising wide streaks versus preserving low-frequency signal components – larger values of K not only result in ZK narrower than the BM3D search neighborhood but also in extremely coarse scales that naturally feature a very high SNR that may lead to overestimating ςK and hence to oversmoothing. For the localized processing of each scale (Section 3.3[link].3[link]), we estimate ςk, [{K \geq k\geq0}] and apply BM3D denoising on 39-pixel wide segments, following the width of BM3D search neighborhood. For the attenuation of extreme streaks, we used a 19 × 19 pixels window.

To consider the computational cost, we note that the full denoising process of Fly (181 × 512 × 512 pixels) run single-threaded on an AMD Ryzen 7 1700 processor takes about one hour, mostly due to the BM3D denoising in CPU. A highly parallel GPU-based implementation is expected to reduce this run time to the scale of seconds (Davy & Ehret, 2020[Davy, A. & Ehret, T. (2020). J. Real-Time Image Process. https://doi.org/10.1007/s11554-020-00945-4.]). The complexity of BM3D is linear with the number of pixels in the noisy image. Thus, the computation time is directly proportional to the sinogram height after vertical binning. As each iteration of the horizontal multiscale denoising halves the number of pixels, the computational cost of BM3D for an extra iteration k > 0 is then a 1/2kth of the cost of k = 0, the total for all [{k\! = \!0,\ldots,K}] being at most twice that of single scale denoising of Z0.

The correlation kernels gc and gB do not depend on the input or scale, and can thus be pre-computed. To compute gB, we use directly the definition (23[link]) through a Monte Carlo simulation of sample standard deviation in the Fourier domain. We note that, as the kernel is vertically constant, it suffices to perform this simulation in 1-D and repeat the kernel mv times in the vertical dimension.

For Münch et al. (2009[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Münch, B., Trtik, P., Marone, F. & Stampanoni, M. (2009). Opt. Express, 17, 8567-8591.] and Vo et al. (2018[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Vo, N. T., Atwood, R. C. & Drakopoulos, M. (2018). Opt. Express, 26, 28396-28412.], we use implementations remove_stripe_fw and remove_all_stripe provided by the tomopy Python library of Gürsoy et al. (2014[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Gürsoy, D., De Carlo, F., Xiao, X. & Jacobsen, C. (2014). J. Synchrotron Rad. 21, 1188-1193.]. The tomograms of each experiment are reconstructed using the xpack library (Marchesini et al., 2020[Marchesini, S., Trivedi, A., Enfedaque, P., Perciano, T. & Parkinson, D. (2020). Proceedings of the 20th International Conference on Computational Science (ICCS 2020), edited by V. V. Krzhizhanovskaya, G. Závodszky, M. H. Lees, J. J. Dongarra, P. M. A. Sloot, S. Brissos and J. Teixeira, 3-5 June 2020, Amsterdam, The Netherlands, pp. 248-261. Cham: Springer International Publishing., ]).

5. Discussion and conclusions

We have presented a model for streak noise in the sinogram domain as locally stationary correlated noise additive in the logarithmic scale. Based on this model, we have described a BM3D-based multiscale denoising procedure removing streak noise, and, consequently, the tomogram ring artifacts. The use of the recently proposed variant (Mäkinen et al., 2020[Mäkinen, Y., Azzari, L. & Foi, A. (2020). IEEE Trans. Image Process. 29, 8339-8354.]) of BM3D is crucial for this work, as we deal with long-range noise correlation which earlier BM3D designs could not handle satisfactorily.

Tested on both synthetic and real data, our denoising procedure achieves state-of-the-art performance in streak removal. Compared with the two popular streak removal algorithms (Münch et al., 2009[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.][Münch, B., Trtik, P., Marone, F. & Stampanoni, M. (2009). Opt. Express, 17, 8567-8591.]; Vo et al., 2018[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Vo, N. T., Atwood, R. C. & Drakopoulos, M. (2018). Opt. Express, 26, 28396-28412.], our procedure achieves superior results both visually and quantitatively in terms of signal-to-noise ratio. Although all tested algorithms manage to successfully remove most streak noise, both Münch et al. (2009[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Münch, B., Trtik, P., Marone, F. & Stampanoni, M. (2009). Opt. Express, 17, 8567-8591.] and Vo et al. (2018[Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911-6930.])[Vo, N. T., Atwood, R. C. & Drakopoulos, M. (2018). Opt. Express, 26, 28396-28412.] tend to create large distortions especially where the intensities of the underlying sinogram columns vary significantly. This type of artifact hinders interpretation of the results and subsequent analysis such as segmentation. In comparison, the proposed algorithm offers similar or better streak removal without further distorting the sinogram.

The proposed multiscale framework, here using basic pixel binning and debinning, could be extended to more sophisticated scale decompositions, such as steerable pyramids, contourlets, or dual-tree wavelets (Kovacevic & Chebira, 2008[Kovacevic, J. & Chebira, A. (2008). An Introduction to Frames. Now Publishers Inc.]). In this work, we have considered each sinogram as a separate input for BM3D for simplicity, but the same mechanism can be used for simultaneous filtering of the entire 3-D stack of sinograms. In this way, the similarities between consecutive sinograms of the stack could be utilized within the collaborative filter.

As the sinograms are commonly processed after a logarithmic transform, we have not discussed inversion of the logarithm needed for denoising. Although the exponential function is naturally the inverse of the logarithm, the non­linearity of the logarithm causes bias in a denoised sinogram if inverted this way. Hence, if the sinogram should be reverted back from the logarithmic domain after denoising, an exact unbiased inverse (Mäkitalo & Foi, 2010[Mäkitalo, M. & Foi, A. (2011). IEEE Trans. Image Process. 20, 99-109.]; Mäkitalo et al., 2010[Mäkitalo, M., Foi, A., Fevralev, D. & Lukin, V. (2010). Proceedings of the 2010 International Conference on Mathematical Methods in Electromagnetic Theory (MMET 2010), 6-8 September 2010, Kyiv, Ukraine, pp. 53-56.]) should be used instead.

Footnotes

1As Praw is a 3-D array, the pixels of IB and ID are replicated through the angle dimension for the operations in (12[link]).

2Extreme streak noise arising from defective detectors is addressed separately in Section 3[link].4[link].

3Adopting the kernels (25[link]) with arbitrary values of ςk corresponds to assuming ηk approximately white with variance [\varsigma_{k}^{2}], which may differ from 2kstd(η0) against (21[link]). This assumption becomes increasingly appropriate as k grows for non-white η0 featuring mild local horizontal correlations, as the binning in [{\cal B}^{\,k}_{{\rm h}}(\eta_{0})] is tantamount to a convolution and decimation, leading to a flattening of the PSD.

Acknowledgements

We are very thankful to Dr Sheraz Gul of Sigray Inc. for providing us with the experimental Fly data.

Funding information

Funding for this research was provided by: Academy of Finland (grant No. 310779); Sigray Inc. partial support by Hong Kong RGC 14302920.

References

First citationAbu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). Phys. Med. Biol. 55, 6911–6930.  PubMed Google Scholar
First citationArtul, S. (2013). BMJ Case Rep. 2013, bcr-2013–201379.  Google Scholar
First citationAzzari, L. & Foi, A. (2016). IEEE Signal Process. Lett. 23, 1086–1090.  CrossRef Google Scholar
First citationBoas, F. E. & Fleischmann, D. (2012). Imaging Med. 4, 229–240.  CrossRef Google Scholar
First citationCroton, L. C., Ruben, G., Morgan, K. S., Paganin, D. M. & Kitchen, M. J. (2019). Opt. Express, 27, 14231–14245.  CrossRef CAS PubMed Google Scholar
First citationDabov, K., Foi, A., Katkovnik, V. & Egiazarian, K. (2007). IEEE Trans. Image Process. 16, 2080–2095.  Web of Science CrossRef PubMed Google Scholar
First citationDabov, K., Foi, A., Katkovnik, V. & Egiazarian, K. O. (2008). Proc. SPIE, 6812, 681207.  CrossRef Google Scholar
First citationDavidson, D., Fröjdh, C., O'Shea, V., Nilsson, H.-E. & Rahman, M. (2003). Nucl. Instrum. Methods Phys. Res. A, 509, 146–150.  CrossRef CAS Google Scholar
First citationDavy, A. & Ehret, T. (2020). J. Real-Time Image Process. https://doi.org/10.1007/s11554-020-00945-4Google Scholar
First citationDe Carlo, F., Gürsoy, D., Ching, D. J., Batenburg, K. J., Ludwig, W., Mancini, L., Marone, F., Mokso, R., Pelt, D. M., Sijbers, J. & Rivers, M. (2018). Meas. Sci. Technol. 29, 034004.  Web of Science CrossRef Google Scholar
First citationGürsoy, D., De Carlo, F., Xiao, X. & Jacobsen, C. (2014). J. Synchrotron Rad. 21, 1188–1193.  Web of Science CrossRef IUCr Journals Google Scholar
First citationHaibel, A. (2008). Advanced Tomographic Methods in Materials Research and Engineering, Monographs on the Physics and Chemistry of Materials, pp. 141–160.  Google Scholar
First citationHampel, F. R. (1974). J. Am. Stat. Assoc. 69, 383–393.  CrossRef Google Scholar
First citationJha, A. K., Purandare, N. C., Shah, S., Agrawal, A., Puranik, A. D. & Rangarajan, V. (2013). Indian J. Nucl. Med. 28, 232–233.  PubMed Google Scholar
First citationKovacevic, J. & Chebira, A. (2008). An Introduction to Frames. Now Publishers Inc.  Google Scholar
First citationMäkinen, Y., Azzari, L. & Foi, A. (2020). IEEE Trans. Image Process. 29, 8339–8354.  Google Scholar
First citationMäkitalo, M. & Foi, A. (2011). IEEE Trans. Image Process. 20, 99–109.  PubMed Google Scholar
First citationMäkitalo, M., Foi, A., Fevralev, D. & Lukin, V. (2010). Proceedings of the 2010 International Conference on Mathematical Methods in Electromagnetic Theory (MMET 2010), 6–8 September 2010, Kyiv, Ukraine, pp. 53–56.  Google Scholar
First citationMarchesini, S., Trivedi, A., Enfedaque, P., Perciano, T. & Parkinson, D. (2020). Proceedings of the 20th International Conference on Computational Science (ICCS 2020), edited by V. V. Krzhizhanovskaya, G. Závodszky, M. H. Lees, J. J. Dongarra, P. M. A. Sloot, S. Brissos and J. Teixeira, 3–5 June 2020, Amsterdam, The Netherlands, pp. 248–261. Cham: Springer International Publishing.,   Google Scholar
First citationMünch, B., Trtik, P., Marone, F. & Stampanoni, M. (2009). Opt. Express, 17, 8567–8591.  Web of Science PubMed Google Scholar
First citationPaleo, P. & Mirone, A. (2015). J. Synchrotron Rad. 22, 1268–1278.  Web of Science CrossRef IUCr Journals Google Scholar
First citationPelt, D. M. & Parkinson, D. Y. (2018). Meas. Sci. Technol. 29, 034002.  Web of Science CrossRef Google Scholar
First citationSeibert, J. A., Boone, J. M. & Lindfors, K. K. (1998). Proc. SPIE, 3336, 348–354.  CrossRef Google Scholar
First citationSwinehart, D. F. (1962). J. Chem. Educ. 39, 333.  CrossRef Google Scholar
First citationVidal, F. P., Létang, J. M., Peix, G. & Cloetens, P. (2005). Nucl. Instrum. Methods Phys. Res. B, 234, 333–348.  Web of Science CrossRef CAS Google Scholar
First citationVo, N. T., Atwood, R. C. & Drakopoulos, M. (2018). Opt. Express, 26, 28396–28412.  CrossRef PubMed Google Scholar

This is an open-access article distributed under the terms of the Creative Commons Attribution (CC-BY) Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are cited.

Journal logoJOURNAL OF
SYNCHROTRON
RADIATION
ISSN: 1600-5775
Follow J. Synchrotron Rad.
Sign up for e-alerts
Follow J. Synchrotron Rad. on Twitter
Follow us on facebook
Sign up for RSS feeds