Ring artifact and Poisson noise attenuation via volumetric multiscale nonlocal collaborative filtering of spatially correlated noise
X-ray micro-tomography 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 photon-counting detector. This work presents a 3-D 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 2-D denoising and avoids artifacts which arise from individual filtering of sinograms.
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 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). Streak noise in projections causes ring artifacts in reconstructed volumes, which present as centered circles or half-circles (Croton et al., 2019). As the sinogram data are obtained through a photon-counting detector, the statistics of the measured raw data can be further modeled through a Poisson distribution, 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), high-quality scintillator screens and detectors, it is difficult to completely avoid them and therefore achieve highest-quality reconstruction solely by experimental measures, requiring algorithmic processing of the acquisitions.
Popular methods to reduce ring artifacts include 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), 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 sinogram-domain collaborative filtering, presenting a multiscale architecture with a Block-matching and 3-D 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 state-of-the-art 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 2-D data1, 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 3-D 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 general-purpose 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 long-range 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 self-calibration of the filtering strength. We demonstrate the denoising 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).
The optical attenuation through the sample is determined experimentally via bright-field corrections through two separate acquisitions, the bright-field and the dark-field (Seibert et al., 1998). The bright-field is 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. The Beer–Lambert law further relates the X-ray transform through the sample to the optical attenuation by a logarithmic transformation (Swinehart, 1962).
Hence, the raw projections Praw are first normalized as
where ID is the dark-field and IB is the bright-field2, and then log-transformed as
Apart from possible completely defective detectors3 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 photon-counting detector, the statistics of the measured raw data can be further modeled through a Poisson distribution with nonstationary parameters after the bright-fielding.
where A are the underlying noise-free projections, is the normally distributed streak noise component, and π is Poissonian noise which we model as white and zero-mean; all components of (3) are considered as 3-D arrays and multiplications are elementwise.
where the approximation comes from .
The denoising is conducted in two steps. First, we aim to estimate the `streak-free' 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 long-range correlated noise; second, to the noise arising from the Poissonian component π.
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 three-dimensional volumetric domain X and
ν being zero-mean i.i.d. Gaussian noise with unit variance, and denoting 3-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) Ψ,
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 gs, 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 gs 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 3-D. To reduce the signal-to-noise ratio (SNR) to acquire a better noise estimate, we convolve z with a 3-D anisotropic kernel gd that provides either low-pass or high-pass filtering along different dimensions; gd 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
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 low-frequency 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 one-dimensional binning operations along the displacement dimension with corresponding 2-D binning operators and executed across both displacement dimensions. Furthermore, instead of using a direct 3-D extension of the 2-D 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 Z0 = is binned K times through to obtain ZK = . The size of each binned volume is a quarter of the input size. Then, we process each scale in a coarse-to-fine 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 ZK; in this way, we obtain from Z *K = ZK the coarsest estimate , which is taken as initialization for the following recursive steps executed for each scale k = :
(1) Replace the coarser-scale components of Zk by those of the estimate :
(2) Denoise Z *k with BM4D to produce the estimate .
Finally, we replace the coarse angular components of the full-size stack Z with those from the finest scale estimate ,
and = = .
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
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 2k 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 gk,p, and hence = 1. We can then express the residuals of any of the components η0,p as
where the noise correlation kernel corresponding to a component can be written with = as
where is a 2-D kernel across the displacement dimensions characterizing the residual from 2-D binning of white noise. Specifically,
where ηG,2D is a two-dimensional white random field. The field size |XG| 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 non-zero 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 gk,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 gd(p) such that estimates the variance of high-frequency streaks, estimates the variance of horizontally low-frequency streaks, and of vertically low-frequency streaks. For this purpose, we define ϕd as a 1-D Gaussian function along dimension d, and ψd as a 1-D high-pass kernel with Daubechies wavelet `db3' of length 6 along d. Hence, convolution with ϕd realizes low-pass filtering, and ψd realizes a high-pass filter. Then, gd(p) is realized as a tensor product of three one-dimensional kernels across the dimensions d chosen based on the noise statistics through that dimension,
Specifically, with m0, m1, m2 as the pixel sizes of the three dimensions of Zk *, ϕ0 is a 1-D Gaussian function along the angular dimension with standard deviation m0/8, and ϕ1 and ϕ2 are 1-D Gaussian functions along the two displacement dimensions with standard deviations of m1/12 and m2/12, respectively. Through these kernels, we obtain estimates of the three coefficients as described in (10) and (11) with gs as either gk,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 low-frequency streak components. In particular, we have and . To this end, we can formulate a simple non-negative least-squares optimization as
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 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 .
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 pre-processing. To this end, after the bright-fielding and log-transform 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.
A filter for additive noise is not immediately applicable to the approximately white noise of originating from the Poissonian component π. Firstly, the bright-fielding (1) introduces substantial spatial variability in the Poisson model. As a result, for a given optical attenuation, noise in bright-fielded projections can be stronger or weaker in different parts of the detector, for example around edges in cone-beam acquisition. Secondly, while the logarithm (2) effectively makes the streak noise additive, it also changes the typical affine-variance 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.
The Poissonian noise component π originates from a counting process which takes place before bright-fielding (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 log-scale version of a homogeneous Poissonian process; S is thus subject to signal-dependent 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 flux
To turn a model like (25) into (6) we again resort to the use of a VST. As large-flux asymptotics are irrelevant for denoising problems characterized by low signal-to-noise ratio, and to pragmatically accommodate for model uncertainties, we model F as a polynomial with arbitrary data-driven coefficients. The method (Foi, 2009; Azzari & Foi, 2014) simultaneously identifies the coefficients for an arbitrary signal-dependent 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 exact-unbiased 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.
To avoid introducing structured artifacts that are present in the bright-field and dark-field images, we further consider a smoothed version of the field component,
where gI is a 2-D Gaussian kernel, and medfilt denotes a 2-D 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 bright-field 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 three-dimensional binning and debinning functions as = and = , and obtain KPoi 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 noise-free volume at each scale. Then, starting from , we combine only the denoised volumes of each scale by recursively replacing the low-scale components of , , by those of the lower scale,
where denotes a 3-D 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 low-pass 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 gd = and gs = 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 non-smooth field IL in (27) [and respectively in (30)], any noise or spurious structures present in IL 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 low-count data.
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 bm3d-streak-removal, the proposed algorithm embedding the conventional BM4D denoiser (Maggioni et al., 2012), as well as two leading streak-removal 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 3-D 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 m0×1×1 (equal to g0,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 bright-fielding ( = 1 and = 0). For the streak removal, we consider as the streak-free yet noisy stack Y.
The results of the phantom experiments4 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 cone-beam geometry) collected using a Sigray Prisma X-ray micro-tomography 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, 2-BM beamline; other experimental parameters are X-ray energy of 20 keV, 10 µm LuAG scintillator, and sample-to-detector distance of 15 mm. The detector size is 2160×2560 pixels. Included are ten samples for bright- and dark-fields, which are averaged to obtain a single bright-field and dark-field. 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 2-D 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 back-projection reconstruction to data originally corrupted by Poisson noise, but can also improve quality of iterative TV reconstruction.
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 low-frequency 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.
We have presented a framework for three-dimensional attenuation of streak noise extending the 2-D 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 sinogram-domain denoising within the framework. By operating fully in the 3-D stack of projections, the 3-D 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 real-world 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 single-threaded 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 single-threaded, the algorithm is embarrassingly parallel, and thus a highly parallel GPU-based 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 high-frequency 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.
Collaborative filtering and the BM4D denoising algorithm
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 patches extracted from the input. One of the most popular collaborative filters is the Block-Matching and 3-D filtering (BM3D) (Dabov et al., 2007) denoising algorithm, which performs denoising on groups of blocks extracted from a 2-D image. In the BM4D volumetric denoiser (Maggioni et al., 2012), the patches are 3-D 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 patch-matching.
(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.
Most of the improvements described by Mäkinen et al. (2020) are directly applicable to the 3-D denoiser. In this section, we consider extensions which are not immediate from the inclusion of an extra dimension.
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 2-D input PSD; for simplicity, we adopt directly also the pre-computed features and parameters computed for a set of 2-D PSDs. To utilize this system with a 3-D PSD, we include a simple procedure which obtains a 2-D projection of the 3-D 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 2-D PSD by Mäkinen et al. (2020) for the estimation of suitable λ and μ2.
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 Nf×Nf×Nf and compute exactly only the Kf first volumes of the 4-D 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 2-D case.
As noted by Mäkinen et al. (2020), even with exact modeling of the collaborative transform-domain 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 three-dimensional spectra are not exempt from these limitations, and as such we adopt the global Fourier thresholding and refiltering procedure through a 3-D FFT.
1Despite its name, BM3D is a filter for 2-D images, where the third dimension originates from the nonlocal prior leveraged by the algorithm.
2As Praw is a 3-D array, the pixels of IB and ID are replicated through the angle dimension for the operations in (1).
3Extreme streak noise arising from defective detectors is addressed separately in Section 4.2.
4For 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.
We are very thankful to Dr Sheraz Gul of Sigray Inc. for providing us with the experimental Fly data.
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.DE-AC02-76SF00515. Partial support by Hong Kong RGC 14302920.
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, bcr-2013–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: Poissonian-Gaussian noise estimation and removal for single-image raw-data, 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 Real-Time 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 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.