## research papers

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

^{a}Tampere University, Finland, and ^{b}SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025, USA^{*}Correspondence e-mail: ymir.makinen@tuni.fi

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.

Keywords: BM3D; noise; image processing; tomography; artifact.

### 1. 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, 2008), leveraging the recent inclusion of exact transform-domain noise variances (Mäkinen *et al.*, 2020), 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).

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

#### 2.1. Correlated noise model

We consider a noisy input to be a combination of underlying data *y* and additive stationary spatially correlated noise η to be filtered,

where 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 . An equivalent way of representing correlated noise is by its power spectral density (PSD) Ψ,

with 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.

#### 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) 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 locally to each block, then a 1-D transform through the `stack' of grouped blocks. Denoting by a group of *M* blocks of *N* pixels extracted from *z* at coordinates *x*_{1},…, *x*_{M}, we obtain the spectrum coefficients as *s*_{i}^{xt} = , for *i* = 1,…, *N, t* = 1,…, *M*, where is the *i*th basis function of . The 3-D spectrum coefficients are calculated through the direct tensor product of the and transforms, as

where is the *t*th element of the *j*th basis function of , ⊗ 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 play a key role; we denote them by . 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 is the reference block, *z*_{xj} is a potential match, *v*_{i,2}^{ xR, xj} is the *i*th coefficient of the block-pair transform-domain variance corresponding to block difference, and . 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 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 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 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 , 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 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 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,

where *Q*^{2D} is the inverse transform of , and is the *j*th transform basis function of the inverse of .

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 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 and they can all be distinct.

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

where is a block-specific weight and *W*_{xj} is a windowing function over blocks at position *x*_{j}. The weights , 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.

### 3. 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.

#### 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). 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 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, 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 . Hence, we have

where the approximation comes from . The additive noise component in (15) corresponds to the streak noise to be denoised. As here we only aim to attenuate the streak noise, through denoising we estimate ; the embedded noise term although not i.i.d. is nevertheless (approximately) white and does not present streaks.

Individual sinograms, each of which is defined as a *P*_{log}, are denoted as

where *Y* denotes the underlying streak-free sinogram, and is the corresponding of . The sinograms *Z* are used as the input for the processes in the following Section 3.2.

#### 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 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), 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 through a binning operator . 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 : *Z*_{k} = = = , . Denoting by *n* the width of *Z* and *Z*_{0}, *Z*_{k} has width ⌈2^{−k}*n*⌉; with every binning, the streak width also gets halved. The multiscale denoising is operated in a coarse-to-fine fashion, where progressively for each we obtain an estimate of . 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 , which is taken as initialization for the following recursive steps executed for each scale :

(1) Replace the horizontal coarser-scale components of *Z*_{k} by those of the estimate ,

(2) Denoise *Z*^{ *}_{k} with BM3D to produce the estimate .

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

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

#### 3.3. Multiscale noise model

For BM3D denoising, we regard *Z*^{ *}_{k} of each scale *k* as *z* of the model (1), as

where

and .

This definition for , , follows from considering the coarser-scale estimate as perfectly denoised. Similar to (3), is treated as correlated noise with PSD

where *g*_{k}^{ *} is a correlation kernel and . As per (4), the kernel *g*_{k}^{ *} can be defined as

##### 3.3.1. Multiscale PSD of white streak noise

Let be horizontally white and vertically constant streak noise like in Figure 1. Under this assumption for η_{0}, we have that also all for are horizontally white and vertically constant, with variance . 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 . We can hence rewrite (18) as

This together with (20) means that we can characterize through a kernel *g*_{k}^{ *} obtained by scaling either of two basic two kernels *g*_{c} and ,

by a factor

as

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, with . Example noise , , and the corresponding kernel and PSD (19) are shown in Figure 7.

*Estimation of ς _{0}.* To estimate , we first convolve

*Z*

_{0}with a 2-D kernel 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}, offers a lower signal-to-noise ratio (SNR), which facilitates the estimation of noise statistics; an example of

*Z*

_{0}and the corresponding are shown in Figure 8 (top). One can compute an estimate of the standard deviation of via its median absolute deviation (Hampel, 1974),

where smed denotes the sample median and the factor 1.4826 calibrates the estimate with respect to a normal distribution of the noise. As = , an estimate of ς_{0} can be obtained through

##### 3.3.2. Adapting the model to non-white streak noise

The above model (21)–(25) assumes that the streak noise 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 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}.

*Estimation of ς _{k}, .* For

*k*=

*K*, ς

_{K}can be estimated from

*Z*

_{K}=

*Z*

^{ *}

_{K}by trivial substitutions of 0 with

*K*in (26) and (27). Although also for 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 offers an even lower SNR than . An example of

*Z*

^{ *}

_{k}and the corresponding are shown in Figure 8 (bottom). Similar to (26) and for any , the standard deviation of the noise in can be estimated as

Noting that , we then estimate ς_{k} as

##### 3.3.3. Horizontal nonstationarity of η_{Z}

Variance of the streak noise may differ across the sinogram due to changes in *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 into overlapping, full-height segments. We apply BM3D separately on each segment of *Z*^{ *}_{k}, using a PSD scaled by estimated on the corresponding segment of , *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 .

#### 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 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 of , we fit a bivariate cubic polynomial ℘_{x} to a window of centered at *x*. Then, consistent with Gaussian modeling of , we mark the center pixel defective if , where sstd denotes the sample standard deviation; each marked pixel in 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.

### 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): Münch *et al.* (2009) and 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 is generated with = 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 , thus defining the Poissonian noise π as the difference between this Poisson variate and . Furthermore, we include experiments with (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 and we consider 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 µ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 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 µ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, and the corresponding tomograms in Figure 15.

#### 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.2.

For the multiscale denoising procedure, we performed vertical binning with = ≃ 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 , 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}, 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). 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/2^{k}th of the cost of *k* = 0, the total for all being at most twice that of single scale denoising of *Z*_{0}.

The correlation kernels *g*_{c} and do not depend on the input or scale, and can thus be pre-computed. To compute , 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 Münch *et al.* (2009) and Vo *et al.* (2018), we use implementations *remove_stripe_fw* and *remove_all_stripe* provided by the *tomopy* Python library of Gürsoy *et al.* (2014). The tomograms of each experiment are reconstructed using the *xpack* library (Marchesini *et al.*, 2020).

### 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) 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; 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.* (2009) and Vo *et al.* (2018) 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). 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.

### Footnotes

^{1}As *P*_{raw} is a 3-D array, the pixels of *I*_{B} and *I*_{D} are replicated through the angle dimension for the operations in (12).

^{2}Extreme streak noise arising from defective detectors is addressed separately in Section 3.4.

^{3}Adopting the kernels (25) with arbitrary values of ς_{k} corresponds to assuming η_{k} approximately white with variance , which may differ from 2^{k}std(η_{0}) against (21). This assumption becomes increasingly appropriate as *k* grows for non-white η_{0} featuring mild local horizontal correlations, as the binning in 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

Abu Anas, E. M., Lee, S. Y. & Hasan, M. K. (2010). *Phys. Med. Biol.* **55**, 6911–6930. PubMed Google Scholar

Artul, S. (2013). *BMJ Case Rep.* **2013**, bcr-2013–201379. Google Scholar

Azzari, L. & Foi, A. (2016). *IEEE Signal Process. Lett.* **23**, 1086–1090. CrossRef Google Scholar

Boas, F. E. & Fleischmann, D. (2012). *Imaging Med.* **4**, 229–240. CrossRef Google Scholar

Croton, L. C., Ruben, G., Morgan, K. S., Paganin, D. M. & Kitchen, M. J. (2019). *Opt. Express*, **27**, 14231–14245. Web of Science CrossRef CAS PubMed Google Scholar

Dabov, K., Foi, A., Katkovnik, V. & Egiazarian, K. (2007). *IEEE Trans. Image Process.* **16**, 2080–2095. Web of Science CrossRef PubMed Google Scholar

Dabov, K., Foi, A., Katkovnik, V. & Egiazarian, K. O. (2008). *Proc. SPIE*, **6812**, 681207. CrossRef Google Scholar

Davidson, 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

Davy, A. & Ehret, T. (2020). *J. Real-Time Image Process.* https://doi.org/10.1007/s11554-020-00945-4. Google Scholar

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. Web of Science CrossRef Google Scholar

Gürsoy, D., De Carlo, F., Xiao, X. & Jacobsen, C. (2014). *J. Synchrotron Rad.* **21**, 1188–1193. Web of Science CrossRef IUCr Journals Google Scholar

Haibel, A. (2008). *Advanced Tomographic Methods in Materials Research and Engineering, Monographs on the Physics and Chemistry of Materials*, pp. 141–160. Google Scholar

Hampel, F. R. (1974). *J. Am. Stat. Assoc.* **69**, 383–393. CrossRef Google Scholar

Jha, 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

Kovacevic, J. & Chebira, A. (2008). *An Introduction to Frames.* Now Publishers Inc. Google Scholar

Mäkinen, Y., Azzari, L. & Foi, A. (2020). *IEEE Trans. Image Process.* **29**, 8339–8354. Google Scholar

Mäkitalo, M. & Foi, A. (2011). *IEEE Trans. Image Process.* **20**, 99–109. PubMed Google Scholar

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. Google Scholar

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 (https://doi.org/10.1007/978-3-030-50371-0_18).

Münch, B., Trtik, P., Marone, F. & Stampanoni, M. (2009). *Opt. Express*, **17**, 8567–8591. Web of Science PubMed Google Scholar

Paleo, P. & Mirone, A. (2015). *J. Synchrotron Rad.* **22**, 1268–1278. Web of Science CrossRef IUCr Journals Google Scholar

Pelt, D. M. & Parkinson, D. Y. (2018). *Meas. Sci. Technol.* **29**, 034002. Web of Science CrossRef Google Scholar

Seibert, J. A., Boone, J. M. & Lindfors, K. K. (1998). *Proc. SPIE*, **3336**, 348–354. CrossRef Google Scholar

Swinehart, D. F. (1962). *J. Chem. Educ.* **39**, 333. CrossRef Google Scholar

Vidal, 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

Vo, N. T., Atwood, R. C. & Drakopoulos, M. (2018). *Opt. Express*, **26**, 28396–28412. Web of Science 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.