alexa Comparative Analysis between a Variational Method and Wavelet Method PURE-LET to Remove Poisson Noise Corrupting CT Images | OMICS International
ISSN: 0974-7230
Journal of Computer Science & Systems Biology
Like us on:
Make the best use of Scientific Research and information from our 700+ peer reviewed, Open Access Journals that operates with the help of 50,000+ Editorial Board Members and esteemed reviewers and 1000+ Scientific associations in Medical, Clinical, Pharmaceutical, Engineering, Technology and Management Fields.
Meet Inspiring Speakers and Experts at our 3000+ Global Conferenceseries Events with over 600+ Conferences, 1200+ Symposiums and 1200+ Workshops on
Medical, Pharma, Engineering, Science, Technology and Business

Comparative Analysis between a Variational Method and Wavelet Method PURE-LET to Remove Poisson Noise Corrupting CT Images

Kaies G1*, Noureddine E1, Chedly FM1 and Ines MM3

1El Manar University, LR-SITI, ENIT, Tunis, Tunisia

2Marzouk Moussa Ines, Diagnostic and Interventional Medical Imaging, University Hospital of Mongi Slim Marsa, Tunisia

3LR-12ES01 ( laboratory of general and digestive surgery), Faculty of medicine of Tunis, El Manar University, Tunis, Tunisia

*Corresponding Author:
Kaies G
El Manar University, LR-SITI
ENIT, Tunis, Tunisia
Tel: +21621590476
E-mail: [email protected]

Received Date: July 01, 2017 Accepted Date: August 14, 2017 Published Date:September 04, 2017

Citation: Kaies G, Noureddine E, Chedly FM, Ines MM (2017) Comparative Analysis between a Variational Method and Wavelet Method PURE-LET to Remove Poisson Noise Corrupting CT Images. J Comput Sci Syst Biol 10: 079- 084. doi:10.4172/jcsb.1000253

Copyright: © 2017 Kaies G, et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Visit for more related articles at Journal of Computer Science & Systems Biology


Poisson noise removal has been implemented by several methods and several approaches. This type of noise generally affects several types of images and in the whole case medical images especially those rebuilt after an X radiation. About it, as tests images are CT ones, further, we will present two methods derived from two different approaches to poisson noise removal. one that is variational, which is a restoring process that takes into account the actual appearance of this noise which is multiplicative, the other methods derived from wavelet approach which is among the growing methods for the Poisson noise removal mainly in medical dataset. This method is called PURE-LET, simply to say Poisson Unbiased Risk Estimation-Linear Expansion of thresholds. This latter is based on of an optimization of statistical tool Mean Square Error (MSE) using an unnormalized Haar decomposition wavelet transform (DWT) and exploiting the property of the orthogonality of these latter. An original comparison at the end of this paper allows us to assess the reliability and robustness of either the variational method or the PURE-LET one against the annoying factors in medical imaging such as artifact, spatial resolution, etc. Finally, we come to the conclusion that the Variational method excels at the Pure-Let wavelet method.


Poisson noise; CT images; Calculus of variations; Variational approach; Multiplicative noise; Image restoration; Wavelets; PURE-LET


Image restoration is a wide item study in the applied mathematics community as the same as the processing image community; it’s collaboration with the real case study and its parameters and the rigorous mathematical solution that fits. Without losing generalities; the test images chosen are the CT ones which are affected by the Poisson noise. This type of noise as the name suggests is considered a random variable that follows a Poisson distribution; which have the particularity of dissociating into distribution sum of Poisson random variable; isomorphic to the number of dissociation package; in addition; and against if one part goes to assemble the information each affected by a Poisson noise; such the image formation issue photons emission by an X-ray beaming; it still is a dataset that is affected by a noise resulting from sum of Poisson distributions outcome for each information. This feature is one of the important keys to the implementation of two algorithms. in the rest of work; we will briefly recall the formation of CT images; its need for a post treatment; here will be a restoring process; then we will implement the two methods; firstly; the variational approach; I’ll call him VAMPN-CT; just to say Variational Approach for Multiplicative Poisson noise in CT images which is implemented by myself; the second being PURE-LET which is implemented by their authors; you can consult [1] for more details. We will be pleased to describe here in brief the main tags of both methods. Further; a comparison between the robustness against noise and artifact of two methods and their reliability; which can be summarized in three points; PSNR (Peak in Signal to Noise Ratio); the presence of artifacts; and the algorithms’ runtime. Finally, we will draw conclusions after the useful comparisons.

CT Images’ Formation and their Quality

The principle of CT images formation resides in the emission of X-Rays on the object body screened under different rotation angles. Each angle allows us to generate some of measures which depends; on the absorption coefficients of the objects which constitute the body under examination and the thickness of the par of body under examination is constituted by several thickness volumes xxi having each one an absorption coefficient μi that means in reality:

μ. x=μ1. x2+μ2.x2+ ... +μn. xn (1)

To know the xxi and μi make up the CT (=Computed Tomography or Tomodensitometry) images, we have to make at least (n+1) measures; n being the number of unknowns xxi under each X-ray beaming angle. The μi are already known; it’s clear. We should make more than n measure; at least n+1 ones to remove the ambiguity which can be caused by the presence of several solutions generated by one system of n equations and unknowns values too. These measures aren’t directly exploitable; we have to rebuilding algorithms which are grouped into two families: the Analytics Reconstruction and the Iterative Reconstruction. For more details; the reader may consult with profit [2]. The quality of CT images can be evaluated by several parameters such as: spatial resolution; density resolution (the CT images are less contrasting) and the temporal resolution (present in the new scanner resolution where the patient is in movement). In addition; the artifacts (defects presented into image rebuilt; often appeared as blurs or stains which have varied origins) bore the quality of images. Several solutions are implemented by the radiologist to improve such or such parameter quoted before or to avoid artifacts according to the organ examined (soft or rigid) [2]. Some improvements which affect parameters confront some technical limits and hygienic (e.g., the dosage of medium contrast must be according to the patient body state and not to the contrast level which we want reach). Others degradations can be appeared caused by the presence of noises due several factors such the sampling of emitted signal and the arrondissement of calculus to form image matrix; but these noise are very negligible in front of the quantum noise [3]. So the CT images need as well a post-pretreatment. If we consider that the image pixels are mutually noisily independent; the noise would be multiplicative. Besides; according to the study explained into [3]; we can say that the quantum noise draws a POISSON process. Finally; the images treated are medical; hence; each almost details can influence the later result into medical work (e.g., diagnostic); thereafter a restoration of these images (slices) will be recommended. Therefore; the result of this study will be devoted to solve the mathematical model for the restoration of CT images in the presence of multiplicative noise that follows a POISSON process [4] (in this case; the successive photons’ emission follows each one a Poisson distribution at regular intervals which tend to zero).

Mathematical Study of Restoration Model Matching with VAMPN-CT

The observation of one photon emission ai outputted from examined body to contribute to form one-pixel image is stochastic and follows a Poisson distribution as the same as the noise which affects the formed data given by the formula above [3]:

equation (2)

ai: Number of photons outputted to form the ith pixel.

ki: Number of photons emitted by the source beaming to contribute to form the ith pixel.

λi: The Poisson distribution parameter which describe the formation of ith pixel.

We should add that the formation of each pixel is a result of photon’s emission issue several angles of X-Ray beaming. Poisson distribution is stable; that means if and are two independent random variables follows each one a Poisson distribution respectively Px(X) and Py(Y), so we have:

X → Px (X) and Y → Px (Y) then

X+Y → Px+y (X+Y) [5] (3)

The global nature of data or noise follows as in the same way Poisson distribution. In brief; it’s about a restoration of some CT images affected by a multiplicative noise which follows a Poisson distribution. We will look at the mathematical model of restoration.

Let f =u.v

where f the observed image; u>0 the image to recover and v the Poisson noise distribution. We consider that f and v are instances of some random variables F, U and V. In the following; we denote the density function matching with the random variable. The noise follows a Poisson distribution has λ as parameter with λ=n.p giving equation in our case study, indeed, n, given the number of measure drawn necessary to have one solution without ambiguity; that’s means; we need at least 1 +N.M=1+N.N with N and M (here M=N).

The image dimensions equations to calculate the unknowns xxi [2]; for the simple reason that a system formed by n equations and n unknowns can lead us to several n-tuple solutions or one n-tuple respond to real case. p, given the probability to determine one image pixel; all the images are equiprobable [6,7].


So that we have a generalized solution it’s necessary and sufficient that we take the equation to continuous case. As a result, in the following we denote gx the density function matching with the random variable X. In addition, view λ<2, we can’t in any case converge the Poisson distribution to Gaussian one (requires λ ≥ 5 or more). In brief, gv is the exponential distribution [8] with:

equation (4)

We suppose also that our image u is distributed as per a GIPPS distribution namely:

equation (5)

where z is a normalizing constant, μ is a constant and ξ is a positive application which describes texture of Image. In addition, let U, V two random variables independent of which their densities functions given respectively by gu and gv and let F=U.V then we have:

equation (6)

Our goal live in maximizing P (U|F) (trivial, since we have to recover u as from f), therefore, we need to know P (F|U) and gF|U ; thereafter, what precede, leads us to the classical Maximum a Posteriori (=MAP) estimator. From Bayes rule, we have

equation (7)

With U and V are two independent random variables. Maximizing P (U|F) get back to minimize the log-likelihood

−log(P(U | F) = −Log(P(F | U)) − log(P(U)) + P(F ) (8)

Thereafter, if we consider, our image as a pixel set S, moreover, the noise v is identically distributed on each pixel s∈ S and affects them independently in twos, with density gv then we have:


According to the formula (5) and (6), we can deduce:

equation (8)

Or Z=constant, therefore, our problem comes back to minimize the following functional:

equation (9)

In the continuous setting, we have to search: (9)


Where Ω the natural space in which computes our solution is named BV (=Bounded Variations); BV is the function space at bounded variations, it’s the necessary and sufficient to compute mathematically our solution matching to the previous functional. This space involves the bounded functions which present discontinuities, that’s why BV is ready to well support the mathematical description of the image in general. equation [8].

VAMPN-CT and its Euler-Lagrange Associated

The variational model which leaded, presents two difficulties, the first one comes from that the proposed model is no nonconvex on account of log function; therefore, we have to look for how to prove the existence and the uniqueness of solution. The second point is the non-differentiability of BV space which causes that the minimum of functional don’t verify an associated Euler-Lagrange equation [9]; therefore, we should look how to get round this difficulty. Without loss of generality, we take back the variational model:


and inspired from Rudin works [10] and [11], we pose ϕ(u)=|Δu|; we thus propose the following restoration model:

equation (10)

where S (Ω)={u∈ BV ((Ω)} with u>0 and λ being a regularization parameter [8]. We will compute equation as a distribution measure. As BV isn’t differentiable, we pose yet, as a distribution measure. As BV isn’t differentiable, we pose yet, equation and were then,

equation is a lebesgue space [12], we have so:

equation (11)

We pose equation and then, we use the Green’s formula applied to measurable functions [13], we obtain as the following:


where n is the outward normal to ∂Ω.

So equation we lead

equation (12)

Therefore, it’s clear that: We assume yet, the Neumann boundary condition [14] equation then, we can deduce that min J(u) comes back to resolve the following Euler-Lagrange equation:

equation (13)

n being the outward normal to the edge.

Existence and Uniqueness of the Solution


We remind: equation andequation where L1 is the dual space of L∞ orequation

With equation is a regular part.

equation is the singular part.

HN-1 is the Hausdorff measure with a dimension N-1.

equation is the classical gradient ofequation

equation s the jump part.

Cf is the Cantor part of Dsf

Sf is the jump set: equation

We have and equation it represents a negligible part. We assume it equal to zero, Cf={0}; with the Hausdorff dimension of the support of Cf is strictly greater than N-1 [8]. Moreover, f+ is the right discontinuity and f- is the left discontinuity; we have f+=f- except one set Sf having size N-1 as Hausdorff measure. Hence, we can deduce the following writing:


Δf: presents the real contour of image in question.

equation Presents a rectifiable curves which fit to image curve, we can lead to calculus of total variation DF, namely:

equation (14)

This proves the existence of a minimizer categorically.


Since we have proved the existence of a solution that goes along with the previous Euler-Lagrange equation, and since BV involves discontinued functions, we will assume that there are two solutions u1, u2 and we lead that these two solutions coincide i.e., u1, u2. We would resort to the following proposition:

Proposition 2: let equation

And equation We assume that f1<f2. We denote by u1 (resp. u2) a solution of (13) for f=f1 (resp. u2), then we have u1 ≥ u2. We consult [9] for further details. We prove that the set {u1>u} has zero as its Lebesgue measure for f1<f2 therefore f1=f=f2, we will necessarily have u1=u2.

Algorithm and Numerical Scheme


In this section and to numerically compute a solution to a problem equation (with λ is a regularization parameter). We use (13) and we embed it into a stationary equation which is a classical done in image analysis.

Discrete gradient and divergence operators

The most commonly used method to discretize the gradient operator and thereafter to calculate the divergence operator is the finite difference method; it comes from the structure of digital image which are formed by a set of pixels uniformly distributed. However, it exists several method to calculate the discreet gradient and thereafter the divergence operators, but, the method has been chosen, prove to be the most appropriate [8]. Let u the considered image and ui,j its value at pixel (i,j). The space step h is set to 1. The calculating of gradient operator by the finite difference method gives, applying the Neumann boundary condition, we have so to extend image on its borders (e.g., by symmetry), we have finally:

equation (15)

With N*N the image dimension

equation (16)

With equation The gradient norm chosen as the standard square:

equation With β is a small fixed parameter.

It would be fixed at 0.1 and this arrangement is just to avoid division by zero and this would be in any case where the division by zero appears. We also introduce a discrete version of the divergence operator and this is by analogy to the continuous case by div=-∇*, where ∇* the adjoint of is ∇ thus:

equation (17)

equation (18)

We note that we have led to previous result by using integration by parts and calculate the derivative distribution measure.Orequation with ∇ is the Laplace operator.That, we would suggest calculatingequation, then we apply the discreet divergence operator posingequation Was thenequation always with the remark to add a small fixed parameter (β=0.1) to avoid division by zero in the order of the matching indices i, j.

Numerical scheme and mainly parameters

No universal method to calculate λ (parameter of regularization), even if there some try in the literature to make the study of this parameter methodical and subjected to theoretical formulas. The experimental study remains the most efficacious to release λ optimum due the complexity of function u (=which describe the image) in BV space and rarely we can find an inverse to this function [8]. We assume a fixed regularization parameter lambda (is chosen to be equal to 0.1), which always gives high quality images and in the majority of cases this value is the optimal.

It should be noted that the correlation coefficient (which is related to the calculation of the matrix product within the algorithm, view, that the range of intensities is of the order of [a<=-2000b>=2000] (in CT images) is of the order of 10-7 (is the best) or 10-6, (this is trivial since the product of two matrices showing intensities of the order of 107, lower or upper, gives degraded images.

The Pure-Let Approach


The main idea of Pure-Let based on an optimization of a statistical tool MSE which measures the risk between the noiseless (restored) image and the noisy image. Assumption and indeed that is the reality, the biomedical images are affected by the Poisson noise, as at least to my knowledge for images acquired by the X-ray emission. This estimator is optimized independently for each sub-band by adopting an unnormalized Haar DWT (Haar discrete wavelet transform), exploiting mainly the property of the orthogonality decomposition of these latter and the stability of Poisson distribution (the sum of independent Poisson random variable is also a Poisson random variable whose intensity is the sum of original intensities). We can easily switch to 2D, exploiting the fact of separability of Haar wavelets. Finally, we lead to a fast algorithm which requires reduced memory consumption and low complexity and according to their authors gives a reliable results manifested by a suitable and acceptable denoising image.

Pure-Let and mainly steps

In this paragraph, we are content to remind roughly Pure- Let major lines and stages of implementation, for the curious and interested readers will be referred to reference [1] for further details and explanation and elsewhere this article is performed by the authors of PURE-LET. We will quote the steps and some mathematical properties that will be the keys of this approach.

First step: We apply the unnormalized Haar discrete wavelet transform (DWT) to our dataset to obtain the scaling coefficients sj at the measurement m for the order of scaling j=1 J (J will be chosen upper than 2 in the implementation of algorithm) and the associated wavelet coefficients dj, we denote that sj, dj ∈ RNj where Nj=N/2j and N is the signal dimension; we assume that the signal is divisible by 2j setting s0=m, these coefficients are obtained from the following sums and differences as the following:

equation (19)

The original sequence m=s0 is simply recovered by computing:

equation (20)

Similarly, we denote by σj and δj scaling and wavelet coefficients of the original signal μ at a given scale j. Note by linearity of the wavelet transform, we have equation then we have thanks to the properties of unnormalized Haar DWT the following results:

• It’s an orthogonal transform, so we can split the MSE into subband-specific error terms:

equation (21)

This implies that we can minimize the MSE for each subband independently, while ensuring a global signal-domain MSE minimization.

• As it’s explained before, that at a given scale j, the scaling coefficients of an input vector of independent Poisson variables are also independent Poisson random variables.

Second step: In this section, we will compute the estimate equation andequation depends on the associated wavelet coefficients dj and the scaling coefficients dj too. This leads to the following functional:

equation In practice and as usual the low pass residual is not processed i.e., equation hence theequation and nextequation which will be minimized independently for each wavelet subband. The authors of PURE-LET lead to prove equation that is an estimate of the noise-free wavelet coefficients δ =δ j and then the random variable εj given by the following formula:

equation (22)

Is an unbiased estimate of the MSE for the subband under consideration, equation andequation

The vector equation are considered as the canonical basis of equation Finally, we lead at unbiased μ details (not affected by the Poisson noise).

It should be mentioned that we have to extend the later process into 2D which is an easy process, thanks the separable Haar DWT.

Third step: In this part, we will introduce more elaborate functions, we propose to consider a wavelet estimator that is formulated as a linear expansion of thresholds (LET), i.e.,

equation (23)

where equation are generic estimators that we can be chosen arbitrarily (in this study, we will specify more elaborate equation y taking into account the redundant information across scales and it will be shown that the choices highlighting improving denoising quality. The MSE is quadratic (in our case) with respect to the parameters equation therefore, and thanks to the linear parameterization in this case, its minimization boils down to the resolution of a linear system of equations with small dimension k as the following: equation where:

For 1 ≤ k,l ≤ K (24)


With equation and equation are similar as above.

Besides, we propose a linearly parameterized thresholding function with K=2 parameters (a1 and a2), whose n-th component is defined by:

equation (25)

In this formula the linear parameters a1 and a2 define a compromise between two regimes: either the wavelet coefficient dn is kept as is (signal preservation) or it’s shrunk towards zero (noise suppression). The exponential function has the advantage of being smooth, which reduces the variance of estimator; that is, to be closer to the original signal (2D of course). Several improvements can be introduced to denoising process such as adding a proportional term of the interscale predictor into the thresholding function quoted above; Note that to construct an interscale predictor of the wavelet coefficient dn. we simply compute the difference between the two scaling coefficients that surround equation Further improvements can be obtained by grouping wavelets coefficients of similar magnitudes [15] to increase the robustness of similar noise. For more details, we refer the reader once again to reference [1]. What remains is to experimentally adjust or from equations described well the key parameters of the algorithm such as PURE-LET order of scale or LET as well as the algorithm denoising parameters and so on. Finally, we should apply the matching inverse transform by the wavelet reconstruction process to get ultimately an optimal estimate of the original image without the most of the Poisson noise, taking into account the particularities of CT images as the pixel record size and the range of scattering intensities, etc.


We will implement in rotas, the two algorithms, the first being the variational approach as I called VAMPN-CT (just to say Variational Approach for the Multiplicative Poisson Noise in CT images) with Matlab R2011a, as the characteristics of the machine being 2.00 GHZ the frequency of the microprocessor, 2 CORE DUO and 2GB RAM. Likewise, for the second algorithm which is The PURE-LET approach which based on wavelet harmonic analysis. We have tested several Dicom slices (Figures 1-3). We have leaded to results which could be resumed by the Table 1. The PSNR, the runtime and the presence of artifacts are the comparison parameters.


Figure 1: Observed image slice of hepatic nodule.


Figure 2: Restored image by VAMPN-CT method.


Figure 3: Restored image by PURE-LET method.

PSNR (db) 22,1O22 52,848
Run time (s) 77,05 17,35
Artifacts None Star-type artifact

Table 1: Comparative analysis between VAMPN-CT and PURE-LET.


First of all, it’s clear that the variational approach VAMPN-CT excels at the wavelets approach PURE-LET. We must focus on the presence of artifacts for the second approach (PURE-LET), which is very troublesome and unacceptable in most time by the Radiologists’ community. The VAMPN-CT leads to highest PSNR in the CT images case study, therefore, a compelling contrast resolution and even spatial resolution. The runtime has no meaning if the image quality is bad and especially the difference between the two don’t suggests to us to choose the faster; unlike and according the reference [1]. We should add that the artifacts are insidious in many experiences related with CT images and elsewhere in our case. Without forget that the PURE-LET method not stable; giving different PSNR for a given image whenever we repeat the implementation of its algorithm. We would compare our method to the PLATELETS [16] one and we would be ambitious that our approach VAMPN-CT excels on PLATELETS in image quality and even in runtime according to [1], but that it must be demonstrated experimentally.


Select your language of interest to view the total content in your interested language
Post your comment

Share This Article

Relevant Topics

Recommended Conferences

Article Usage

  • Total views: 223
  • [From(publication date):
    September-2017 - Dec 15, 2017]
  • Breakdown by view type
  • HTML page views : 201
  • PDF downloads : 22

Post your comment

captcha   Reload  Can't read the image? click here to refresh

Peer Reviewed Journals
Make the best use of Scientific Research and information from our 700 + peer reviewed, Open Access Journals
International Conferences 2017-18
Meet Inspiring Speakers and Experts at our 3000+ Global Annual Meetings

Contact Us

Agri & Aquaculture Journals

Dr. Krish

[email protected]

1-702-714-7001Extn: 9040

Biochemistry Journals

Datta A

[email protected]

1-702-714-7001Extn: 9037

Business & Management Journals


[email protected]

1-702-714-7001Extn: 9042

Chemistry Journals

Gabriel Shaw

[email protected]

1-702-714-7001Extn: 9040

Clinical Journals

Datta A

[email protected]

1-702-714-7001Extn: 9037

Engineering Journals

James Franklin

[email protected]

1-702-714-7001Extn: 9042

Food & Nutrition Journals

Katie Wilson

[email protected]

1-702-714-7001Extn: 9042

General Science

Andrea Jason

[email protected]

1-702-714-7001Extn: 9043

Genetics & Molecular Biology Journals

Anna Melissa

[email protected]

1-702-714-7001Extn: 9006

Immunology & Microbiology Journals

David Gorantl

[email protected]

1-702-714-7001Extn: 9014

Materials Science Journals

Rachle Green

[email protected]

1-702-714-7001Extn: 9039

Nursing & Health Care Journals

Stephanie Skinner

[email protected]

1-702-714-7001Extn: 9039

Medical Journals

Nimmi Anna

[email protected]

1-702-714-7001Extn: 9038

Neuroscience & Psychology Journals

Nathan T

[email protected]

1-702-714-7001Extn: 9041

Pharmaceutical Sciences Journals

Ann Jose

[email protected]

1-702-714-7001Extn: 9007

Social & Political Science Journals

Steve Harry

[email protected]

1-702-714-7001Extn: 9042

© 2008- 2017 OMICS International - Open Access Publisher. Best viewed in Mozilla Firefox | Google Chrome | Above IE 7.0 version