ASSESSMENT OF RESTORATION METHODS OF X-RAY IMAGES WITH EMPHASIS ON MEDICAL PHOTOGRAMMETRIC USAGE

Nowadays, various medical X-ray imaging methods such as digital radiography, computed tomography and fluoroscopy are used as important tools in diagnostic and operative processes especially in the computer and robotic assisted surgeries. The procedures of extracting information from these images require appropriate deblurring and denoising processes on the preand intra-operative images in order to obtain more accurate information. This issue becomes more considerable when the X-ray images are planned to be employed in the photogrammetric processes for 3D reconstruction from multi-view X-ray images since, accurate data should be extracted from images for 3D modelling and the quality of X-ray images affects directly on the results of the algorithms. For restoration of X-ray images, it is essential to consider the nature and characteristics of these kinds of images. X-ray images exhibit severe quantum noise due to limited X-ray photons involved. The assumptions of Gaussian modelling are not appropriate for photon-limited images such as X-ray images, because of the nature of signal-dependant quantum noise. These images are generally modelled by Poisson distribution which is the most common model for low-intensity imaging. In this paper, existing methods are evaluated. For this purpose, after demonstrating the properties of medical X-ray images, the more efficient and recommended methods for restoration of X-ray images would be described and assessed. After explaining these approaches, they are implemented on samples from different kinds of X-ray images. By considering the results, it is concluded that using PURE-LET, provides more effective and efficient denoising than other examined methods in this research. * Corresponding author


INTRODUCTION
There are different types of medical X-ray images such as digital radiography, computed radiography, computed tomography and fluoroscopy which are widely applied for operating planning, diagnosis and treatment processes.For achieving more reliable results in medical applications, the captured images should be denoised by considering the nature of noises.This issue is very important for photogrammetry usage and 3D reconstruction applications from X-ray images because of the limited number of X-ray photons.Therefore, for achieving more accurate results, in this research, algorithms for reducing the noise will be assessed.It is often suggested to approximate the combined effect of the various noise sources by a Gaussian distribution.The variance of a Gaussian noise is constant.In image processing, a common assumption is that the noise is not only Gaussian, but additive white Gaussian noise (AWGN) (Mäkitalo, 2013).AWGN is signal-independent and normally distributed.Since this method completely ignores the signal-dependent noise, it is not sufficient and accurate for X-ray medical images.In order to denoise and restore degraded medical X-ray images adequately, the nature and characteristics of images should be considered.Noise in X-ray imaging comes from various sources such as quantum noise, electronic noise, sampling noise, and anatomical noise.X-ray imaging operates in the photon limited realm due to the low photon counts available over a reasonable exposure time.X-ray images exhibit severe quantum noise due to limited X-ray photons involved.The most important noise in X-ray images, quantum noise, is a signaldependant noise.Quantum noise is produced by the random manner of production and interaction of the photons in a detector.Gaussian modelling assumptions are accurate when the sources of error are signal-independent.Therefore, they are not suitable for photon-limited images such as X-ray images, because of the nature of signal-dependant quantum noise.X-ray images are generally modelled by Poisson distribution which successfully models the photon counting statistics of imaging detectors and it is the most common model for low-intensity imaging.Since the noise variance depends on the intensity value, Poisson noise is signal dependent.Despite of similar nature and properties of noise in X-ray images, various types of X-ray images such as fluoroscopy, digital radiography, film-screen radiography, and computed tomography have different characteristics and levels of noise.Therefore, the restoration methods should be applied by considering their noise behaviour.In fluoroscopy, there is a continuous beam of radiation.The relatively low exposure in fluoroscopy produces images with high quantum noise, due to the more restricted number of X-ray photons.In normal fluoroscopy, an average of several frames, not one frame, is seen at a time.Digital radiography such as direct radiography and computed radiography has lower quantum noise than fluoroscopy.In Digital radiography there is not a fixed sensitivity for receptors, in contrast to film-screen receptors.Digital receptors have a wide exposure dynamic range.For computed radiography, it was shown that a coherent data distribution where the noise variance increases linearly with the image intensity which is the behavior of images with Poisson noise.In film-screen radiography, the film must first be scanned before noise can be studied on it and the scanning process introduces photon, read-out and quantization noise in addition to grain noise already present in the film (Gravel et al, 2004).(Gravel et al, 2004)expressed that the scanner quantization noise dominates over all other types of noise in the film-based radiography and the photon noise which is present during the X-ray exposure and recorded in the grain distribution is negligible when compared to the scanning noise.In Computed tomography (CT), Gaussian rather than Poisson distributions of local pixel intensities on real and synthetic CT scans series are observed.The image reconstruction algorithms involve weighted linear combinations of intensities that transform the noise PDF from Poisson to Gaussian according to the central limit theorem.The noise PDF on such an image has an asymptotic Gaussian distribution.Therefore, they are considered Gaussian despite of their signal-dependent variance.By considering these issues, it seems essential to evaluate denoising algorithms and their effects on each type of X-ray images.Although many researches have been done on methods of denoising images with Gaussian noise and different algorithms have been introduced for this purpose, there are relatively a few methods offered for restoring images properly with Poisson noise, such as radiographs.Denoising of X-ray images should be considered as an important pre-processing step in 3D reconstruction, more specifically to facilitate and increase the accuracy of the segmentation process which has high importance for registration step.Hence, the existing denoising methods should be evaluated in order to find appropriate methods that improve the results of image processing procedures.In this paper, after demonstrating the properties of medical X-ray images, existing methods for restoration of X-ray images will be described and assessed.First, various approaches such as methods based on Variance Stabilising transformation, Total Variation and Tikhonov regularization with different data fit functions, methods, PDE and complex diffusion processes, Block-matching and 3D filtering methods, and PURE-LET are explained.Next, these methods will be applied and implemented on image samples from three kinds of X-ray images which are radiography, CT, and fluoroscopic images with different levels of noise.Finally, the results of the explained methods are compared.For this purpose, different metrics such as PSNR (peak signal to noise ratio) are explored and applied for the assessment.

METHODS
As mentioned, X-ray images contains high quantum noise because of their nature.Quantum noise is a signal dependent noise, in contrast to Gaussian noise which is signal independent.In images with signal-dependent noise, the noise variance is typically not constant and varies with the expectation of the pixel value.Therefore, they are generally modelled by Poisson distribution.Denoising of these images can be done by considering the noise behaviour and statistics directly such as PURE-LET method or indirectly such as variance-stabilizing transform based methods.In the following, various important methods for this issue are explored.We will review several outperforming denoising algorithms (VST-based, BM3D, PURE-LET and Regularization and gradient-based methods).It should be mentioned that BM3D is a strong denoising algorithm for images with Gaussian noise but it is widely used in combination of other methods for denoising images with signal-dependent noise.

Variance Stabilizing Transformations
Many of the denoising algorithms are designed to reduce noise in transform domain.Poisson statistics are more difficult to handle in transformed domain.Hence, a solution for denoising images with Poisson noises is variance-stabilizing and rendering the noise variance constant throughout the image, which is called Gaussianizing, in order to remove the signaldependency.Gaussianizing the Poisson measurements is done by nonlinear mapping.This method was first presented by (Anscombe, 1948) and it was named as Anscombe variancestabilizing transform (VST).(Donoho, 1993) applied this method in denoising problems for the first time.Today, Anscombe transformation is widely used in denoising applications.In variance-stabilizing transform (VST) denoising procedure, first a nonlinear variance stabilizing transformation is applied on the noisy data.The achieved transformed data can be considered to have an approximately Gaussian noise distribution with a constant variance.After the stabilization step, any denoising algorithm designed for the removal of Gaussian noise can be used.Finally, the estimate of the unknown noise-free image is obtained by applying an inverse VST to the denoised data.There are other variance-stabilizing transformations such as Anscombe taransformation (Anscombe, 1948), GAT (Murtagh et al, 1995), MS-VST (Zhang et al, 2008), OPT-VST (Foi, 2008), and the exact unbiased inverse (Mäkitalo et al, 2011) of the Anscombe Taransformation.The classical Anscombe transformation has been the most common VST applied for the Poisson noise because of its simplicity and low expense.However, some images are degraded by photon and readout noises.In this case, the reading noise from a CCD detector generated by the thermal fluctuations.Therefore, these images contain some amount of Gaussian noise, too.For denoising images corrupted by mixed Poisson-Gaussian noise, the generalized Anscombe transformation (GAT) was proposed by (Murtagh et al, 1995) in order to stabilize the variance of Poisson-Gaussian noise.MS-VST is another denoising algorithm which is proposed by (Zhang et al, 2008) as an extension of the Anscombe transform applied on a filtered discrete Poisson process, yielding a near Gaussian process with asymptotic constant variance (Zhang et al, 2008).This method uses of a well-designed VST allowing to Gaussianize and stabilize a filtered Poisson process.Then, it is combined with wavelets, ridgelets and curvelets yielding MS-VSTs.It has been shown that the MS-VST approach is effective in recovering important structures of various shapes within lowcount images.For the inverse of MS-VST which is needed to be applied on the denoised coefficients, the inversion is formulated as a convex sparsity-driven minimization problem, solved by iterative steepest descent.This method is efficient in even very low-count situations.Since exact stabilization and even achieving some approximate stabilization for many of the distributions such as Poisson is challenging, Foi (2008) approached the variance stabilization problem as an explicit optimization problem.In the proposed method, the discrepancy between the standard deviation of the transformed variables and a fixed desired constant is measured by a nonlinear stabilization functional, which is then iteratively minimized.For applying variance stabilization, a properly constructed forward VST and also a suitable inverse transformation are needed.However, all practical transformations are approximately or asymptotically stabilizing and normalizing (Mäkitalo, 2013).Therefore, (Mäkitalo et al, 2011) proposed an exact unbiased inverse of the Anscombe transformation that is equally applicable regardless of the values.The same work was done for the GAT.(Mäkitalo, 2013) showed the importance of a properly designed inverse transformation, and constructed exact unbiased inverses for these two VSTs.They proposed an exact unbiased inverse for the Anscombe transformation and similarly for the GAT.

PURE-LET
PURE-LET method was designed by (Luisier et. al., 2010) as a non-bayesian framework to estimate Poisson intensities in the unnormalizied Haar wavelet domain, particularly for Poisson noise.PURE-LET's concept is built on the idea of estimating the "risk" (MSE) between the noise-free and the denoised images.This algorithm is based on: (1) the minimization of an unbiased estimate of the MSE for Poisson noise, (2) a linear parametrization of the denoising process and, (3) the preservation of Poisson statistics across scales within the Haar DWT (Luisier et al, 2010).This method is similar to SURE-LET approach for Gaussian denoising.For a Gaussian distribution, this MSE can be estimated with Stein's unbiased risk estimate (SURE).However, for the Poisson noise, Luisier et al derive an interscale Poisson unbiased risk estimate (PURE), which minimizes the MSE in the Haar wavelet domain.Then, they minimize the estimated MSE over several denoising processes, in order to find the one providing the best SNR.In particular, these denoising processes are wavelet estimators expressed as a linear expansion of thresholds (LET); thresholds in this context mean arbitrary elementary estimators with unknown weights (Mäkitalo, 2013).Thus, the PURE estimate can be minimized by solving a low-dimensional system of linear equations.The Poisson-Gaussian noise model is considered in (Luisier et al, 2011) and PURE-LET was extended as a Poisson denoising method.Specifically, they generalize their method from the Haar wavelet domain to a general redundant transform domain, and update the PURE estimate to take into account both the Poisson and Gaussian noise components (Mäkitalo, 2013).Finally, they show that PURE allows for the global optimization of a LET spanning several different transform domains, such as undecimated wavelet transform (UWT) and block DCT.This method has low computational complexity and memory.In this method, all of the parameters of the algorithm are adjusted completely automatically.Besides, this method has low computational complexity and memory.

Block-Matching and 3D Filtering
Block-Matching and 3D Filtering (BM3D) is a strong denoising method proposed for the problem of attenuating AWGN.However, due to its high quality performance and fast execution speed, it is applied widely for reducing even signal-dependent noise in combination with VST.Therefore, it is considered in this paper.BM3D procedure is based on an enhanced sparse representation in transform-domain.The enhancement of the sparsity is achieved by grouping similar 2D image patches into 3D groups (Dabov, et al, 2007).Collaborative filtering is a special procedure developed to deal with these 3D groups.BM3D was proposed as an extension of Non Local Means and wavelet domain transform filtering.This method exploits the intra-and inter-patch correlation resulting from smoothness and self-similarity of the images jointly.This is an edgepreserving denoising method.As shown in Figure 1, BM3D algorithm has two major filtering steps.In both stages collaborative filtering is utilized (Elahi et al, 2014).Collaborative filtering itself has four stages: 1) finding the image patches similar to a given image patch and grouping them in a 3D block, 2) 3D wavelet transformation of each stack of patches, 3) denoising the wavelet coefficients (thresholding or Wiener filtering) and, 4) inverse 3D transformation.In the first step, denoising would be done by hard thresholding, but, in the second step, BM3D denoises the patches by Wiener filter.The result is a 3D estimate that consists of the jointly filtered grouped image patches.By attenuating the noise, the collaborative filtering reveals even the finest details shared by the grouped patches at the same time, it preserves the essential unique features of each individual patch.The filtered patches are then returned to their original positions.Since these patches overlap, for each pixel we obtain many estimates which need to be combined.Aggregation is a particular averaging procedure used to take advantage of this redundancy.The first collaborative filtering step is much improved by a second step using Wiener filtering.This second step mimics the first step, with two differences.The first difference is that it compares the filtered patches instead of the original patches.The second difference is that the new 3D group (built with the unprocessed image samples, but using the patch distances of the filtered image) is processed by Wiener filtering instead of a mere threshold.The final aggregation step is identical to those of the first step (Lebrun, 2012).

Regularization and gradient-based methods
Denoising algorithms based on gradient dependent regularizers such as PDE-based processes, Tikhonov and total variation regularization methods are widely applied for denoising purposes.Therefore, in this paper some of the well known proposed methods in this group are explored and evaluated.Partial differential equations (PDEs) have been commonly used for denoising of images with Gaussian noise.PDE-based nonlinear diffusion processes originate from variational calculus (Aubert et al, 2002).The diffusion-like PDE is derived by a functional minimization process.The primary PDE-based methods had the disadvantage of smoothing edges.(Perona et al, 1990) introduced a PDE-based non-linear adaptive diffusion method, called anisotropic diffusion, in which an anisotropic diffusion coefficient decreases the smoothing effect near the edges. .(Gilboa et al, 2004) proposed an adaptive PDE filter which combines the Perona-Malik anisotropic diffusion and a complex shock filter, using a sharpness factor in the direction of the image gradient.To regularize the shock filter, a complex diffusion term was added and the imaginary value was used to control the direction of the flow instead of the second derivative.They improved the Perona and Malik method by applying the diffusion equations in the complex domain.They generalized the linear scale spaces in the complex domain, by combining the diffusion equation with the free Schrödinger equation (Gilboa et al, 2004).In this method, the imaginary value can serve as an edge detector (smoothed second derivative scaled by time) when the complex diffusion coefficient approaches the real axis (Gilboa et al, 2004).This method is often used for Gaussian noise, but it is applied by several authors such as (S.Kadoury et al, 2010) for denoising images with Poisson noise.Therefore, in this paper, the performance of this method is evaluated for Denoising images corrupted with Poisson noise.Regularization approaches are obtained by formulating the image restoration problem as a non-negatively constrained minimization problem: Min J (x) = J0 (x) + η JR (x) (1) s.t x≥0 η is the regularization parameter.The choice of a suitable fit function J0 (x) mainly depends on the statistics of the noise affecting the recorded data, while the choice of the regularization function JR (x) is related to the features of the image to be restored as it is shown in (Landi et al,2013).When the data are affected by Gaussian noise, the Least-Squares (LS) function can be the best statistical estimator.But, the generalized Kullback-Leibler divergence can be the optimal data fit term when the data are affected by Poisson noise.Tikhonov regularization function can be used when the image to be restored is made by diffused objects, such as certain images from microscopy.(Landi et al, 2012(Landi et al, , 2013) ) proposed a deblurring and denoising method for images with Poisson noise by minimizing the objective function which is the sum of the Kullback-Leibler divergence, used to express fidelity to the data in the presence of Poisson noise, and a Tikhonov and TV regularization term.Although the TV regularization is well known for recovering sharp edges of an image, it has been shown that the TV norm transforms a smooth signal into piecewise constants.In general, Tikhonov regularization has smoothing effect on the images, which is not appropriate for images which have edges.In contrast, the total variation (TV) based regularization technique preserves the edges in the restored image.Total variation (TV) regularization is commonly used in the denoising of medical images, such as X-ray images.

EXPERIMENTS
In this paper the performance of the denoising algorithms is evaluated by considering their results on test images from various types of X-ray medical images such as Digital Radiography, Computed Tomography (CT), and Fluoroscopy.For the first step, the Poisson noise is simulated and original images are corrupted with Poisson noise.Then, different denoising algorithms are implemented on them.The algorithms which are applied in this research are as followed: (1) PURE-LET method, (2) MS-VST, (3) Exact unbiased inverse of the VST applying BM3D which is designed by (Mäkitalo et al, 2012), (4) Total variation, and (5) a PDE-based method proposed by (Gilboa et al,2004).The results of implementing these denoising algorithms are shown in figures 2, 3,and 4 for various kinds of X-ray images.It should be mentioned that the noise in the CT scan image has a different behaviour from other mentioned images, and this issue is considered in the tests.The comparison between achieved results of these five methods is done by applying PSNR, and MSE metrics.This assessment shows that PURE-LET method works more efficiently in relatively all X-ray images with simulated Poisson noise.In addition, this algorithm adjust all the parameters automatically, without any need for user input.That is very important and efficient for real data.Also, this method takes less time to reach the acceptable results than VST-based methods.

DISCUSSION AND CONCLUSIONS
In this research, after exploring the properties of noise in radiographs and explaining characteristics of these images, more suggested denoising algorithms for images with Poisson noise have been described.In medical photogrammetry and 3D reconstruction applications from images, the edges acts important role, particularly for segmentation.Therefore, the edge preserving methods are considered to enhance and denoise images in a pre-processing step in these applications.In denoising methods designed for AWGN the signal-dependent noise is ignored.So they are not suitable and accurate for denoising X-ray images, because of their quantum noise.For choosing an appropriate denoising method in medical photogrammetric applications, the performance of the denoising algorithms designed for Poisson noise beside some strong denoising methods not specifically proposed for this kind of noise, has been evaluated.Poisson noise and low-count photon images.In addition, in this algorithm, it is not necessary to adjust the parameters by users, since they are configured automatically.This causes reproducibility and makes the method more applicable for real data.Also, this method needs computationally less time to reach the acceptable results comparing to VST-based methods.It has low computational complexity and memory, so it can be used for even large data sets.Therefore, it is suggested to be applied for denoising X-ray images and any other images with Poisson noise.
The algorithms which have been utilized for denoising in this research are PURE-LET, MS-VST, Exact unbiased inverse transformation applying BM3D, Total variation, and a PDE-based method.The results of these methods on sample images from various types of X-ray medical images such as Digital Radiography, CT, and Fluoroscopy which corrupted with simulated Poisson noise, indicate that PURE-LET, which specifically designed for Poisson noise, produces more efficient and reliable results, in most of the Xray images with simulated Poisson noise even with high