 TECHNICAL ADVANCE
 Open Access
 Open Peer Review
 Published:
A sinogram denoising algorithm for lowdose computed tomography
BMC Medical Imaging volume 16, Article number: 11 (2016)
Abstract
Background
From the viewpoint of the patients’ health, reducing the radiation dose in computed tomography (CT) is highly desirable. However, projection measurements acquired under lowdose conditions will contain much noise. Therefore, reconstruction of highquality images from lowdose scans requires effective denoising of the projection measurements.
Methods
We propose a denoising algorithm that is based on maximizing the data likelihood and sparsity in the gradient domain. For Poisson noise, this formulation automatically leads to a locally adaptive denoising scheme. Because the resulting optimization problem is hard to solve and may also lead to artifacts, we suggest an explicitly local denoising method by adapting an existing algorithm for normallydistributed noise. We apply the proposed method on sets of simulated and real conebeam projections and compare its performance with two other algorithms.
Results
The proposed algorithm effectively suppresses the noise in simulated and real CT projections. Denoising of the projections with the proposed algorithm leads to a substantial improvement of the reconstructed image in terms of noise level, spatial resolution, and visual quality.
Conclusion
The proposed algorithm can suppress very strong quantum noise in CT projections. Therefore, it can be used as an effective tool in lowdose CT.
Background
Computed tomography (CT) is one of the most widely used imaging modalities in medicine and its usage has been consistently growing over the past decades [1]. Although CT is an indispensable tool in modern medicine, one of its drawbacks is that the radiation used in this type of imaging can be harmful to the patient health. Reducing the radiation dose requires reducing the number of projection measurements or reducing the radiation intensity, which in turn results in noisier measurements. Reconstructing a highquality CT image from such measurement is a great challenge. The traditional image reconstruction methods in CT are based on filtered backprojection (FBP). Although FBPbased methods are fast and relatively easy to implement, they perform very poorly when the number of projections is reduced or when the measurements are very noisy. With increased awareness of the harmful effects of excessive radiation, recent years have witnessed a new interest in statistical reconstruction methods in CT. These methods promise reconstructing a highquality image from undersampled or noisy projections. Statistical image reconstruction methods in CT are based, as their name implies, on including the statistics of the detected photon counts in the reconstruction process. They can be classified into three categories [2, 3]:

(A)
Methods that work on the raw data (sinogram) these are mostly smoothing techniques that aim at reducing the Poisson noise in the projection data. The processed sinogram data are then used to reconstruct the image using an iterative or filtered backprojection (FBP) method.

(B)
Methods that work on the reconstructed image these are also denoising techniques, but operate in the image domain.

(C)
Methods that use statistical modeling of the imaging process to devise iterative algorithms for image reconstruction.
Whereas the methods in category (C) usually lead to better image quality, the first two types of methods are usually easier to implement, faster, and independent of the image reconstruction method. Among the first two categories of methods, category (A) is usually much more effective. This is because while we have a good understanding of the properties of the noise in the raw (i.e. projection) data, the statistics of the noise in the reconstructed image are not as well understood. Hence, methods that filter the reconstructed image are postprocessing algorithms that cannot use the photon statistics and their performance is in general inferior to filtering in the raw data domain.
Previous studies in sinogram denoising are numerous and they present a diverse set of possible approaches to the problem. Because the variance of the Poisson noise changes with the signal amplitude, adaptive filters are a natural choice and they have been researched extensively [4, 5]. Bayesian approaches, usually maximizing a weighted sum of the data likelihood and a regularity term, have also been proposed in several studies [6, 7]. These algorithms can be very elaborate because the resulting optimization problem is usually hard to solve. At the other end of the spectrum of algorithmic complexity, there exist very simple methods involving shiftinvariant lowpass filters that reduce the highfrequency variations in the sinogram [8, 9]. As one might expect, these methods are much less effective because they do not consider the signaldependent nature of the noise. Another class of algorithms that have been used for sinogram denoising include multiresolution methods involving shorttime Fourier and wavelet transforms [10, 11]. These methods are based on a standard signal processing approach of separating the noise from the signal in the transform domain. A number of studies have focused on denoising the projection data in the Radon space (i.e. after the logarithm transformation). Even though the photon counts follow a Poisson distribution, the noise in the Radon space follows an approximately Gaussian distribution [12, 13]. Therefore, some studies have proposed algorithms for noise suppression in the Radon space based on the assumption that the noise is normally distributed [14, 15]. The Gaussianity assumption significantly simplifies the problem and widens the range of tools that can be employed in algorithm development. However, as the photon counts decrease, the distribution of the noise in the Radon space can significantly depart from a Gaussian distribution [13]. Therefore, these methods become less efficient in lowdose setting, which is exactly where noise suppression is mostly needed.
In the past decade, image denoising has witnessed a growing interest in patchbased methods. The two main classes of patchbased denoising methods include dictionarybased denoising methods and nonlocalmeans methods. Dictionarybased denoising methods are based on the assumption that small patches of natural images have a sparse representation in a (usually overcomplete) dictionary that can be learned from training data [16, 17]. As opposed to more traditional denoising methods that use an offtheshelf basis (e.g., a wavelet basis), these dictionarybased methods adapt the dictionary to the specific class of images at hand. If the dictionary is welltrained, genuine image features will have a sparse representation in the dictionary, whereas noise will not have such a sparse representation. Therefore, sparse representation of image patches will lead to denoising. Dictionarybased denoising methods have shown to be highly successful in denoising of general natural images [16, 17] as well as medical images [18, 19]. The nonlocalmeans methods, on the other hand, exploit the selfsimilarities in natural images. To estimate the denoised value of each pixel, they consider a small patch centered on that pixel and search the image for patches that are similar to this patch. The true (i.e., denoised) value of the pixel is estimated by some collaborative filtering of the found similar patches. This collaborative filtering was in the form of a weighted averaging in the original nonlocalmeans algorithm [20, 21], but has much more elaborate forms in more recent algorithms [22]. Although the origins of these denoising methods go back approximately 10 years, they have been extended and improved in many ways and they present some of the best available image denoising methods. One of the limitations of patchbased methods, especially for large 3D images that are common in medical applications, is their high computational demands.
In this study, we suggest smoothing the noisy sinograms by minimizing a cost function that consists of a data likelihood term and a regularization term that promotes gradient sparsity. The resulting optimality condition suggests an adaptive filter that carries out a stronger denoising where the signal intensity is higher, which is the expected outcome under Poisson distribution. Instead of solving the resulting optimization problem directly, we suggest modifying an existing algorithm such that it approximates the exact solution locally. Therefore, our approach computes the denoised value of each sinogram pixel by solving a local optimization problem in a small neighborhood around that pixel. We will evaluate the performance of the suggested approach by applying it to noisy simulated projections and lowdose projections acquired from a microCT scanner.
Methods
Given a noisy image, v, the goal is to find an image u that ideally preserves all important features in the image while reducing the unwanted noise. Since this is an illposed inverse problem, additional information needs to be introduced to regularize the problem. The denoised image can then be found by minimizing a cost function of the form:
The first term in E _{ λ }(u) reflects the requirement that the denoised image should be close to the noisy image and the second term is the regularization. In the context of denoising problem formulated as (1), a proper regularization is necessary because without a regularizing term the noisy image itself is the only solution of the problem; i.e., if λ=0 then u=v is the only minimizer of E _{ λ }(u). A famous example of variational denoising methods is the RudinOsherFatemi (ROF) model that has the following form [23]:
where Ω denotes the image domain and λ is the regularization parameter. The second term, where ∇u is the gradient of u, is the regularization term. The choice of ℓ _{2}norm in the first term stems from the assumption of additive Gaussian noise. For the case of Poisson distribution considered in this study, this term has to be modified.
Let u(i,j) and v(i,j) denote pixels of discretized versions of u and v. For an arbitrary pixel location (from the probability mass function of a variable with Poisson distribution):
assuming the pixel values are independent, for the entire image we will have:
We ignore the denominator, which is independent of u. Since we want to find a functional to minimize, we consider the negative logarithm of the numerator:
With this modified fidelity term, the new energy functional to be minimized will have the following form [24]:
Now, compare the optimality conditions for the two models (obtained from the EulerLagrange equation):
where p is a subgradient of \(\int _{\Omega } \nabla u\). The only difference between the two equations is the dependence of the regularization parameter on u in the new model. The new model suggests that a stronger smoothing must be applied where signal has larger values. This outcome agrees with what we expect since under Poisson distribution noise level is proportional to signal intensity.
Minimization of the new functional in (6) is not a straightforward problem. One approach is to first replace ∇u with \(\sqrt {\nabla u^{2}+\epsilon }\) for a small ε>0 and then to apply a gradient descent iteration [24]. Another approach, suggested in [25], is to use a Taylor’s expansion of the data fidelity term and to minimize this approximate model. In this paper, we follow an algorithm that was developed for the original ROF model [26]. However, we denoise each sinogram pixel separately by minimizing E _{ λ }(u) in a small neighborhood around that pixel, and with a regularization parameter inspired by the new model.
Total variation denoising has proved to be very effective when applied to piecewiseconstant images. On more complex images, however, it can lead to undesired effects. Most importantly, on images with piecewiselinear features, i.e. regions with smooth change in intensity, ROF’s original model leads to staircase artifacts. This is because, by design, ROF’s model favors piecewiseconstant solutions. We believe that this can be a major problem when applying TV denoising to sinogram images because even if the imaged object is piecewiseconstant, its projections can be very far from piecewiseconstant. This is easy to visualize because a feature with uniform intensity in the imaged object will have a piecewiseconstant projection in the sinogram only if different rays from the xray source to the detector plane cut the same length through that feature, which is a very unlikely scenario. Hence, the projections are very likely to contain features of higher orders of regularity (i.e., piecewiselinear, piecewisequadratic, etc.) that would suffer from a staircase effect when treated with the ROF model. A heavily researched approach to reducing the staircase effect is to replace the ℓ _{1} norm of the gradient with the ℓ _{1} norm of higherorder differential operators [27, 28]. A less sophisticated approach, but one that has a trivial implementation, is to perform the total variation minimization locally. This approach has also been shown to alleviate the staircase effect [29]. Moreover, with a local minimization strategy, if the size of the neighborhood considered in minimization is small enough, one can safely assume that the signal intensity and noise level are approximately constant. Therefore, a solution based on the ROF’s original model will be a good approximation to the solution of the new model based on Poisson noise. This way, we can utilize efficient existing algorithms for the ROF model while avoiding the staircase artifacts.
Since our approach is based on Chambolle’s famous algorithm [26], we briefly describe this algorithm here. If we denote by X and Y the space of the image u and its gradient, ∇u, respectively, then an alternative definition of total variation of u is:
Chambolle introduced the discrete divergence operator as the dual of the gradient operator, i.e. 〈p,∇u〉_{ Y }=〈−div p,u〉_{ X }. In the discrete image domain:
Because of the duality of the gradient and divergence operators, total variation can also be written as:
The minimizer of the energy functional in (2) is then obtained by projecting v onto the set λ K:
which is equivalent to minimizing the Euclidian distance between v and λ div p, and this can be achieved via the following iteration for computing p:
where τ>0 is the step size. For a small enough step size, τ≤1/8, the algorithm is guaranteed to converge [26].
Instead of a global solution, we minimize the energy functional (6) in a small neighborhood of each pixel. To this end, let us denote by ω the set of indices that define the desired neighborhood around the current pixel. For example, for a square neighborhood of size (2m+1)×(2m+1) pixels that we used in this study, \(\omega =\left \lbrace (i,j): \, i,j=m:m \right \rbrace \). We also consider a normalized Gaussian weighting function on this neighborhood:
The local problem will then become that of minimizing the following functional:
where \(._{W}^{2}\) denotes the weighted norm with weights W and v _{ ω } denotes the noisy image restricted to the window ω around the current pixel. It must be clear that here u ^{′} is a function defined only on the neighborhood ω. The solution of this local optimization problem will be similar to Chambolle’s algorithm described above [29]. The only difference is in the update formula for p:
where D is a diagonal matrix whose diagonal elements are the values of W.
The regularization parameter, λ ^{′}, must be chosen according to (7). The simplest approach is to set λ ^{′}= λ v(i,j), where λ is a global regularization parameter and v(i,j) is the value of the current pixel in the noisy image. Since v(i,j) is noisy, a better choice is to use a weighted local average as the estimate of the intensity of the true image at the current pixel (note that the maximumlikelihood estimate of the mean of a Poisson process from a set of observations is the arithmetic mean of the observations). Therefore, we suggest the following choice for the local regularization parameter.
There are several parameters in the proposed algorithm. The global regularization parameter λ controls the strength of the denoising. It should be set based on the desired level of smoothing. Parameter m sets the size of the neighborhood considered around each pixel, which in this study was chosen to be a square window of size (2m+1)×(2m+1). Numerical experiments in [29] have shown that in total variation denoising, the influence map of a pixel is usually limited to a radius of approximately 10 pixels for typical values of the regularization parameter. Therefore, a good value for m would be around 10, which is the value we used for all experiments reported in this paper. The width of the Gaussian weighting function W is adjusted through h. We used h=2m which we found empirically to work well. Similarly, a and h ^{′} in (16) determine the size of the window and the weights used for determining the local regularization parameter. These have to be set based on the noise level in the image; larger values should be chosen when noise is stronger. In this study, we used a=4 and h ^{′}=2a.
A simple implementation of the proposed algorithm can be computationally intensive because it will involve solving a minimization problem, though very small, for every individual pixel in the sinogram. This will be a major drawback because a big advantage of sinogram denoising methods, compared to iterative image reconstruction methods, is the shorter computational time. To reduce the computational time, after minimizing the local cost function (14) around the current pixel, we will replace the value of all pixels in the window of size (2a+1)×(2a+1) around the current pixel, instead of just the center pixel, and then shift the window by (2a+1). Our extensive numerical experiments with simulated and real projections showed that with this approach the results will be almost identical to the case where only one pixel is denoised at a time. This is the approach that we followed in all experiments reported in this paper.
Results and discussion
We evaluated the proposed denoising algorithm on simulated projections and two sets of real lowdose projections of a physical phantom and a rat obtained using a microCT scanner. We compared the performance of our proposed algorithm with two other methods:

1.
A bilateral filtering algorithm proposed in [14]. This bilateral filtering algorithm works by minimizing a cost function of the following form:
$$ E(u)= \sum\limits_{k} \sum\limits_{k' \in \Omega_{k}} P_{1}(k,k') P_{2}(k,k') $$where k denotes the index of a pixel, Ω _{ k } is a neighborhood around this pixel, and P _{1} and P _{2} are two cost functions in terms of the spatial distance and difference in pixel values, respectively, both having Gaussian forms:
$$ P_{1}(k,k')= \exp \left( \frac{\left(kk'\right)^{2}}{2 d^{2}} \right) $$$$ P_{2}(k,k')= 1 \exp \left( \frac{\left(u_{k}u_{k'}\right)^{2}}{2 \sigma^{2}} \right) $$The algorithm parameters include d and σ that control the range of weighting for P _{1} and P _{2}. The authors in [14] suggest a fixed filter length of w=5 and fixed d=w/6=5/6 (that they found empirically to be a good choice) and try a range of values for σ between 0.7 and 2.8. In this study, we applied the bilateral filtering for several values of σ in this range and apply the proposed algorithm for several values of the regularization parameter, λ.

2.
A nonlocal principal component analysis (NLPCA) algorithm proposed in [30]. In this method, patches of the image are first clustered using the KMeans algorithm. For all patches in a cluster a Poisson PCA (also known as exponential PCA) is performed to denoise them. The PCA problem is solved via alternating Newton’s iterations. The denoised patches are returned to their original locations and averaged (to account for the patch overlaps) in order to form the denoised image. Patchbased methods are computationally very intensive. Therefore, with this algorithm we used parameter settings that resulted in a reasonable computational time.
Simulated data
We simulated 360 noisy conebeam projections, from 0° to 359° with 1° increments from a 3D SheppLogan phantom according to the following model [12, 31]:
where \({N_{d}^{i}}\) and \({N_{0}^{i}}\) denote, respectively, the number of detected and incident photons for the ray that extends from the source to the detector i and \(\int _{i} \mu ds\) is the line integral of the attenuation coefficient along that ray. We assumed \({N_{0}^{i}}\) to be constant for all i, i.e., no bowtie filtration, and \({N_{d}^{i}}\) to be a Poissondistributed random variable with the expected value given by (17). We used two values of \({N_{0}^{i}}= 500\) and 2000 to simulate two sets of projections, which we will call highnoise and lownoise, respectively. The phantom size was 512×512×512 voxels and the projections were each 700×700 pixels in size.
Figure 1 shows onedimensional profiles of the noisy and denoised projections. The plots in this figure show that the proposed TVbased denoising significantly removes the noise and seems to be superior to bilateral filtering and NLPCA. For quantitative comparison, because we have access to the true noisefree projections, we use the following criteria:

(a)
Root Mean Square of the Error (RMSE), where error is defined as the difference between the denoised and the true (i.e., noisefree) projections.

(b)
Mutual information (MI), which treats the projections as stochastic processes [32, 33]:
$$ \text{MI}(u^{*},\hat{u})= \sum\limits_{i=1}^{h} \sum\limits_{j=1}^{h} p(u^{*}_{i}, \hat{u}_{j}) \text{log} \left(\frac{p(u^{*}_{i}, \hat{u}_{j})}{p(u^{*}_{i}) p(\hat{u}_{j})} \right) $$Here, u ^{∗} and \(\hat {u}\) represent the true and denoised projections, respectively. We used histograms of u ^{∗} and \(\hat {u}\) for estimating \(p(u^{*}_{i})\) and \(p(\hat {u}_{j})\) and their joint histogram for estimating \(p(u^{*}_{i},\hat {u}_{j})\), and h is the number of bins in the histograms. We normalized the computed \(\text {MI}(u^{*},\hat {u})\) by dividing it by MI(u ^{∗},u ^{∗}).
Figure 2 shows the plots of RMSE and MI for the proposed algorithm, bilateral filtering, and NLPCA. For the proposed algorithm, we have plotted these values for 10 logarithmicallyspaced values of λ in the range [0.01,1], which we found to give the best denoising results. For bilateral filtering, following [14], we have plotted these values for 10 linearlyspaced values of σ in the range [0.5,3.2]. From these plots it is clear that the proposed algorithm has achieved significantly better denoising results than bilateral filtering and NLPCA. Best results with the proposed algorithm are achieved with λ values around 0.1 and the denoising is too strong for λ>1. For bilateral filtering, we found that best denoising results were usually obtained for values of σ close to 3.0 and the performance did not improve or slightly deteriorated when σ was increased beyond 3.2. The solid squares on these plots show the optimum value of the corresponding parameter (i.e., lowest RMSE or highest MI). The phantom profiles shown in Fig. 1 for the proposed algorithm and bilateral filtering were obtained with the parameter values that resulted in the lowest RMSE.
Real conebeam projections of a physical phantom
Conebeam projections were acquired from a physical phantom using a Gamma Medica eXplore CT 120 microCT scanner. The scanner had a flat panel detector. The distance from the source to the detector and to the axis of rotation were 449.29 m m and 397.04 m m, respectively. We generated two scans of the same phantom:

1.
Lownoise scan. Consisting of 720 projections between 0° and 360° at 0.5° intervals. The tube voltage, tube current, and exposure time were 70 kV, 40 mA, and 25 ms, respectively.

2.
Highnoise scan. Consisting of 240 projections between 0° and 360° at 1.5° intervals. The tube voltage, tube current, and exposure time (per projection) were 50 kV, 32 mA, and 16 mAs, respectively.
Since we do not have access to the true (i.e., noisefree) projections, we compare the performance of the denoising algorithms in terms of the quality of the reconstructed images. We used the lownoise scan to reconstruct a highquality image of the phantom using the FeldkampDavisKress (FDK) algorithm [34]. We will refer to this as “the reference image”. To evaluate the denoising algorithms, we applied them on the highnoise projections, reconstructed the image of the phantom from the denoised projections using the FDK algorithm, and compared the reconstructed image with the reference image. Similar to the experiment with the simulated projections, we performed the denoising for 10 linearlyspaced values of σ in the range [0.5,3.2] for bilateral filtering. Similarly, we ran the proposed algorithm with 10 logarithmicallyspaced values of λ in the range [0.001,0.1]. Each projection was 875×568 pixels in size and the size of the reconstructed image of the phantom was 880×880×650 voxels, with isotropic voxels of 0.1×0.1×0.1 m m ^{3}.
In order to assess the overall quality of the reconstructed images, we use the following two criteria:

(a)
Root Mean Square of the Error (RMSE). Where error is defined as the difference between the image reconstructed from denoised projections and the reference image.

(b)
Structural similarity index (SSIM) between the reconstructed image and the reference image. SSIM is used as a measure of the overall closeness of two images and is defined as [35]:
$$ \text{SSIM}(x,\hat{x})= \frac{\left(2\mu_{x} \mu_{\hat{x}}+C_{1}\right)\left(2\sigma_{x \hat{x}}+C_{2}\right)}{\left({\mu_{x}^{2}}+\mu_{\hat{x}}^{2}+C_{1}\right)+\left({\sigma_{x}^{2}}+\sigma_{\hat{x}}^{2}+C_{2}\right)} $$where μ _{ x } and σ _{ x } represent the mean and standard deviation of the image x, \(\sigma _{x\hat {x}}\) is its covariance with image \(\hat {x}\), and C _{1} and C _{2} are constants.
The plots of RMSE and SSIM are shown in Fig. 3. Compared with both bilateral filtering and NLPCA, the image reconstructed from projections denoised using the proposed algorithm has a significantly lower RSME and higher SSIM. Best results in terms of SSIM with the proposed algorithm are obtained with λ=0.0129 and for bilateral filtering algorithm with σ=2.6.
The plots in Fig. 3 reflect the overall closeness of the images reconstructed from the denoised projections and the reference image. The imaged phantom had different modules that allowed for a more detailed evaluation of the quality of the reconstructed images [36]. A set of fine coils inside the phantom allow for visual inspection of the spatial resolution in the reconstructed image. Figure 4 shows two of these coils in the images reconstructed from noisy and denoised projections. These coils have thicknesses of 500 μ m and 200 μ m, corresponding to spatial resolutions of 1 and 2.5 line pairs per mm, respectively. The image shown for the proposed algorithm corresponds to λ=0.0129 and the image shown for bilateral filtering corresponds to σ=2.6. As we mentioned above, these parameter values led to highest SSIM. The images show a marked improvement in the image quality via sinogram denoising. It also seems that the proposed algorithm leads to a smoother image without affecting the spatial resolution. In Fig. 5 we have shown a profile through the center of the 500−μ m coil for the images reconstructed from noisy and denoised projections and also the difference between them and the reference image for a closer comparison. It is clear from these profiles that the image reconstructed from the projections denoised using the proposed algorithm are closer to the reference image.
In order to compare the denoising algorithms in terms of the tradeoff between noise and spatial resolution, we followed an approach similar to that suggested in [14]. Specifically, we computed the following two numbers as measures of spatial resolution and noise level in the reconstructed image of the phantom: Measure of spatial resolution. The imaged phantom included a slanted edge that consisted of a plasticair boundary, specially designed for accurate estimation of the modulation transfer function (MTF). We used the method proposed in [37] for estimation of the MTF over the range of spatial frequencies between 0 and 5 m m ^{−1}. As also suggested in [14], we use the spatial frequency at which the normalized MTF reaches a value of 0.10 as a representative number for spatial resolution. Measure of noise level. A uniform polycarbonate disk is included in the phantom for the purpose of assessing the noise level in the reconstructed image. We selected five cubes, each 10×10×10 voxels, at different locations within this disk and computed the standard deviation of the voxel values in each cube. We use the average standard deviation of voxel values in these cubes as a measure of noise level.
We computed the above two values for bilateral filtering algorithm with 10 linearlyspaced values of σ in the range [0.5,3.2] and for the proposed TVbased algorithm with 10 logarithmicallyspaced values of λ in the range [0.001,0.1]. In Fig. 6, we have shown plots of these two values for the three denoising algorithms. Note that a high spatial resolution and a low noise level are desirable. Therefore, all three denoising algorithms have improved the quality of the reconstructed image for the range of parameter values used (except for λ=0.1 with the proposed algorithm). Moreover, the proposed algorithm has achieved better results than bilateral filtering and NLPCA. Specifically, for λ∈[0.0077,0.0359] the proposed algorithm has achieved both higher spatial resolution and lower noise than bilateral filtering (for any parameter value) and NLPCA. In Fig. 6, we have also shown plots of the MTF obtained with the three denoising algorithms. All three sinogram denoising algorithms have led to an improvement in the spatial resolution in the reconstructed image. The proposed algorithm has resulted in a higher MTF than bilateral filtering and NLPCA for all spatial frequencies.
Real conebeam projections of a rat
We obtained a fresh rat carcass from our institutional animal facility and used the same microCT scanner described above to obtain post mortem images of the rat thorax. Since the animal was obtained post mortem, no ethical approvals were required for the microCT scans. Because the internal organs of the rat constantly moved, we were unable to create two identical scans with different noise levels as we did for the phantom. Therefore, we scanned the rat only once. The scan consisted of 720 projections between 0° and 360° at 0.5° intervals with the tube voltage, tube current, and exposure time equal to 70 kV, 32 mA, and 16 ms, respectively.
Similar to our approach in the experiment with the physical phantom, since we do not have access to the true projections, we evaluate the performance of the denoising algorithms in terms of the quality of the reconstructed images. To create a highquality reference image from the full set of 720 projections, we first reconstructed an initial image using the FDK algorithm. Then, we used five iterations of MFISTA algorithm [38] to improve the quality of the FDKreconstructed image. The resulting image had a very high quality and we used it as the reference image for evaluating the performance of the denoising algorithms. We applied the denoising algorithms on a subset of 240 projections of the same scan (projections at 1.5° intervals) and reconstructed the image of the rat using the FDK algorithm.
Similar to the physical phantom experiment, we use RMSE and SSIM as a measure of the overall closeness of the reconstructed images to the reference image. Figure 7 shows these criteria for the three sinogram denoising algorithms. From this figure, denoising of the projections with the proposed algorithm has lead to superior results in terms of RMSE and SSIM compared to bilateral filtering and NLPCA.
For a visual comparison, Fig. 8 shows a typical 2D slice of the reconstructed images of the rat. For the proposed algorithm and bilateral filtering, the images shown in this figure are obtained using the parameter values that resulted in the lowest SSIM, i.e., λ=0.0129 and σ=2.6 (see Fig. 7). The window of the linear attenuation coefficient, μ, used to display these slices was [0,0.55]. To allow a better visual comparison, we have selected two regions of interest (ROI) within this slice and have shown them in zoomedin views and with narrower μwindows. The ROI shown on the top left of each slice contains fat surrounded with soft tissue; this ROI is shown with a magnification factor of 1.5 and with a μwindow of [0.15,0.20]. The ROI shown on the top right of each slice contains bone surrounded with soft tissue; this ROI is shown with a magnification factor of 2.0 and with a μwindow of [0.18,0.50]. These images show a strong positive effect for sinogram denoising in terms of the visual quality of the reconstructed image. Moreover, denoising with the proposed algorithm seems to have resulted in a higherquality image, especially in the softtissue ROI.
In order to compare the denoising algorithms in terms of the tradeoff between noise suppression and spatial resolution, we again followed an approach similar to that proposed in [14]. Specifically, we selected an ROI shown in Fig. 9 and computed the following measures of noise level and spatial resolution: Measure of spatial resolution. We compute the maximum absolute value of the gradient (i.e., slope) along the line L marked in the ROI shown in Fig. 9 as a measure of spatial resolution. A sharper slope indicates a higher spatial resolution and it will give a larger gradient. Measure of noise level. We consider a cube of size 50×50×50 voxels, the crosssection of which is shown in the displayed ROI. From the reference image, we identified this cube as being highly uniform. Therefore, we computed the standard deviation of the voxel values in this cube as a measure of noise level.
The results are plotted in Fig. 9. This plot is very similar to the plot shown for the physical phantom experiment in Fig. 6. The main observations are that all three sinogram denoising algorithms have improved the quality of the reconstructed image in terms of spatial resolution and noise level, and that the proposed algorithm can outperform the bilateral filtering algorithm and NLPCA with the right selection of the regularization parameter. Specifically, with λ∈[0.0129,0.0359], the proposed algorithm has resulted in lower noise and better spatial resolution than bilateral filtering (with any choice of σ) and NLPCA.
Computation time
In order to compare the computational time of the proposed algorithm with that of bilateral filtering and NLPCA, we considered the denoising of 240 projections of the rat. As we mentioned above, each projection was 875×568 pixels. The proposed TVbased algorithm implemented in Matlab version R2012b and executed on a Windows 7 PC with 16 GB of memory and 3.4 GHz Intel Core i7 CPU needed approximately 6 minutes to denoise all 240 projections. In comparison, bilateral filtering and NLPCA needed 8.5 minutes and 42 minutes, respectively, for the same denoising task. In general, patchbased denoising methods such as NLPCA require long computational times. Nevertheless, sinogram denoising methods are in general much less computationally intensive than iterative reconstruction methods. A single forward or back projection using fast algorithms such as the separable footprints algorithm [39] or the distancedriven algorithm [40] takes approximately 2 hours on the same computer; each iteration of an iterative reconstruction algorithm needs one forwardprojection and one backprojection. Of course, it is becoming very common to implement iterative reconstruction methods on GPU, but the same can also be done for sinogram denoising algorithms. Sinogram denoising is by nature highly parallelizable because each projection can be denoised independently of others.
Conclusions
Sinogram denoising can be a very effective approach to improving the quality of lowdose CT images. In this paper, we presented a fast and efficient method for denoising of lowdose CT projections. The performance of the proposed method on simulated and real CT projections shows that this method can lead to a substantial reduction in the noise level without degrading the spatial resolution. Our results show that the proposed algorithm is superior to bilateral filtering and NLPCA in terms of denoising performance and computational speed. The proposed method is based on the assumption that the true image has a sparse gradient. Although algorithms based on this assumption have proved successful in image processing, CT projections are likely to not fit this model very well. This is because projections of piecewise constant objects are not in general piecewise constant, rather they are piecewise smooth. It is wellknown that TVdenoising of piecewise smooth images can result in staircase artifacts in the denoised image. Our approach of local denoising reduces the risk of staircase artifacts. Another approach to reducing these artifacts is to consider higherorder differentials in the model. Such a formulation may lead to results comparable to those reported in this paper, but the optimization algorithms involved will be significantly more complicated. We consider this approach as a future work.
Abbreviations
 CT:

Computed tomography
 PCA:

Principal component analysis
 FBP:

Filtered backprojections
 ROF:

RudinOsherFatemi
 TV:

Total variation
 NLPCA:

Nonlocal principal component analysis
 RMSE:

Root mean square of error
 MI:

Mutual information
 FDK:

FeldkampDavisKress
 SSIM:

Structural similarity index
 MTF:

Modulation transfer function
 ROI:

Region of interest
References
 1
Berrington de González A, Mahesh M, Kim K, et al.Projected cancer risks from computed tomographic scans performed in the United States in 2007. Arch Intern Med. 2009; 169(22):2071.
 2
Beister M, Kolditz D, Kalender WA. Iterative reconstruction methods in Xray CT. Physica Med. 2012; 28(2):94.
 3
Fessler JA. Statistical image reconstruction methods for transmission tomography. Handb Med Imaging. 2000; 2:1.
 4
Hsieh J. Adaptive streak artifact reduction in computed tomography resulting from excessive xray photon noise. Med Phys. 1998; 25(11):2139.
 5
Kachelriess M, Watzke O, Kalender W. Generalized multidimensional adaptive filtering for conventional and spiral singleslice, multislice, and conebeam CT. Med Phys. 2001; 28(4):475.
 6
Forthmann P, Köhler T, Begemann PGC, Defrise M. Penalized maximumlikelihood sinogram restoration for dual focal spot computed tomography. Phys Med Biol. 2007; 52(15):4513.
 7
LaRivière PJ. Penalizedlikelihood sinogram smoothing for lowdose CT. Med Phys. 2005; 32(6):1676.
 8
Kak AC, Slaney M. Principles of Computerized Tomographic Imaging. Philadelphia, PA: Society for Industrial and Applied Mathematics; 2001.
 9
Natterer F. The Mathematics of Computerized Tomography. Philadelphia, PA: Society for Industrial and Applied Mathematics; 2001.
 10
Jiao C, Wang D, Lu H, Zhang Z, Liang JZ. Multiscale noise reduction on lowdose CT sinogram by stationary wavelet transform. In: Nuclear Science Symposium Conference Record, 2008. NSS ’08. IEEE: 2008. p. 5339–44.
 11
Sahiner B, Yagle A. Reconstruction from projections under timefrequency constraints. Med Imaging IEEE Trans. 1995; 14(2):193.
 12
Macovski A. Medical Imaging Systems. Upper Saddle River, NJ: Prentice Hall; 1983.
 13
Wang J, Lu H, Liang Z, Eremina D, Zhang G, Wang S, et al.An experimental study on the noise properties of xray CT sinogram data in Radon space. Phys Med Biol. 2008; 53(12):3327.
 14
Manduca A, Yu L, Trzasko JD, Khaylova N, Kofler JM, McCollough CM, et al.Projection space denoising with bilateral filtering and CT noise modeling for dose reduction in CT. Med Phys. 2009; 36(11):4911.
 15
Wang J, Li T, Lu H, Liang Z. Penalized weighted leastsquares approach to sinogram noise reduction and image reconstruction for lowdose Xray computed tomography. Med Imaging IEEE Trans. 2006; 25(10):1272.
 16
Elad M, Aharon M. Image denoising via sparse and redundant representations over learned dictionaries. Image Process IEEE Trans. 2006; 15(12):3736.
 17
Mairal J, Elad M, Sapiro G. Sparse representation for color image restoration. IEEE Trans Image Process. 2008; 17(1):53.
 18
Chen Y, Yin X, Shi L, Shu H, Luo L, Coatrieux JL, et al. Improving abdomen tumor lowdose CT images using a fast dictionary learning based processing. Phys Med Biol. 2013; 58(16):5803.
 19
Li S, Fang L, Yin H. An efficient dictionary learning algorithm and its application to 3D medical image denoising. Biomed Eng IEEE Trans. 2012; 59(2):417.
 20
Buades A, Coll B, Morel JM. A nonlocal algorithm for image denoising. In: Computer Vision and Pattern Recognition, 2005. CVPR 2005. IEEE Computer Society Conference on, vol 2. IEEE: 2005. p. 60–5.
 21
Buades A, Coll B, Morel JM. A review of image denoising algorithms, with a new one. Multiscale Model Simul. 2005; 4(2):490.
 22
Maggioni M, Katkovnik V, Egiazarian K, Foi A. A nonlocal transformdomain filter for volumetric data denoising and reconstruction. IEEE Trans Image Process. 2013; 22(1):1057.
 23
Rudin LI, Osher S, Fatemi E, Physica D. Nonlinear total variation based noise removal algorithms. Nonlinear Phenomena. 1992; 60(14):259.
 24
Le T, Chartrand R, Asaki T. A variational approach to reconstructing images corrupted by Poisson noise. J Math Imaging Vis. 2007; 27(3):257.
 25
Sawatzky A, Brune C, Müller J, Burger M. Total variation processing of images with Poisson statistics In: Jiang X, Petkov N, editors. Computer Analysis of Images and Patterns, Lecture Notes in Computer Science, vol. 5702. Berlin, Heidelberg: Springer: 2009. p. 533–40.
 26
Chambolle A. An algorithm for total variation minimization and applications. J Math Imaging Vis. 2004; 20(12):89.
 27
Lefkimmiatis S, Bourquard A, Unser M. Hessianbased norm regularization for image restoration with biomedical applications. Image Process IEEE Trans. 2012; 21(3):983.
 28
Lysaker M, Tai X. Iterative image restoration combining total variation minimization and a secondorder functional. Int J Comput Vis. 2006; 66(1):5.
 29
Louchet C, Moisan L. Total variation as a local filter. SIAM J Imaging Sci. 2011; 4(2):651.
 30
Salmon J, Harmany Z, Deledalle CA, Willett R. Poisson noise reduction with nonlocal PCA. J Math Imaging Vis. 2014; 48(2):279. doi:10.1007/s1085101304356. http://dx.doi.org/10.1007/s1085101304356.
 31
Kundel HL, Van Metter RL, Beutel J. Handbook of Medical Imaging, Volume 1. Physics and Psychophysics. Bellingham, WA: SPIE Publications; 2000.
 32
Pluim J, Maintz J, Viergever M. Mutualinformationbased registration of medical images: a survey. Med Imaging IEEE Trans. 2003; 22(8):986.
 33
Bian J, Siewerdsen JH, Han X, Sidky EY, Prince JL, Pelizzari CA, et al.Evaluation of sparseview reconstruction from flatpaneldetector conebeam CT. Phys Med Biol. 2010; 55(22):6575.
 34
Feldkamp LA, Davis LC, Kress JW. Practical conebeam algorithm. J Opt Soc Am A. 1984; 1(6):612.
 35
Wang Z, Bovik A, Sheikh H, Simoncelli E. Image quality assessment: from error visibility to structural similarity. Image Process IEEE Trans. 2004; 13(4):600.
 36
Du LY, Umoh J, Nikolov HN, Pollmann SI, Lee TY, Holdsworth DW. A quality assurance phantom for the performance evaluation of volumetric microCT systems. Phys Med Biol. 2007; 52(23):7087.
 37
Buhr E, GüntherKohfahl S, Neitzel U. Simple method for modulation transfer function determination of digital imaging detectors from edge images. In: Medical Imaging 2003. Bellingham, WA, USA: International Society for Optics and Photonics: 2003. p. 877–84.
 38
Beck A, Teboulle M. Fast gradientbased algorithms for constrained total variation image denoising and deblurring problems. Image Process IEEE Trans. 2009; 18(11):2419.
 39
Long Y, Fessler J, Balter J. 3D forward and backprojection for xray CT using separable footprints. Med Imaging IEEE Trans. 2010; 29(11):1839.
 40
De Man B, Basu S. Distancedriven projection and backprojection in three dimensions. Phys Med Biol. 2004; 49(11):2463.
Acknowledgements
MicroCT imaging was performed in the Centre for HighThroughput Phenogenomics at the University of British Columbia, a facility supported by the Canada Foundation for Innovation, British Columbia Knowledge Development Foundation, and the UBC Faculty of Dentistry. Nancy L. Ford would like to acknowledge the funding from her NSERC Discovery Grant.
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
DK and PD conceived of the study. PD and NF performed the study design and data collection. Algorithm development and data analysis was carried out by DK and RW. DK prepared and submitted the manuscript. All authors cooperated in the interpretation of the results. All authors have read and approved the final manuscript.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 Computed tomography
 Sinogram denoising
 Poisson noise
 Total variation
 Lowdose CT