research papers
Ring artifact and Poisson noise attenuation via volumetric 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 email: ymir.makinen@tuni.fi
Xray microtomography systems often suffer from high levels of noise. In particular, severe ring artifacts are common in reconstructed images, caused by defects in the detector, calibration errors, and fluctuations producing streak noise in the raw sinogram data. Furthermore, the projections commonly contain high levels of Poissonian noise arising from the photoncounting detector. This work presents a 3D multiscale framework for streak attenuation through a purposely designed collaborative filtering of correlated noise in volumetric data. A distinct multiscale denoising step for attenuation of the Poissonian noise is further proposed. By utilizing the volumetric structure of the projection data, the proposed fully automatic procedure offers improved feature preservation compared with 2D denoising and avoids artifacts which arise from individual filtering of sinograms.
Keywords: BM4D; noise; image processing; tomography; artifact.
1. Introduction
Computed tomography is commonly affected by streak noise in measured raw sinogram data (Jha et al., 2013; Artul, 2013; Boas & Fleischmann, 2012), which can be caused by miscalibration of detector linear response, beam fluctuations, beam hardening, or dusty or damaged scintillator screens (Haibel, 2008; Vidal et al., 2005; Anas et al., 2010). Streak noise in projections causes ring artifacts in reconstructed volumes, which present as centered circles or halfcircles (Croton et al., 2019). As the sinogram data are obtained through a photoncounting detector, the statistics of the measured raw data can be further modeled through a which may result in high levels of Poissonian noise, commonly attenuated within the reconstruction process through iterative approaches (Mohan et al., 2014; Venkatakrishnan et al., 2013).
Although ring artifacts can be reduced by scanning protocols (Pelt & Parkinson, 2018), highquality scintillator screens and detectors, it is difficult to completely avoid them and therefore achieve highestquality reconstruction solely by experimental measures, requiring algorithmic processing of the acquisitions.
Popular methods to reduce ring artifacts include waveletFFT 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), smoothing filters with segmentation in the tomogram domain (Massimi et al., 2018), ring removal in the tomogram domain upon polar coordinates transformation (Sijbers & Postnov, 2004; Li et al., 2021), and iterative algorithms (Paleo & Mirone, 2015) that combine regularized reconstruction with denoising. Recently, in Mäkinen et al. (2021), we proposed effective ring artifact attenuation through sinogramdomain collaborative filtering, presenting a multiscale architecture with a Blockmatching and 3D filtering (BM3D) image denoiser for correlated noise (Dabov et al., 2008; Mäkinen et al., 2020) at the core of the process. To the best of our knowledge, Mäkinen et al. (2021) offers stateoftheart results in ring attenuation. In particular, it does not cause new artifacts around strong signal features, common to other popular ring removal algorithms. However, being based on a filter for 2D data^{1}, applied to individual sinograms, it may cause discontinuities across the third dimension.
In this work, we address both streak reduction and Poissonian noise removal from volumetric stacks of projections. The contribution of this work is threefold: (1) We propose a multiscale streak denoising framework for denoising of volumetric data. In particular, this framework can be seen as an extension of Mäkinen et al. (2021) for filtering of 3D volumes. (2) After streak noise removal, and before reconstruction, we embed a distinct multiscale denoising step to attenuate the Poissonian noise component of the projections. This allows to apply the reconstruction process using milder regularization and improve the tradeoff between noise reduction and artifact suppression. (3) As a generalpurpose algorithmic contribution, the filter used at the core of the multiscale denoising process is an improved version of the BM4D (Maggioni et al., 2012) volumetric denoising algorithm. The included enhancements, discussed in Appendix A, allow the longrange noise correlation which characterizes the streaks to be dealt with.
The proposed filtering procedure for both streaks and Poissonian noise is fully automatic and includes selfcalibration of the filtering strength. We demonstrate the denoising performance of the proposed approach on real data from the tabletop Prisma XRM microCT at Sigray, and from the synchrotronbased microCT at the Advanced Photon Source (APS) in Argonne, available through Tomobank (De Carlo et al., 2018).
2. Brightfield normalization
The following normalization of the raw projections and the streak model upon a logarithmic transformation follow that of Mäkinen et al. (2021).
The optical attenuation through the sample is determined experimentally via brightfield corrections through two separate acquisitions, the brightfield and the darkfield (Seibert et al., 1998). The brightfield is obtained by the imaging procedure with no sample, and the darkfield is obtained with no beam; both are 2D arrays the size of effective pixels of the detector. The Beer–Lambert law further relates the Xray 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 darkfield and I_{B} is the brightfield^{2}, and then logtransformed as
2.1. Noise model for normalized projections
Apart from possible completely defective detectors^{3} we treat the variation in detector response as normally distributed; as such, the streak noise will follow a normal distribution. Furthermore, we model the streak noise as locally stationary, meaning that the variance is presumed constant within the support of the denoising filter. Note that this does not mean that the noise is i.i.d. or white, as it is instead characterized by very long range correlation presenting as streaks.
As the data are obtained through a photoncounting detector, the statistics of the measured raw data can be further modeled through a
with nonstationary parameters after the brightfielding.Given the premises of normally distributed streak noise and Poissonian noise, noise in projections normalized by (1) can be formally written as
where A are the underlying noisefree projections, is the normally distributed streak noise component, and π is Poissonian noise which we model as white and zeromean; all components of (3) are considered as 3D arrays and multiplications are elementwise.
The natural logarithm of (2) acts as a variancestabilizing transformation (VST) for the multiplicative noise component . Hence, we have
where the approximation comes from .
3. Correlated noise
The denoising is conducted in two steps. First, we aim to estimate the `streakfree' projections
which are corrupted by white Poissonian noise. Then, as a separate denoising step, we consider the attenuation of the remaining noise originating from π.
Throughout this work, we will represent the volume to be filtered according to the correlated noise model presented in the following subsection. This model will assume different meaning at different parts of the algorithm. First, applying locally to the streaks as a type of longrange correlated noise; second, to the noise arising from the Poissonian component π.
3.1. Correlated noise model
We consider the 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 threedimensional volumetric domain X and
ν being zeromean i.i.d. Gaussian noise with unit variance, and denoting 3D 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 3D Fourier transform, and X denoting the cardinality (i.e. number of elements) of X. Equivalently, a kernel g satisfying (7)–(8) can be defined from Ψ as
3.2. Estimation of noise standard deviation
When applying the above model to noisy data, it is essential to have knowledge of either the kernel g or, equivalently, the PSD Ψ, as they fully characterize the noise. Assuming g in (7) is known modulo a scaling factor ς from a known kernel g_{s}, the noise estimation can be simplified to estimating ς. In particular, in order to model the streak and Poissonian noise components arising from the particular composition of noise given in (3), the kernels g_{s} should induce either very long range correlation or near white noise across each dimension d. The estimation procedure can be built as a direct extension of the one adopted by Mäkinen et al. (2021) to 3D. To reduce the signaltonoise ratio (SNR) to acquire a better noise estimate, we convolve z with a 3D anisotropic kernel g_{d} that provides either lowpass or highpass filtering along different dimensions; g_{d} is designed based on the noise statistics so that it preserves the noise component of interest while attenuating signal contrast. Specific instances are given in Section 4.1.1 and Section 5.3. One can then 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 ς can be obtained through
4. Multiscale streak filtering
In the following, we treat the first dimension of the stack of projections as the angular dimension, and the second and third as the horizontal and vertical displacement dimensions.
Because the streaks are inherently lowfrequency with respect to the angle, they are filtered entirely at a coarse angular scale; for this task, we extend the multiscale procedure of Mäkinen et al. (2021). The main changes in the proposed procedure arise from replacing the onedimensional binning operations along the displacement dimension with corresponding 2D binning operators and executed across both displacement dimensions. Furthermore, instead of using a direct 3D extension of the 2D streak PSD, we adjust the streak model to account for possible long correlation also along the displacement dimensions.
In detail, the multiscale streak attenuation procedure proceeds as follows. We begin by an angular binning . The result of the angular binning Z_{0} = is binned K times through to obtain Z_{K} = . The size of each binned volume is a quarter of the input size. Then, we process each scale in a coarsetofine fashion, where progressively for each k = , we obtain an estimate of = . We start by taking as noisy input Z^{ *}_{K} of BM4D the smallest binned volume 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 k = :
(1) Replace the coarserscale components of Z_{k} by those of the estimate :
(2) Denoise Z^{ *}_{k} with BM4D to produce the estimate .
Finally, we replace the coarse angular components of the fullsize stack Z with those from the finest scale estimate ,
4.1. Multiscale noise model
For BM4D denoising, we regard Z^{ *}_{k} of each scale k as z of the model (6), as
where
and = = .
This definition for , , follows from considering the coarserscale estimate as perfectly denoised. Similar to (8), is treated as correlated noise with PSD,
where g_{k}^{ *} is a correlation kernel and X_{k} is the pixel size of Z_{k}. As per (9), the kernel g_{k}^{ *} can be defined as
4.1.1. Adaptive parametric model of
We note that, in addition to the approximately white streak noise, the sinograms may contain streaks with very long range correlation across the displacement dimensions. As this correlation is aligned along the detector axis, it is not clearly observable in individual sinograms, but may create significant noise structure in the full volume. Hence, we approximate the streak noise η_{0} through three angularly constant streak noise components distinct in the displacement,
where η_{0,w} is streak noise white across both displacement dimensions, η_{0,u} is streak noise constant across horizontal displacement, and η_{0,v} is streak noise constant across vertical displacement.
Let us denote by a noise component of η_{0}. For each η_{0,p}, we can define a respective scaled correlation kernel
where = , and = 1. Example realizations, kernels, and PSDs for each of these components as well as η_{k} are shown in Fig. 1 (top).
We note that each = is characterized by a kernel = . In particular, this property arises from the noise structure of the corresponding components: as operates through summation and elimination of adjacent pixels, the operation preserves both noise whiteness and constant noise. The ratio 2^{k} follows from the summation along two dimensions, meaning that the variance of the coarser scale is four times that of the finer scale.
Disregarding the specific support size of their actual finite realizations, we can identify the stationary random fields as
where η_{G,p} is noise characterized by g_{k,p}, and hence = 1. We can then express the residuals of any of the components η_{0,p} as
Then,
where the noise correlation kernel corresponding to a component can be written with = as
where is a 2D kernel across the displacement dimensions characterizing the residual from 2D binning of white noise. Specifically,
where η_{G,2D} is a twodimensional white random field. The field size X_{G} is included only for the normalization of the Fourier transform, and the formula holds for an arbitrary size.
Then, the PSD of , , can be written as
As any is constant in angle, is nonzero only across the DC plane with respect to the angular dimension. Example realizations, kernels, and PSDs for the residual components are shown in Fig. 1 (bottom).
Although (23) allows for modeling of very long range correlation, the streak noise is likely to contain minor correlation along the displacement not accounted for by this model. To adapt to such deviations, we allow the scaling parameters 0 for each noise component to vary with each scale k by estimating them individually at each scale, effectively accounting for mild local correlation in the noise.
Estimation of ς_{k,w}, ς_{k,u}, and ς_{k,v}. Based on (22) and (23), the PSD is completely determined by the values assumed by the three parameters ς_{k,w}, ς_{k,u}, and ς_{k,v} and the known kernels g_{k,p} and . To adaptively obtain the parameters, we begin by obtaining three noise variance estimates , , and . For each estimate, we define a corresponding filtering kernel g_{d}^{(p)} such that estimates the variance of highfrequency streaks, estimates the variance of horizontally lowfrequency streaks, and of vertically lowfrequency streaks. For this purpose, we define ϕ_{d} as a 1D Gaussian function along dimension d, and ψ_{d} as a 1D highpass kernel with Daubechies wavelet `db3' of length 6 along d. Hence, convolution with ϕ_{d} realizes lowpass filtering, and ψ_{d} realizes a highpass filter. Then, g_{d}^{(p)} is realized as a tensor product of three onedimensional kernels across the dimensions d chosen based on the noise statistics through that dimension,
Specifically, with m_{0}, m_{1}, m_{2} as the pixel sizes of the three dimensions of Z_{k}^{ *}, ϕ_{0} is a 1D Gaussian function along the angular dimension with standard deviation m_{0}/8, and ϕ_{1} and ϕ_{2} are 1D Gaussian functions along the two displacement dimensions with standard deviations of m_{1}/12 and m_{2}/12, respectively. Through these kernels, we obtain estimates of the three coefficients as described in (10) and (11) with g_{s} as either g_{k,p} ( k = K) or ().
We note that these three components do not directly correspond to ς_{k,w}, ς_{k,u}, and ς_{k,v}, as the frequencies of the white streak component partly overlap with those of and , i.e. includes also some lowfrequency streak components. In particular, we have and . To this end, we can formulate a simple nonnegative leastsquares optimization as
Finally, we construct the PSD through (22) and (23) with = .
4.1.2. Nonstationarity of
Variance of the streak noise may differ across the sinogram due to changes in Z^{ *}_{k} assuming an equal for all spatial positions without either oversmoothing or leaving noise artifacts in some areas. To adapt to nonstationarity, we further relax the streak model allowing the PSD to vary within each scale k. In particular, before noise estimation and denoising, we split Z^{ *}_{k} into overlapping, volumetric segments. We apply BM4D separately on each segment of Z^{ *}_{k}, using a PSD scaled by parameters estimated from the same segment, i.e. we consider each segment as a separate noisy volume z with a corresponding Ψ. After denoising, the segment estimates produced by BM4D are recombined with a windowing function to form the full estimate .
or noise in the brightfield. Thus, it may not be possible to denoise4.2. 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 that reasonably produced by the distribution of and therefore require a specific preprocessing. To this end, after the brightfielding and logtransform and before the multiscale denoising procedure, we apply the simple extreme streak attenuation procedure as described by Mäkinen et al. (2021), which applies median filtering on extreme streak values detected through local polynomial fit of angular medians.
5. Poisson denoising
A filter for additive noise is not immediately applicable to the approximately white noise of originating from the Poissonian component π. Firstly, the brightfielding (1) introduces substantial spatial variability in the Poisson model. As a result, for a given optical attenuation, noise in brightfielded projections can be stronger or weaker in different parts of the detector, for example around edges in conebeam acquisition. Secondly, while the logarithm (2) effectively makes the streak noise additive, it also changes the typical affinevariance model of the Poissonian noise to a nonlinear one where the variance is not constant, but asymptotically inversely proportional to the mean. In order to model the noise in through (6), we take care of these two issues as follows.
5.1. Reducing nonstationarity induced by brightfielding
The Poissonian noise component π originates from a counting process which takes place before brightfielding (1), and specifically before the division by , which introduces a spatially variant scaling of the variances. To undo this scaling, we consider
where = . Then, S can be treated as the logscale version of a homogeneous Poissonian process; S is thus subject to signaldependent noise where the variance of the noise can be expressed as a smooth nonnegative function of the underlying signal,
where the same F applies to each pixel. In particular, it can be shown that asymptotically for large
5.2. Stabilization of variance
To turn a model like (25) into (6) we again resort to the use of a VST. As largeflux asymptotics are irrelevant for denoising problems characterized by low signaltonoise ratio, and to pragmatically accommodate for model uncertainties, we model F as a polynomial with arbitrary datadriven coefficients. The method (Foi, 2009; Azzari & Foi, 2014) simultaneously identifies the coefficients for an arbitrary signaldependent noise model where the variance is a positive power of an unknown polynomial, and returns the associated variance stabilizing transformation f as well as the corresponding exactunbiased inverse VST . An example of an estimated standard deviation function and the corresponding VST f are illustrated in Fig. 2, where the effectiveness of the stabilization can be deduced by the estimates of being scattered around 1.
5.3. Multiscale denoising of the stabilized Poisson noise
To avoid introducing structured artifacts that are present in the brightfield and darkfield images, we further consider a smoothed version of the field component,
where g_{I} is a 2D Gaussian kernel, and medfilt denotes a 2D median filter. The median filter is adopted in order to remove extreme outliers (e.g. from broken pixels), and the convolution with the Gaussian ensures a smooth result. Then, can be used for approximate correction for the brightfield induced nonstationarity with .
The stabilized noisy stack can then be written as
where corresponds to the stabilized noise and to the signal upon stabilization.
We consider π white, and assume the streak denoising procedure to remove all streak noise frequencies, including those of π. Hence, we treat as missing the streak frequencies, i.e. with a PSD,
where c is a constant such that = .
For multiscale denoising of the Poisson component, we define threedimensional binning and debinning functions as = and = , and obtain K_{Poi} scales of binned noisy volumes as = , . Then, unlike the progressive denoising of the streaks, we begin by BM4D denoising of of each scale k; at each scale, we model the noise through a PSD of the form (28). This way, we obtain an initial estimate of the corresponding noisefree volume at each scale. Then, starting from , we combine only the denoised volumes of each scale by recursively replacing the lowscale components of , , by those of the lower scale,
where denotes a 3D Gaussian kernel. Although the low frequencies are obviously denoised more effectively in the coarser scale, the higher frequencies of the coarser scale are commonly estimated worse than the respective estimate of the finer scale (Facciolo et al., 2017). As such, realizes a lowpass filtering which selects only low frequencies of the coarser estimate to be used in the full estimate.
To account for possible remaining nonstationarity and slight correlation of the noise, we perform the denoising in segments similar to as described in Section 4.1.2 for streak noise, and estimate a separate scaling parameter in construction of the PSD at each scale. In particular, we estimate as described in Section 3.2 with g_{d} = and g_{s} = defining the unscaled noise correlation kernel, and finally construct the PSD through (28) with c = .
The final estimate of the underlying stack of projections can be obtained by applying to the finest scale estimate and then removing the field ,
As (30) negates the field correction, we note that had we used the nonsmooth field I_{L} in (27) [and respectively in (30)], any noise or spurious structures present in I_{L} could be introduced into , as they might have been denoised by BM4D and hence not preserved in .
Upon variance stabilization, Poissonian data become asymptotically normal (Curtiss, 1943). Due to the additional Gaussianization induced by the binning and by the linear transformations operated by the filter, the assumption of normality in (7) can be adopted for denoising of the Poissonian component in this work even for lowcount data.
6. Experiments
We test our pipeline on synthetic data as well as two real acquisitions displaying ring artifacts and Poisson noise.
As a comparison, we show results for Mäkinen et al. (2021) available on PyPI as bm3dstreakremoval, the proposed algorithm embedding the conventional BM4D denoiser (Maggioni et al., 2012), as well as two leading streakremoval procedures from the tomopy Python 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. These streak denoising algorithms are run with the default parameters provided by the software library. To evaluate the benefit of the proposed Poisson denoising procedure with reconstruction which includes further regularization of the data, we include experments with the iterative Total Variation (TV) reconstruction (Goldstein & Osher, 2009) of Marchesini et al. (2020).
For the synthetic experiments, we replicate the noise generation setup of Mäkinen et al. (2021) on a stack of projections ( 238×181×238 pixels) obtained from a 3D BrainWeb phantom (Cocosco et al., 1997) obtained through a padding and Radon transform upon a sign change and an exponential transformation. Specifically, we regard this stack as the underlying projections A and generate noise according to (3) with g as a constant of size m_{0}×1×1 (equal to g_{0,w} of Fig. 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 [2560, 5120] (higher SNR), [1280, 2560], 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 π = 0 (infinite SNR), thus resulting in a total of 16 combinations of Poisson and streak noise strengths. We do not simulate extreme streaks or the brightfielding ( = 1 and = 0). For the streak removal, we consider as the streakfree yet noisy stack Y.
The results of the phantom experiments^{4} for streak attenuation are collected in Table 1, and, for full denoising, evaluating the reconstructed volumes, in Table 2 using iterative regularized TV reconstruction with optimized regularization parameter strength r. The experiments for both streak and Poisson denoising are illustrated in Figs. 4 and 5. All reconstructions are performed upon a sign change.

The Fly dataset consists of 180 projections with 50 s exposure (detector pixel size 27 µm, demagnified to 15.7 µm by conebeam geometry) collected using a Sigray Prisma Xray microtomography instrument at 34 kV; the detector size is 512×512 pixels. The denoising results for two different sinograms, as well as a corresponding tomogram after streak attenuation, are shown in Fig. 6. A comparison of denoising on a vertical slice of the stack of tomograms is shown in Fig. 7, and a comparison for fully denoised reconstructions is shown in Fig. 8.
We also test the algorithm on a soft tissue sample 00072 displaying severe ring artifacts freely available in TomoBank (De Carlo et al., 2018). The data contain 1500 projections with 1.43 µm pixels, obtained at the Advanced Photon Source, 2BM beamline; other experimental parameters are Xray energy of 20 keV, 10 µm LuAG scintillator, and sampletodetector distance of 15 mm. The detector size is 2160×2560 pixels. Included are ten samples for bright and darkfields, which are averaged to obtain a single brightfield and darkfield. A sinogram and a corresponding tomogram from the denoising results for streak removal are shown in Fig. 9, and slices of the stack of tomograms are compared in Fig. 10. Reconstructions of fully denoised projections are further compared in Fig. 11.
The proposed method achieves superior SNR values in streak removal in all simulated noise experiments. Although the difference to the 2D implementation of Mäkinen et al. (2021) is not immediately visually obvious from individual sinograms or tomograms, the displayed vertical slices of the reconstructed objects show clear improvement in both signal preservation and avoiding discontinuity between different tomograms. Compared with Münch et al. (2009) and Vo et al. (2018), the proposed method avoids creation of shadow artifacts around strong signal features. Furthermore, performing the Poisson denoising through the proposed framework allows application of standard filtered backprojection reconstruction to data originally corrupted by Poisson noise, but can also improve quality of iterative TV reconstruction.
6.1. Parameters
For streak attenuation, we calculate K following the formula of horizontal binning from Mäkinen et al. (2021), using as the base the size of the smallest displacement dimension. As a result, we use K = 5 for 00072, K = 3 for Fly, and K = 2 for the phantom. These values were found to offer a reasonable compromise between denoising wide streaks versus preserving lowfrequency signal components. Other processing parameters are adjusted for the smaller block size and processing neighborhood of BM4D. For angular binning, we use = 32 pixels, where m is the original angular size and m_{α} the output size; the resulting size is half of that used by Mäkinen et al. (2021). For segmentation of the streak denoising, we use a window of size pixels. For the Poisson denoising, we use K = 1 and m×19×19 segments. For variance stabilization, we use the implementation ClipPoisGaus (Azzari & Foi, 2015) of Foi (2009) and Azzari & Foi (2014), and use a quadratic polynomial for the variance model F.
7. Discussion and conclusions
We have presented a framework for threedimensional attenuation of streak noise extending the 2D framework of Mäkinen et al. (2021), as well as a BM4D denoiser utilizing the algorithmic improvements of Mäkinen et al. (2020). Furthermore, we have included a denoising step for Poisson noise in the sinogram domain through BM4D and the adaptive variance stabilization of Foi (2009) and Azzari & Foi (2014).
We test the algorithm on both synthetic and real data, demonstrating superior SNR compared with other popular streak removal algorithms, and showing improvements in streak attenuation over Mäkinen et al. (2020). Furthermore, we compare the results with those which use the conventional BM4D for correlated noise, demonstrating that the included improvements in BM4D for correlated noise are essential for successful streak attenuation. The included Poisson denoising allows for full sinogramdomain denoising within the framework. By operating fully in the 3D stack of projections, the 3D structure of the data can be leveraged for more accurate noise removal. The proposed procedure is fully automatic and does not require extra input parameters.
To compare different methods under their own ideal conditions, we have specifically selected the TV regularization parameter values that provide the best reconstruction quality. However, in realworld applications, these values cannot be identified precisely, and too small or too large parameter values may lead to residual noise or excess smoothing of the reconstructions. Inclusion of the proposed Poisson denoising step allows for weaker regularization, but notably also reduces the effects of relative shifts in the parameter values, meaning that the reconstruction can be safely deployed even when the regularization cannot be precisely tuned.
To consider the computational cost, we note that both denoising steps of Fly ( 181×512×512 pixels) run singlethreaded on an AMD Ryzen 7 1700 processor each take about one hour. The computational cost is mostly due to the BM4D denoising in CPU. Although the adopted implementation is singlethreaded, the algorithm is embarrassingly parallel, and thus a highly parallel GPUbased implementation is expected to reduce the total run time to the scale of seconds (Davy & Ehret, 2020).
The Poissonian noise attenuation can also be performed without the preceding ring reduction step on data which do not display ring artifacts. In such case, should be replaced by a flat PSD, as the Poissonian noise is approximately white prior to streak attenuation, whereas (28) considers the streak noise frequencies removed. Running the full denoising procedure in the absence of either streak or Poisson noise will lead to very small estimates for the corresponding noise components, meaning that no significant denoising will be performed for that noise.
We note that although we have focused on the full denoising of the projections, typical reconstruction pipelines, such as the iterative TV, provide further noise attentuation. For best results in combining the proposed denoising procedure with such pipelines, it may be necessary to adjust the filter strength for the denoising of Poissonian noise, e.g. for reduced attenuation of highfrequency noise, as it is further attenuated within the reconstruction. Likewise, integration of the proposed procedure within an iterative alternating reconstruction is left for future study.
APPENDIX A
Collaborative filtering and the BM4D denoising algorithm
A1. Collaborative filtering
The rationale of transformdomain 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 nonlinear 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 patches extracted from the input. One of the most popular collaborative filters is the BlockMatching and 3D filtering (BM3D) (Dabov et al., 2007) denoising algorithm, which performs denoising on groups of blocks extracted from a 2D image. In the BM4D volumetric denoiser (Maggioni et al., 2012), the patches are 3D volumes extracted from the volumetric data.
All operations of collaborative filters are made with regard to a reference patch moving through the volume. For each position of the reference patch, the following steps are executed:
(1) Collect similar patches into a group through patchmatching.
(2) Obtain a group transform spectrum by collectively transforming the group of patches.
(3) Perform shrinkage.
(4) Transform the shrunk spectra back to patch estimates and aggregate them to the original locations from which they were collected.
For details about the algorithm in arbitrary dimensionality, we refer the reader to Mäkinen et al. (2020). In the following section, we describe special considerations for the implementation of the algorithmic improvements in Mäkinen et al. (2020) for volumetric denoising.
A2. Improved BM4D for correlated noise
Most of the improvements described by Mäkinen et al. (2020) are directly applicable to the 3D denoiser. In this section, we consider extensions which are not immediate from the inclusion of an extra dimension.
A2.1. Shrinkage parameters
For the selection of shrinkage parameters λ and μ^{2}, we embed the parameter selection subroutine of Mäkinen et al. (2020), which is based on the 2D input PSD; for simplicity, we adopt directly also the precomputed features and parameters computed for a set of 2D PSDs. To utilize this system with a 3D PSD, we include a simple procedure which obtains a 2D projection of the 3D PSD by preserving the two largest principal components, aiming to preserve the characterizing features of the PSD shape. This projection is then used to compute features as described for a 2D PSD by Mäkinen et al. (2020) for the estimation of suitable λ and μ^{2}.
A2.2. Fast implementation
We consider the fast implementation as suggested by Mäkinen et al. (2020). In particular, we perform all operations on a downscaled PSD of size N_{f}×N_{f}×N_{f} and compute exactly only the K_{f} first volumes of the 4D spectrum and approximate the rest using the conventional variances. Furthermore, Fourier symmetries and sparsity of the transformed arrays can be exploited to reduce computational cost similar to the 2D case.
A2.3. Refiltering
As noted by Mäkinen et al. (2020), even with exact modeling of the collaborative transformdomain noise spectrum, the accuracy of collaborative filtering is limited by the systemic factors arising from the used transforms, both in size and possible symmetries of the transform spectrum which may limit the modeling of the global PSD. As a result, the denoising may attenuate excess signal, leading to oversmoothing in some frequencies; Mäkinen et al. (2020) proposes the mitigation of these systemic issues through an extra filtering step performed on the denoising residual. The threedimensional spectra are not exempt from these limitations, and as such we adopt the global Fourier thresholding and refiltering procedure through a 3D FFT.
Footnotes
^{1}Despite its name, BM3D is a filter for 2D images, where the third dimension originates from the nonlocal prior leveraged by the algorithm.
^{2}As P_{raw} is a 3D array, the pixels of I_{B} and I_{D} are replicated through the angle dimension for the operations in (1).
^{3}Extreme streak noise arising from defective detectors is addressed separately in Section 4.2.
^{4}For the purpose of visualization and objective SNR comparison, we have corrected the intensity response of every estimate through a cubic polynomial such that it matches the ground truth. This is done particularly to improve the results of those methods such as Münch et al. (2009) whose intensity response significantly deviates from the ground truth.
Acknowledgements
We are very thankful to Dr Sheraz Gul of Sigray Inc. for providing us with the experimental Fly data.
Funding information
This work was supported by the Academy of Finland (project No. 310779) and by the Stanford Synchrotron Radiation Lightsource, SLAC National Accelerator Laboratory, which is supported by the US Department of Energy, Office of Science, under Contract No.DEAC0276SF00515. 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. Web of Science PubMed Google Scholar
Artul, S. (2013). BMJ Case Rep. 2013, bcr2013–201379. Google Scholar
Azzari, L. & Foi, A. (2014). 2014 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP 2014), 4–9 May 2014, Florence, Italy, pp. 5357–5361. IEEE. Google Scholar
Azzari, L. & Foi, A. (2015). ClipPoisGaus: PoissonianGaussian noise estimation and removal for singleimage rawdata, Matlab code, https://webpages.tuni.fi/foi/sensornoise.html#ref_software. Google Scholar
Boas, F. E. & Fleischmann, D. (2012). Imaging Med. 4, 229–240. CrossRef Google Scholar
Cocosco, C. A., Kollokian, V., Kwan, R. K.S., Pike, G. B. & Evans, A. C. (1997). NeuroImage, 5, S425. 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
Curtiss, J. H. (1943). Ann. Math. Stat. 14, 107–122. CrossRef 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
Davy, A. & Ehret, T. (2021). J RealTime Image Proc. 18, 57–74. CrossRef 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
Facciolo, G., Pierazzo, N. & Morel, J.M. (2017). SIAM J. Imaging Sci. 10, 1603–1626. CrossRef Google Scholar
Feldkamp, L. A., Davis, L. C. & Kress, J. W. (1984). J. Opt. Soc. Am. A, 1, 612–619. CrossRef Web of Science Google Scholar
Foi, A. (2009). Signal Process. 89, 2609–2629. CrossRef Google Scholar
Goldstein, T. & Osher, S. (2009). SIAM J. Imaging Sci. 2, 323–343. 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, pp. 141–160. Oxford University Press. 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
Li, Y., Zhao, Y., Ji, D., Lv, W., Xin, X., Zhao, X., Liu, D., Ouyang, Z. & Hu, C. (2021). Phys. Med. Biol. 66, 105011. CrossRef Google Scholar
Maggioni, M., Katkovnik, V., Egiazarian, K. & Foi, A. (2012). IEEE Trans. Image Process. 22, 119–133. CrossRef PubMed Google Scholar
Mäkinen, Y., Azzari, L. & Foi, A. (2020). IEEE Trans. Image Process. 29, 8339–8354. Google Scholar
Mäkinen, Y., Marchesini, S. & Foi, A. (2021). J. Synchrotron Rad. 28, 876–888. CrossRef IUCr Journals Google Scholar
Marchesini, S., Trivedi, A., Enfedaque, P., Perciano, T. & Parkinson, D. (2020). Lecture Notes Comput. Sci. 12137, 248–261. CrossRef Google Scholar
Massimi, L., Brun, F., Fratini, M., Bukreeva, I. & Cedola, A. (2018). Phys. Med. Biol. 63, 045007. Web of Science CrossRef PubMed Google Scholar
Mohan, K. A., Venkatakrishnan, S., Drummy, L. F., Simmons, J., Parkinson, D. Y. & Bouman, C. A. (2014). 2014 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP 2014), 4–9 May 2014, Florence, Italy, pp. 6909–6913. IEEE. Google Scholar
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
Sijbers, J. & Postnov, A. (2004). Phys. Med. Biol. 49, N247–N253. Web of Science CrossRef PubMed Google Scholar
Swinehart, D. F. (1962). J. Chem. Educ. 39, 333. CrossRef Google Scholar
Venkatakrishnan, S. V., Bouman, C. A. & Wohlberg, B. (2013). 2013 IEEE Global Conference on Signal and Information Processing (GlobalSIP 2013), 3–5 December 2013, Austin, Texas, USA, pp. 945–948. IEEE. 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 openaccess article distributed under the terms of the Creative Commons Attribution (CCBY) Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are cited.