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

Here, streak noise is modeled in the sinogram domain as correlated noise, which is attenuated by the Block-Matching and 3-D (BM3D) filtering algorithm applied across multiple scales. The proposed procedure is fully automatic and attenuates streak noise and the corresponding ring artifacts without creating major distortions common to other streak removal algorithms.


Introduction
Ring artifacts are ubiquitous in computed tomography (Jha et al., 2013;Artul, 2013;Boas & Fleischmann, 2012); they originate from angular streak noise in measured raw sinogram data used to reconstruct a tomographic volume (Croton et al., 2019) 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;Vidal et al., 2005;Anas et al., 2010).
Minimization of ring artifacts by using adequate scanning protocols (Pelt & Parkinson, 2018), 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), combinations of polynomial smoothing filters and careful calibration of the detector response function (Vo et al., 2018;Croton et al., 2019), or iterative algorithms (Paleo & Mirone, 2015) 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 et al., , 2008, leveraging the recent inclusion of exact transformdomain noise variances (Mä kinen et al., 2020), which allow for accurate modeling of long noise correlation within the jointly transformed blocks. ISSN 1600-5775 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).

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 for the more general case of real-world streak noise.

Correlated noise model
We consider a noisy input z : X ! R to be a combination of underlying data y and additive stationary spatially correlated noise to be filtered, zðxÞ ¼ yðxÞ þ ðxÞ; x 2 X; where x 2 X & Z 2 is the coordinate in the finite two-dimensional image domain X (representing angles and displacements when z is a sinogram) and with being the zero-mean independent and identically distributed (i.i.d.) Gaussian noise with unit variance, and ' Ã ' 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 kgk 2 ¼ stdðÞ. An equivalent way of representing correlated noise is by its power spectral density (PSD) É, 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)-(3) can be defined from É as 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). 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), when setting g as a vertically constant line. This simple kernel as well as the corresponding PSD É are demonstrated in Figure 1; an example of real streak noise viably approximated through this model is shown in Figure 2.
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, a complete processing pipeline is further proposed to allow modeling more complex cases of streak noise through (1)-(4), enabling their attenuation through the collaborative filter.

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) of BM3D for correlated noise denoising where the Example noise ¼ Ã g, Á ð Þ $ N 0; 1 ð Þ, 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
Example of streak noise in a sinogram (a fragment of Fly) that can be seen as well approximated by the model in Figure 1. 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 blockmatching.
(ii) Obtain the 3-D transform spectrum by collectively transforming the obtained blocks.
(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 T 2D locally to each block, then a 1-D transform T NL through the 'stack' of grouped blocks. Denoting by z x 1 ; . . . ; z x M È É a group of M blocks of N pixels extracted from z at coordinates x 1 , . . . , x M , we obtain the T 2D spectrum coefficients as s i is the ith basis function of T 2D . The 3-D spectrum coefficients are calculated through the direct tensor product of the T 2D and T NL transforms, as where b NL j ðtÞ is the tth element of the jth basis function b NL j of T NL , denotes the tensor product, and Á ; . . . ; Á ½ denotes the stacking of the blocks along the third dimension. The indexing and notation of the transform spectrum coefficients are illustrated in Figure 3.
In each step of the algorithm, the variances of s x 1 ;...; x M i;j play a key role; we denote them by v x 1 ;...; x M i;j . For their calculation, we refer the reader to Mä kinen et al. (2020). 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 where z x R is the reference block, z x j is a potential match, v is the ith coefficient of the block-pair transform-domain variance corresponding to block difference, and 2 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) compensates for bias in the ranking caused by noise correlation. Specifically, with ¼ 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 > 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 ¼ 3 as proposed by Mä kinen et al.
(2020) 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 where i, j is a shrinkage attenuation factor which depends on s x 1 ;...; x M i; j , 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, where ! 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 whereŷ y HT is the estimate of y obtained from the hardthresholding stage, and 2 is a scaling factor included due to Notation and indexing of patch coordinates x l , patches z x l , and coefficients s x l i and s x 1 ;...;x M i;j in the the corresponding T 2D spectra and T 3D spectrum, reproduced with permission from Mä kinen et al. (2020). The illustration is for a group of three blocks of size 2Â2 at coordinates 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, where Q 2D is the inverse transform of T 2D , and q NL j is the jth transform basis function of the inverse of T NL .
Hence, an estimate is produced for all blocks included in the group. As a new group is built for every position x R of the reference block, there is a large amount of block estimates providing a highly redundant covering of the image. Let X R be the set of coordinates of all reference blocks and denote by fŷ y j¼1 the set of estimates (10) for the group of blocks matched to the reference block at position x R . Then, the block estimates of an image are [ j¼1 and they can all be distinct.
We aggregate all block estimates at their respective positions into the image through an adaptive weighted average, where !
x R x j is a block-specific weight and W x j is a windowing function over blocks at position x j . The weights !
x R x j , 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.

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

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). 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).
Hence, the raw projections P raw are first normalized as where I D is the dark-field and I B is the bright-field, 1 and then log-transformed as Bright-fielding (12) provides a partial, but not thorough, correction of the streak noise (Davidson et al., 2003); 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 detectors 2 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) can be modeled as where A are the noise-free projections, P is the normally distributed streak noise component, and is (approximately) white Poissonian noise with zero mean; all components of (14) are considered as 3-D arrays and multiplications are elementwise.
We note that the natural logarithm of (13) acts as a variance-stabilizing transformation for the multiplicative noise component ; (f ) the denoised group of blockŝ y y HT x 1 ; . . . ;ŷ y HT x 8 , and (g) the denoising result of hard-thresholdingŷ y 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.
where the approximation comes from lnð1þ P Þ % P . The additive noise component P in (15) corresponds to the streak noise to be denoised. As here we only aim to attenuate the streak noise, through denoising we estimate ln½A þ =ð1 þ P Þ; the embedded noise term =ð1 þ 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 P log , are denoted as where Y denotes the underlying streak-free sinogram, and Z is the corresponding cross section of P . The sinograms Z are used as the input for the processes in the following Section 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 Z is characterized by very long-range correlation. In particular, because vertically there are no highfrequency 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), which is based on spline upsampling. The multiscale denoising process is illustrated in Figure 5 and proceeds as follows.
We begin with a single vertical binning of the full noisy sinogram Z of height m to a sinogram Z 0 of height m v m through a binning operator B v . On Z 0 , we perform all consequent horizontal operations and denoising.
After vertical binning, the sinogram Z 0 is progressively halved in size K times through a horizontal binning operator Denoting by n the width of Z and Z 0 , Z k has width d2 Àk ne; 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; . . . ; 2; 1; 0 we obtain an estimateŶ We start by taking as noisy input Z Ã K of BM3D the smallest binned sinogram Z K ; in this way, we obtain from Z Ã K ¼ Z K the coarsest estimateŶ Y k , which is taken as initialization for the following recursive steps executed for each scale k ¼ KÀ1; . . .  Flowchart of the multiscale denoising process, starting from the noisy sinogram Z and resulting in the estimateŶ Y of the streak-free sinogram, both of size nÂm. First, Z is vertically rescaled into a sinogram Z 0 of size nÂm v , m v m , through the binning operator B v . Then, by repeated horizontal binning B h , Z 0 is progressively downscaled into a series of sinograms Z k ¼ B h Z kÀ1 ð Þ, k ¼ 1; . . . ; K, each of size d2 Àk neÂm v . The coarsest scale noisy input Z Ã K ¼ Z K is denoised with BM3D to producê Y Y K . Then, recursively for k ¼ KÀ1; . . . ; 0, the noisy input is denoised by BM3D to produceŶ Y k ; this definition of Z Ã k means that the coarse-scale horizontal components of Z k are replaced byŶ Y kþ1 . The PSD for each scale is estimated as described in Section 3.3.2. The resulting estimateŶ Y 0 of the horizontal debinning (size nÂm v ) similarly replaces the coarse-scale vertical components of Z to obtain the full-size estimateŶ Y.
(1) Replace the horizontal coarser-scale components of Z k by those of the estimateŶ Y kþ1 , (2) Denoise Z Ã k with BM3D to produce the estimateŶ Y k . The resultŶ Y 0 of the last denoising step is the fully denoised estimate the size of Z 0 . To produce the full-size estimateŶ Y of Y, we replace the vertical coarse-scale components of Z with those ofŶ Y 0 , similar to the Step (1) above, Figure 6 illustrates the sinograms over the various stages of the multiscale denoising process.

Multiscale noise model
For BM3D denoising, we regard Z Ã k of each scale k as z of the model (1), as Ã . This definition for Ã k , k < K, follows from considering the coarser-scale estimateŶ Y k as perfectly denoised. Similar to (3), Ã k is treated as correlated noise with PSD where g Ã k is a correlation kernel and jX k j ¼ d2 Àk nem v . As per (4), the kernel g Ã k can be defined as 3.3.1. Multiscale PSD of white streak noise.
be horizontally white and vertically constant streak noise like in Figure 1. Under this assumption for 0 , we have that also all k ¼ B k h 0 ð Þ for 0 k K are horizontally white and vertically constant, with variance varð k Þ ¼ 2 k varð 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 where W is a white streak noise like in Figure 1 with varð W Þ ¼ 1. We can hence rewrite (18) as This together with (20) means that we can characterize Ã k through a kernel g Ã k obtained by scaling either of two basic two kernels g c and g B , by a factor We note that although the equalities (23) formally depend on the realization size |X|, in practice this term only renormalizes the kernel with regard to the Fourier transforms; hence (23) can be computed for an arbitrary support. The kernel g c is single-pixel wide and vertically constant like in Figure 1 Figure 7. Estimation of & 0 . To estimate stdð 0 Þ¼ & 0 , we first convolve Z 0 with a 2-D kernel g d ¼ where is a 1-D column Gaussian function of length m v /2 and standard deviation m v /12 and is a horizontal high-pass Daubechies wavelet 'db3' of length 6, hence convolution with g d realizes low-pass filtering in the vertical and high-pass filtering in the horizontal. Thus, compared with Z 0 , Z 0 Ã g d offers a lower signal-to-noise ratio (SNR), which facilitates the estimation of noise statistics; an example of Z 0 and the corresponding Z 0 Ã g d are shown in Figure 8 (top). One can compute an estimate of the standard deviation of 0 Ã g d via its median absolute deviation (Hampel, 1974), Multiscale denoising of Fly. Left: the noisy sinogram Z. Center and right: three scales of the multiscale denoising process, each showing Z k , Z Ã k , andŶ Y k . The full-size estimateŶ Y is displayed in Fig. 12.
where smed denotes the sample median and the factor 1.4826 calibrates the estimate with respect to a normal distribution of the noise. As std 0 3.3.2. Adapting the model to non-white streak noise. The above model (21)-(25) assumes that the streak noise 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) of & k and allow the scaling parameter & k ! 0 to vary with each scale k, while assuming the kernels as in (25) for simplicity. 3 In this way, to adaptively model the PSD (19), we require only the estimation of & k on each Z Ã k . (26) and (27). Although also for k < K one could estimate & k similarly from Z k , a more accurate estimate can be obtained using Z Ã k as this leverages the denoising of the coarser scales and thus Z Ã k Ã g d offers an even lower SNR than Z k Ã g d . An example of Z Ã k and the corresponding Z Ã k Ã g d are shown in Figure 8 (bottom). Similar to (26) and for any K ! k ! 0, the standard deviation of the noise in Z Ã k Ã g d can be estimated as Noting that std Ã k Ã g d 3.3.3. Horizontal nonstationarity of g 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 Ã g d into overlapping, full-height segments. We apply BM3D separately on each segment of Z Ã k , using a PSD scaled by& & k estimated on the corresponding segment of Z Ã k Ã g 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Ŷ Y k .

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 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 P log 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 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 ofP P, we fit a bivariate A fragment of Z 0 of Fly and the corresponding Z 0 Ã g d (top), and a fragment of Z Ã 0 with the corresponding Z Ã 0 Ã g d (bottom). For Z 0 Ã g d and Z Ã 0 Ã g 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 Ã g d and Z Ã 0 Ã g d look very similar, careful visual inspection reveals slight differences similar to those between ¼ k in Figure 1 and Ã k in Figure 7.

Figure 7
Example of noise Ã k of Z Ã k , k < K, the correlation kernel g Ã k ¼ & k g B where g B is produced by B and B À1 , and the corresponding PSD É Ã 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 Ã k compared with the white streak noise in Figure 1. cubic polynomial } x to a windowP P x ofP P centered at x. Then, consistent with Gaussian modeling of P , we mark the center pixelP PðxÞ defective if jP PðxÞÀ} x ðxÞj > 4 sstdfP P x À} x g, where sstd denotes the sample standard deviation; each marked pixel inP 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 P log identical to its input.
The complete streak noise attenuation procedure is illustrated in Figure 9. The procedure is fully automatic, requiring as an input only the raw projections and the bright-and dark-fields.

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 Vo et al. (2018). In particular, for the latter we combine 'Algorithm 3', 'Algorithm 5', and 'Algorithm 6', which is demonstrated by Vo et al. (2018) 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) with g as a one-pixel wide image-height vertical kernel like the one in Figure 1. To obtain streak noise of different strengths, the streak noise component ð1þ P Þ is generated with stdð 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þ P Þ, thus defining the Poissonian noise as the difference between this Poisson variate and Að1þ P Þ. Furthermore, we include experiments with ¼ 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 log and we consider ln A þ =ð1þ P Þ À Á as the streak-free sinogram Y. The results of the phantom experiments are collected in Table 1, and illustrated in Figures 10 and 11.
The Fly dataset consists of 180 projections with 50 s exposure (detector pixel size 27 mm, demagnified to 15.7 mm 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.4 and relaxing the stationarity assumption by performing PSD estimation and denoising in multiple parts for each sinogram as described in Section 3.3.3. The denoising results of two different sinograms are shown in Figure 12; the corresponding tomograms of the second sinogram of Figure 12 are shown in Figure 13.
We also test the algorithm on a soft tissue sample 00076 displaying severe ring artifacts freely available in TomoBank (De Carlo et al., 2018). The data contain 2000 projections with 2.2 mm 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 mm LuAG scintillator, and sample-to-detector distance as 90 mm. The The complete denoising process, requiring as inputs the noisy projections P raw and the bright-and dark-fields I B , I D , and producing as the output an estimate of the streak-free stack of projections composed of sinogram estimatesŶ Y. 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), with different combinations of stdð P Þ and peak values of A, with peak ¼ 1 being the limiting case for which ¼ 0.
As all of the algorithms aim to remove streak noise only, the SNR values are calculated with Y ¼ ln½A þ =ð1þ P Þ as SNRðŶ YÞ = 10 log 10 ðsvar X fY 2 g= smean X ððŶ 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. 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, and the corresponding tomograms in Figure 15.

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.2.
For the multiscale denoising procedure, we performed vertical binning with m v = dm=dm=64ee ' 64 pixels; this value, being slightly larger than the height of the BM3D search neighborhood of Section 2.2.1 (39 Â 39 pixels), allows our method to attenuate streaks which change slowly across the angle -a larger value of m v might be used to deal with streaks featuring faster variation. The number of horizontal scales K for each sinogram was set as K ¼ blog 2 ðn=40Þc, 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 Z K 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.3), we estimate & k , K ! k ! 0 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 singlethreaded 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   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/2 k th of the cost of k ¼ 0, the total for all k ¼ 0; . . . ; K being at most twice that of single scale denoising of Z 0 . The correlation kernels g c and g B do not depend on the input or scale, and can thus be pre-computed. To compute g B , we use directly the definition (23) 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 m v times in the vertical dimension. For

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) 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    proposed procedure based on BM3D denoising. Although Vo et al. (2018) 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) and Vo et al. (2018) 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. removal. Compared with the two popular streak removal algorithms (Mü nch et al., 2009;Vo et al., 2018), 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.    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 sophisti-cated scale decompositions, such as steerable pyramids, contourlets, or dual-tree wavelets (Kovacevic & Chebira, 2008). 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 nonlinearity 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 et al., 2010) should be used instead.