A Review on CT and X-Ray Images Denoising Methods

In medical imaging systems, denoising is one of the important image processing tasks. Automatic noise removal will improve the quality of diagnosis and requires careful treatment of obtained imagery. Computed tomography (CT) and X-Ray imaging systems use the X radiation to capture images and they are usually corrupted by noise following a Poisson distribution. Due to the importance of Poisson noise removal in medical imaging, there are many state-of-the-art methods that have been studied in the image processing literature. These include methods that are based on total variation (TV) regularization, wavelets, principal component analysis, machine learning etc. In this work, we will provide a review of the following important Poisson removal methods: the method based on the modified TV model, the adaptive TV method, the adaptive non-local total variation method, the method based on the higher-order natural image prior model, the Poisson reducing bilateral filter, the PURE-LET method


Introduction
Image denoising and noise removal with structure preservation is one of important tasks that are integrated in medical diagnostic imaging system, such as X-Ray, computed tomography (CT).X-ray and CT images are formed when an area under consideration of a patient is exposed under X-ray/CT and resulting attenuation is captured [1].The noise density in these systems follows by the Poisson distribution and well known as the Poisson noise, shot noise, photon noise, Schott noise or quantum noise.Although Poisson noise does not depend on temperature and frequency, it depends on photon counters.Poisson noise strength is proportional with the pixel intensity growth: Poisson noise at higher intensity pixel is greater than one at less intensity pixel [2].
Nowadays, digitization is an important technique to improve image quality in medical imaging systems and the Poisson noise characteristics needs to be considered to remove it effectively [1].Because the Poisson noise is a type of signal dependent noises, applying the usual denoising methods like for additive noises is ineffective, we need to design specific methods based on its characteristics.
There are many approaches were used to remove the Poisson noise, including total variation, mathematical transforms (wavelets, etc.), Markov random field, principal component analysis (PCA), machine learning etc.This paper mainly focuses on non-learning-based methods, learning technique is just a tiny part of this review that relates to the field of expert image prior model.
The approach that has been widely studied in the past few year and earn many achievements is regularization by total variation.This approach based on the regularization that was developed long time ago.Rudin et al. [3] used the total variation regularization to remove noise on digital images.Basically, they minimized an energy functional based on L 2 norm of image gradient with fixed constraint for noise variance.The proposed model was also known as ROF (Rudin-Osher-Fatemi) model.This work is wellknown and was cited by tens of thousands of times.However, the ROF model focuses on restoring images that are degraded by Gaussian noise.This model is ineffective to process Poisson noise: in the resulting image, the edge is not well preserved; if regularization strength is decreased, the noise in higher intensity-region still remains.
To pass over those limitations of the ROF model, Triet et al. [2] proposed an improved version that can process the Poisson noise well.This model is known as modified ROF model (MROF).However, both of original methods that based on ROF and MROF create an effect: artificial artifacts [1].The artificial artifacts on digital images are misrepresentations of image processing.This effect makes some regions of images get unnatural [4].The artifacts have many types, such as: staircasing, star, halo etc.In medical imaging, these artifacts can cause doctors to mistake for actual pathology.Usually, they need to learn to recognize these artifacts to avoid mistaking.So, during the processing, these artificial regions should not to be created.Prasath [1] proposed an adaptive version of MROF to remove this effect.This method is known as the adaptive total variation method (ATV).
A common problem of both MROF and ATV methods is ineffective to process on photon-limited image.To enhance quality of this type of image in denoising process, Salmon et al. [5] proposed the non-local PCA method.Thereafter, Liu et al. [6] proposed another adaptive nonlocal total variation method (ANLTV).This method increases the information structure of image and gives the better denoising result on photon-limited images.
Non-local approaches like ANLTV are state-of-theart.However, if the local models are combined with training process, we can get the result that is not inferior to other state-of-the-art non-local models.Wensen et al. [7] proposed a local variational model that incorporates the fields of expert prior image that is widely used in image prior and regularization model.This model is known as the higher-order natural image prior model (HNIPM).The HNIPM can remove Poisson noise on both high and low peak images.Although this model is local, since the model is trained on the Anscombe transform domain (very effective for Poisson denoising), it is also a competitive model to compare to other state-of-the-art Poisson denoising models.
However, above methods are performed on iteration and this requires more execution time to remove noise.Kirti et al. [8] proposed a spatial domain filter by modifying bilateral filter framework to remove Poisson noise.The Poisson reducing bilateral filter (PRBF) is non-iterative nature.So, it can treat Poisson noise faster than iterative based approaches.
Another approach is highly expectedwavelet and its modifications.Thierry et al. proposed a denoising method based on image-domain minimization of Poisson unbiased risk estimation: PURE-LET (Poisson Unbiased Risk Estimation -Linear Expansion of Thresholds) [9].This method is performed in a transformed domain: undecimated discrete wavelet transform and can be extended with some other transforms.Zhang et al. [10] also proposed a multiscale variance stabilizing transform (MS-VST) that can be deemed as an extension of Anscombe transform.This transform also can be combined with wavelet, ridgelet, and curvelet [10].Both PURE-LET and MS-VST are competitive relative to many existing denoising methods, in which, the VST based methods are new research trend for CT and X-Ray images denoising [11] [12] [13] [14], because of using VST, Poisson noise can be treated as the additive Gaussian noise.Hence, researchers can reuse the existing Gaussian denoising methods, that get many achievements and it is unnecessary to develop a partial denoising method to treat Poisson noise.
Our paper is organized as follows: in Section 2, a detail about image formation on CT/X-Ray imaging systems and characteristics of Poisson noise are provided; in Section 3, methodology of Poison denoising methods are covered shortly; Section 4 and Section 5 present the discussion about accuracy, performance, advantages/disadvantages of methods and the conclusion.

Image formation in medical imaging systems and Poisson noise
In CT and X-Ray imaging systems, to produce a radiographic image, X-Ray photons must pass through tissue and interact with an image receptor.The process of image formation is a result of differential absorption of the X-Ray beam as it interacts with the anatomic tissue [15].Differential absorption is a process whereby some of the X-Ray beam is absorbed in the tissue and some passes through the anatomic part.Because varying anatomic parts do not absorb the primary beam to the same degree, anatomic parts composed of bone absorb more X-Ray photons than parts filled with air.Differential absorption of the primary X-Ray beam creates an image that structurally represents the anatomic area of interest.Poisson noise is a fundamental form of uncertainty associated with the measurement of light, inherent to the quantized nature of light and the independence of photon detection [16].Its expected magnitude is signal-dependent and causes the dominant source of image noise except in low-light conditions.
Image sensors measure scene irradiance by counting the number of discrete photons incident on the sensor over a given time interval.Because of the photoelectric effect in digital sensors, photons are converted into electrons, whereas film-based sensors rely on photo-sensitive chemical reaction.Then, the random individual photon arrival leads to Poisson noise.
Individual photon detections can be considered as independent events that follow a random temporal distribution.The photon counting is a Poisson process, and the number of photons  measured by a given sensor element over a time interval  is described by the discrete probability distribution where expected number of photons per unit time interval, which is proportional to the incident scene irradiance.
Since the Poisson noise is derived from the nature of signal itself, it provides a lower bound on the uncertainty of measuring light.Any measurement would relate to Poisson noise, even under the ideal conditions of freenoise sources.When Poisson noise is the only significant source of uncertainty, as commonly occurs in bright photon-rich environments, imaging is called photon-limited [16].By the Poisson distribution, to reduce the Poisson noise, need to capture more photons.This requires longer exposures times or increasing the X-Ray intensity beam.However, the number of photons captured in a single shot is limited by the full well capacity of the sensor.Moreover, increasing exposures times or photon intensity beam would be harmful for health of patients.Since this limitation of technology, it is necessary to reduce the Poisson noise by image processing algorithms.
Figure 1 simulates the Poisson noise generation on image.We use the built-in imnoise function of MATLAB to generate the Poisson noise on skull image [17].The Poisson noise in the higher intensity regions is greater than one of the lower intensity regions.
3 Denoising methods on CT and X-Ray images

The modified ROF model
Suppose that a given grayscale image on Ω (a bounded open subset of ℝ 2 , i.e.Ω ⊂ ℝ 2 ), an expected denoising image that closely matches to observed image,  = ( 1 ,  2 ) ∈ Ωpixels.By using total variation regularization, Triet et al. convert the Poisson denoising problem to the following minimization problem: where,  > 0regularization parameter.
To solve this problem, Triet et al. used the gradient descent method that replaces the regularization parameter by function that is suitable to process noise on image regions with both low and high intensity.This manner exactly suits the signal-dependent nature of Poisson noise.

The adaptive Total variation method
The adaptive total variation method is similar with above method.However, the second term in ( 1) is replaced by an adaptive total variation: where, the Gaussian kernel for smoothing with  variance,  > 0contrast parameter, operator * is convolution.
In order avoid staircasing artifacts, Prasath [1] proposed the generalized inverse gradient term incorporating to the local statistics with patches extracted from image.The detail about this term is presented below.
Let  , be the local region centered at  with radius .Consider the local histogram of a pixel  ∈ Ω and its corresponding cumulative distribution function [18]: Where 0 ≤  ≤ , maximum possible pixel value of the image, |•|the number of elements of set (cardinality).The local histogram quantity to quantify local regions of given image is: Finally, the adaptive weight in (2) is defined as: The alternating direction method of multipliers [1] is provided to solve the problem (2).This iterative manner also gives good performance.

The adaptive non-local Total Variation method
In the case of photon-limited image, a lot of useful structure information of original image has been lost.So, the corrupted image is close to the binary image.If we only apply the denoising methods, such as the modified ROF model or the adaptive total variation, the denoising result is not really effective.For this type of images, firstly, we need to enhance image (improve light, contrast, etc.) and after that, perform the denoising process.
The method that Liu et al. [6] proposed is similar with above idea.For first step, they enhance the image detail by using Euler's elastica.In second step, they remove noise by using non-local total variation to aim to preserve the structure information.
The Poisson denoising model based on non-local total variation is provided as follows: where is non-local total variation, (, )the non-local weight to measure the similarity of patches centered at the pixels  and .The denoised version will be restored from (4) by using an inverse Anscombe transform as bellow: The alternating direction method of multipliers is also recommended to solve the models (3) and (4).

The higher-order natural image prior model
The denoising method by the higher-order natural image prior model is based on the fields of expert image prior model that can be presented as follows: where  number of filters,  set of learned linear filters with corresponding weights   > 0, number of image pixels, () = ln (1 +  2 )the potential function, (  * ) a convolution at pixel .The first term is derived from the fields of expert image prior model, the second term (, ) is data fidelity that has various forms.By using model ( 5), Wensen et al. [7] proposed two models that were trained in various transform domains: the first modelis trained in the original image domain with the Poisson noise statistics derived data term; the second modelis trained in the Anscombe transform domain with a quadratic data term.The first model removes Poisson noise on high peak images effectively, but it fails for low peak image.The reason is for the low peak image, there are large regions of image, in which, there are many pixels with zero intensity.This leads to those pixels with zero intensity cannot be updated and fixed at 0 in the iterations.Hence, noise still remains.The second model is powerful to remove noise for low peak images, but the quadratic data term is only effective to treat Gaussian noise.So, Wensen et al. combined the advantages of two models to make a novel model by replacing the quadratic data term in the second model by the Poisson noise statistics of data term in the first model.The resulting model proved its own power to remove the Poisson noise for both cases of high and low peak images.
The iPiano algorithm [19] is recommended to solve the resulting model.It is an efficient algorithm for nonconvex optimization problems.

The Poisson reducing bilateral filter
The bilateral filter was proposed by Tomasi et al. [20] to reduce additive Gaussian noise.This filter was developed based on the geometric and photometric distances in a local window.Kirti et al. [8] modified this filter by replacing the geometric distance by Poisson distribution.Therefore, the mean value is selected as mean of image intensity in a local window.For every mean value in the local window, the expected value is estimated by the maximum likelihood estimation method.
Since the Poisson reducing bilateral filter is non-iterative nature, its performance primarily depends on the maximum likelihood estimation method.

The PURE-LET method
The PURE-LET (Poisson Unbiased Risk Estimation -Linear Expansion of Thresholds) method [9] is extended from SURE-LET method [21].The PURE-LET method is used to reduce Poisson noise.Basically, this denoising method was proposed based on a minimization of Poisson unbiased risk estimation by using the linear expansion of thresholds (LET).Luisier et al. [9] proposed the PURE-LET to reduce Poisson noise without any priori hypotheses on noise-free image.
The main goal of this proposed denoising method is a minimization of the mean squared error of the noise-free image and the denoised image.However, since the noisefree image is unknown, unbiased risk estimation was used that known as the Poisson unbiased risk estimation.This estimation was given in the unnormalized-Haar-discretewavelet-transform domain.In this estimation, an unknown image function used to replace for the noise-free image.
To minimize this estimation, above unknow image function will be expressed in the linear expansion of thresholds.If elementary denoising functions are given, the minimization problem gets to be the problem of finding weight parameters in the linear expansion.Hence, the main task of this PURE-LET method focuses on solving a linear system of equations, in which the variables are weight parameters of the linear expansion.
The linear expansion of thresholds can be presented in transformed domain, such as unnormalized wavelet transform, Anscombe transform and Haar-Fisz transform.
Another important task in the PURE-LET methods is choosing a set of elementary denoising functions (or thresholding functions).These functions need to be satisfied the following minimal properties: differentiability, anti-symmetry, linear behavior for large coefficients.
The PURE-LET method is a competitive method to compare to other state-of-the-art Poisson denoising methods.This method is also easy to be extended to treat other noises, such as Gaussian noise [22] [23] [24] [25], the mixed noise [26] [27] [28] [29] [30].The method performance much depends on the performance of methods of solving linear system of equations, for example, the Gauss-Seidel method.

The multiscale variance stabilizing transform method
The multiscale variance stabilizing transform method is proposed by Zhang et al. [10] to reduce Poisson noise on photon-limited image.This method is based on the variance stabilizing transform (VST) that is incorporated within the multiscale framework offered by the undecimated wavelet transform (UWT).This transform is used because of its translation-invariant denoising.The denoising task comes to finding coefficients of the multiscale variance stabilizing transform.By using these coefficients, we can estimate the noise-free image.
The denoising method involves in the following steps: transformationcomputation of UWT in conjunction with MS-VST; detectiondetection of significant detail coefficients by hypotheses test; estimationreconstruction of the final estimate by using the knowledge of the detected coefficients.Since the signal reconstruction requires inverting the MS-VST-combined UWT, this reconstruction process is formulated as a convex sparsity-promotion optimization problem.This optimization problem can be solved by many iterative methods, such as the iterative hybrid steepest descent method.
The MS-VST method can be combined with wavelet, as well as ridgelet (wavelet analysis in Radon domain) or curvelet.Further, this method can also to be extended to reduce other types of noise.

Adaptive variance stabilizing transform based methods
The Poisson denoising methods by VST-based approach is often performed by three steps: applying the variance stabilizing transform, such as Anscombe transform; applying the denoising methods to resulting image, in which the denoising methods are the one for additive Gaussian noise; using inverse transformation to denoised image to get the Poisson denoised image.Hence, VST-based methods can use state-of-the-art Gaussian denoising methods.By this idea, there are some very effective methods, such as BM3D [31], SAFIR [32], BLS-GSM [33].
For VST-based methods, the choice of inverse transformation is very important.Makitalo and Foi [11] proposed the optimal inverse Anscombe transform.The adaptive variance stabilizing transform-based method of Makitalo et al. can be covered as follows: Step 1: Apply the Anscombe transform to Poisson noisy image to get asymptotically additive Gaussian noisy image.For the observed pixel values obtained through an image acquisition device, the Anscombe transform is ,  = ( 1 , … ,   ),  − pixel numbers.
Step 2: Denoise the transformed images by additive Gaussian denoising method.
Step 3: The denoising of () produces a signal  that considered as an estimate of {()},  = ( 1 , … ,   )pixel values of denoising image, {.}the mean.So, it is necessary to apply inverse transformation to  to obtain the desired estimate of .The inverse transformations can be used include: a) The exact Unbiased inverse b) The ML inverse c) The MMSE inverse where, (|) = Another adaptive VST-based method that has high accuracy and performance to treat Poisson noise was proposed by Azzari and Foi [12].This method is known as the iterative VST-based method.
This method is also handled via three steps as above.However, in step 2, authors proposed another method to remove noise, but they did not use existing additive Gaussian denoising methods.The method is effective and has high performance because it exploited characteristics of Anscombe transformation.
This process loops until  = .
The accuracy and performance of this method are competitive to other state-of-the-art Poisson denoising methods.

Other Poisson denoising methods
Since the Poisson denoising problem has important role not only in medicine, but also in other fields, such material science, astronomy etc., beside above state-of-the-art denoising methods, there are also many other denoising methods are highly assessed, such as: The adaptive BLS-GSM method of Li et al. [34].They proposed this method based on Bayesian least squares method.Basically, this Poisson denoising method is a term of VST-based approach.
The optimized anisotropic Poisson denoising method of Radow et al. [35].This method is proposed based on variational approach and anisotropic regulariser in the spirit of anisotropic diffusion.This method can be considered as a part of the total variation regularization.
The Poisson denoising based on greedy approach of Dupe and Anthoine [36].The goal of this method is combination of a greedy method with Moreau-Yosida regularization of the Poisson likelihood.
The Poisson reduction based on region classification and response median filtering of Kirti et al [37].Their contribution is usage of modified Harris corner point detector to predict noisy pixels and responsive median filtering in spatial domain.
The primal-dual hybrid gradient algorithm [38] is a Poisson denoising method that should be also noticed.This method is based on total variation regularization and primal-dual hybrid gradient.So, this method has very good performance.

Discussion
Firstly, we will discuss on the accuracy of Poisson denoising methods.The MROF, ATV, ANLTV and HNIPM methods based on regularization, their accuracy is good enough to perform in medical imaging systems.Since the HNIPM method is trained on Anscombe transform domain, regardless of its localization, its accuracy is competitive enough to other methods.If we combine the MROF, ATV, ANLTV methods with training process to select optimal parameters in iterative manners, their accuracy might be so far better than the HNIPM method, especially, for the ANLTV method, because it does not change the information structure of image in denoising process.
An effect that reduces the accuracy in denoising process is artificial artifacts.Almost of local usually create this effect.So, we need to perform some techniques to avoid adding artifacts to images, such in the case of the ATV method.For non-local methods, since the information structure of image is preserved, the artifacts will be seldom added.The PRBF method has the lowest accuracy to compare to other denoising methods, including the PURE-LET, MS-VST and adaptive VST-based methods.When filter noise by PRBF, the hallo artifacts will appear in resulting images and the artifacts strength depends on filter parameters.Although we can control these parameters to reduce the hallo artifacts, it is very hard to select optimal values.There are some methods were developed to reduce this type of artifacts [39] [40], but it is still unfinished, especially, on Poisson noise reduction process by bilateral filter.For the PURE-LET, MS-VST and adaptive VST-based methods, the accuracy might be better than local variational based methods without training process, particularly, for the photon-limited images.However, the PURE-LET method is usually unstable.In our test, the denoising result by the PURE-LET method is slightly different in every execution, regardless of unchangeable input setting of parameters and configuration.When we compare the MS-VST method to the PURE-LET method, the MS-VST method has better accuracy, especially, for photon-limited images [10].The adaptive VST-based methods have better accuracy and performance to compare to MS-VST method [11] [12].Both the PURE-LET and MS-VST cause the artifacts.For the adaptive VST-based methods, appearance of the artifacts depends on selection of Gaussian denoising methods.
Secondly, we focus on method performance by assessing the execution time.Poisson denoising methods, such as MROF, ATV, ANLTV, PURE-LET, MS-VST and adaptive VST-based methods are designed on iterative manner, so their execution time is longer than one of the PRBF method.The PRBF method is very fast and this is proven in processing large images.Execution time of both of PURE-LET, MS-VST and adaptive VST-based methods also depends on computation time of transforms.Otherwise, for the PURE-LET method, it also depends on execution time of solving system of linear equations, and for the MS-VST methoddepends on performance of method to solve convex optimization problem, such as the hybrid steepest decent method, and for adaptive VST-based methodsdepends on performance of selective Gaussian denoising methods.For other methods: MROF, ATV, ANLTV, HNIPM, execution time much depend on performance of method to solve optimization problem (convex optimization for the MROF, ATV, ANLTV methods and nonconvex optimization for the HNIPM method).There were some methods are recommended in their proposed works to solve these optimization problems: the gradient descent method, the alternating direction method of multipliers for convex optimization; iPiano for non-convex optimization.However, for the convex optimization, we can use other faster methods, such as: the primal-dual modified extragradient method, the primal-dual Arrow-Hurwitz method, the graph-cut method [41].In work [41], Chambolle et al. showed comparison of execution time of above methods with the alternating direction method of multipliers.Among of these methods, the primal-dual Arrow-Hurwitz method is the fastest, but proof of its convergence is open problem.The primal-dual modified extragradient method is certainly convergent and it is easy to parallelize on GPU.The graph-cut method is very fast and give exact discrete solutions, but an efficient parallelization on GPU is still open problem.For the non-convex optimization, the iPiano method is state-of-the-art algorithm and fast enough to applied in this situation.Parallelization of the iPiano method is still open problem.Hence, the execution time problem of all above methods can be solved by combining with higher performance algorithms and/or parallel processing.
Finally, about methodology, the MROF, ATV, ANLTV and PRBF methods are simple and easy to understand and easy to write program.The HNIPM is slightly more complex and requires the training process.Both of PURE-LET, MS-VST and adaptive VST-based methods are the most complex.They are performed in various transform domains.Their accuracy and performance also depend on calculation of these transforms.
To choose suitable method for Poisson denoising in specific cases, we need to know their advantages and disadvantages.These advantages and disadvantages are listed in Table 1 in terms of denoising capabilities of the reviewed denoising methods here.After decades of denoising research there are no universal denoising method even in the case of additive Gaussian noise.However, by concentrating on the state of the art denoising methods with emphasize of domain specific techniques will pave the way for choosing an optimal denoising method.We believe the overview of Poisson denoising methods based on mathematically well-defined techniques studied here can be used by researchers in developing and utilizing these in various domains.

Conclusion
The denoising on CT/X-Ray images is still a challenge in medical image processing, especially, on the photon-limited images.The state-of-the-art methods cannot solve simultaneously the following tasks: high accuracy on both photon-limited and photon-unlimited images, avoid adding artificial artifacts and the performance.The goal to develop an effective universal method that reduces multiple types of noise is even more difficult challenge.
In this paper, we reviewed on the following methods: MROF, ATV, ANLTV, HNIPM, PRBF, PURE-LET and MS-VST.The PRBF is excellent choice if the execution time is the most important.However, if the accuracy is priority, non-local methods are recommended.If we need to process the photon-limited images, the ANLTV, MS-VST and adaptive VST-based methods are very good choices.If we want to exploit the existing Gaussian denoising methods, we can use adaptive VST-based methods, including MS-VST.
During denoising process is performed, it is necessary to avoid adding artificial structures, and one can choose ATV or ANLTV methods that provide good denoising performance without introducing discernible artifacts.In this case, the VST-based methods can be used if they are combined to the image structure preservation Gaussian denoising methods, such as BM3D [31], SAFIR [32] etc.
By the research trend, the VST-based approach is a novel option by the criteria to create an "universal" method to remove multiple type of noises.This approach has potential if it is possible to expand the VST-based approach to apply to other signal dependent noises.

a b Figure 1 :
The Poisson noise generation: a) The expected noise-free image; b) The noisy image