Comparative analysis between a variational method and wavelet method Pure-Let to remove Poisson noise corrupting CT images Текст научной статьи по специальности «Компьютерные и информационные науки»

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

Ghazouani Kais, Ellouze Noureddine, Marzouk Moussa Ines

Abstract—The 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.

Keywords— Poisson noise, Computed tomography (CT) images, calculus of variations, variational approach, multiplicative noise, image restoration, wavelets, PURE-LET.

I. Introduction

Image restoration is a wide item study in the applied mathematics community as the same as the processing image community, it's a collaboration with the real case study and its parameters and the rigorous mathematical solution that fits.

Ghazouani Kaies, LR-SITI (Laboratory Signal, Image and Technologies of information, National School of Engineers of Tunis, University of EL Manar, Tunisia and phone: +21621590476, email: ghks009@gmail.com. Ellouze Noureddine, LR-SITI, National School of Engineers of Tunis, University of EL Manar, Tunisia, email: N.ellouze@yahoo.fr Marzouk Moussa Ines, LR-12ES01 (Laboratory f gastrointestinal diseases), Medical School of Tunis, University of EL Manar Tunisia, email: Marzouk.ines@gmail.com

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 Signal Noise Ratio), the presence of artifacts, and the algorithms' runtime. Finally, we will draw conclusions after the useful comparisons.

II. 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 x, having each one an absorption coefficient ^ that means in reality:

fl-X = f4.X2 +/fa.x3 + ■■' (1)

To know the and to make up the CT (=Computed Tomography or Tomodensitometry) images, we have to make at least (n+1) measures; n being the number of unknowns under each X-ray beaming angle. The are already known, it's clear. We should making 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).

III. 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]:

nfef) = Ï^a^kil^-Ilt*«) =

-a,-

(2)

X - iiOT) and Y ¿>00 then X + Y ~Px+y(X + V) [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, I/and V. In the following, we denote gx the density function matching with the random variable. The noise follows a Poisson distribution lias A as parameter with A =n.p giving

A =■

BJB

a* 1 in our case study, indeed,?! 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 = the image dimensions equations to calculate the unknowns [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 p = ttt = ~ and

images are equiprobable p = 77777 = 777; ■:=.-.- .'.'..'.' .-. e :■'. 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 A < 2, we can't in any case converge the Poisson distribution to Gaussian one (requires A > 5 or more). In brief, gv is the exponential distribution [8] with:

0 if v < 0

(4)

exp(i0 if v> 0

:- : Number of photons outputted to form the i pixel.

L.' : Number of photons emitted by the source beaming to

contribute to form the ith pixel.

: 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 X and Y are two independent random variables follows each one a Poisson distribution respectively Py (aJ and (1Q, so we have:

We suppose also that our image u is distributed as per a GIPPS distribution namely: 3u(u) =-zBxp{-u<p{u)-) [6] (5)

Where 2 is a normalizing constant, u is a constant and 'p is a positive application which describe Texture of Image. In addition, let U, V two random variables independent of which their densities functions given respectively by .gL.L. ge and let F = U.V then we have:

ftO-; = *i-0T»> m (6)

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

we

have: FÖ/li7) =

p (F|t/).PC[J]

(7).

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

—logffCl/lF) =-hog (P (F|i/3)- Log (f (U)~) + F(F) (8) Ori7 is known that gives: P(F') = 1 => log (F(F| j = 0 . Thereafter, if we consider, our image as a pixel set moreover, the noise v is identically distributed on each pixel .ve S and affects them independently in twos, with density gv, then we have:

PtFlfcO = n^s^^Wli/fc)) , hence:

-logV^ffWIffW))- =

-^[tofftFWI"«)] togfa, ))

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

log(PtU\F')) =

SSEJ ^(uk)) + + ¡t. 0(u(s» + logtä

(8)

Or inconstant, therefore, our problem comes back to minimize the following functional:

, in the continuous

setting, we have to search:

mi*[Jfl G°g(«U) + dx + PJn *(«(*)**] (9)

Where J2 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. BV d1 [8].

IV. VAMPN-CT AND ITS EULER-LAGRANGE ASSOCIATED

J(ü = lim/.

VCu + tvJ - Vu I

+ÀLim_T

h (y + tvj- h(u)

dx

t-io ' " t

h(u + tvj - h(u)

tf p

äx

Viil

= /n Vv. dx + ljah'(ii). v. dx

'a Vli

^^Vr.dx + Jj^.^-ix. (11)

We pose A = j' — Vv.dx and B = ^ .v.di and

Li-/

< a pu ' - — J ci

then, we use the Green's formula applied to measurable functions [13], we obtain as the following:

{ Vu

'il Vu

Where n is the outward normal to 3ÎL So v v e l^fal), wq lead :

Vu Vu\

Therefore, it's clear that:

We assume yet, the Neumann boundary condition [14] Fu

: —. « = 0 (OndS2 ] Then, we can deduce that min )(;ii)

comes back to resolve the following Euler-Lagrange equation:

-div(^) + ),h'(M) =0

ivul

MTul .n = 0

(13)

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:

Min J (Jog [ uixjj ¿j / <p(uG0 )dx and inspired

from Rudin works [10] and [11], we pose ^ (uJ = Idul; we thus propose the following restoration model:

№ = i^fljl^l^ + Sa + £] dx (10)

Where 5 03]= j> e BV CO)} with u > 0 and X being a regularization parameter [8]. We will compute /60 as a distribution measure. As BV isn't differentiable, we pose yet, = login) and were then:

[12], we have:

n being the outward normal to the edge.

V. EXISTENCE AND UNIQUENESS OF SOLUTION

A. Existence

We remind: sftf) = f / e BV(J2) oaf > 0} and BViii) = if G LLtO) at |d/| dx < +<x} where L1 is the dual space ofi".

or Df = Vf. dx + (/ + - /-Jp^JfjSj1 +-Cf

With Vf. dx is a regular part.

is the singular part.

is the Hausdorff measure with a dimension N — 1. Sf if

Vf is the classical gradient of/ (—, —].

- ■:-■-'-.'■;:■ is the jump part. Cf is the Cantor part of Dsf. Sf is the jump set: {x e Q;f~ CO < f+(x)} We liave C;- 1 dx and it represents a negligible part. We assume it equal to zero, C,-= {0}; with the Hausdorff dimension of the support of C,- is strictly greater than N-l [8] .Moreover, /+is the right discontinuity and / "is the left

discontinuity; we have /+=/ except one set Sr- having size N-1 as Hausdorff measure. Hence, we can deduce the following writing:

C r - f-y^K »j1 = tf+- f-yw-1

¡Presents the real contour of image in question. - .-" ■:-■>.'■;;-: Presents a rectifiable curves which fit to image curve, we can lead to calculus of total variation DF,

namely:

ID/Itfl] =jn\Df\

+Ja_sf\Cf\ (14)

- uui if i < N

With N*N the image dimension

(16)

This proves the existence of a minimizer categorically @

B. Uniqueness

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 ui,u2 and we lead that these two solutions coincide i.e. ut = uz . We would resort to the following proposition:

Proposition: let A and f2 E ¿"with InfVi/L > 0 And [iifn/1; > 0. We assume that Ji < f2- We denote by ui (resp.H23 a solution of (13) for/ = fL(resp.f =/2}, then we have > ii2 .

We consult [9] for further details. We prove that the set > lias zero as its Lebesgue measure for f < f2 therefore A =/ =/a , we will necessarily liave =

VI. ALGORITHM AND NUMERICAL SCHEME

A. Introduction

In this section and to numerically compute a solution to a problem -dir (j^q 1 + A. (^r) = 0 (with A is a

regularization parameter). We use (13) and we embed it into a stationary equation which is a classical done in image analysis.

B. Discrete Gradient and Divergence

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 ll the considered image and 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:

Witli:S™cKu) = fTE/}*,.. + (FUyij. The gradient norm chosen as the standard square: IFt/| =^(Vxy + (Fj)2 +p2 ; With jl 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 = —V", where F" the adjoint of is F Thus:

{divfr) i : = (div G0)£J + {div with:

{div GO)!1,

i-j

tj-Pii-Uifl<i<H *£j if i = i -PU-ifi=N

i^ii/1 <j <«

I-kj-i ifJ = N

(17)

(18)

We note that we have led to previous result by using integration by parts and calculate the derivative distribution measure. Or dip ('j^ 1 /(Al/] with A is the Laplace

FU

operator. That, we would suggest calculating(j^7p , then we apply the discreet divergence operator posing jj

f—) = c-

VlFul' N

(—) Vul*

FL? FE?

; always with the

Was then: ,. .. .. . ..

remark to add a small fixed parameter (/? = 0,1) to avoid division by zero in the order of the matching indices i, j.

C. Numerical scheme and mainly parameters No universal method to calculate A (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 A optimum due the complexity of function u (=wliich 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 <= - 2000 b> = 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.

VII. THE PURE-LET APPROACH

A. Introduction 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.

B. 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. 1) First step:

We apply the unnormalized Haar discrete wavelet transform (DWT) to our dataset to obtain the scaling coefficients S* at the measurement m for the order of scaling / = 1 ..../(/ will be chosen upper than 2 in the implementation of algorithm) and the associated wavelet coefficients riJ, we denote that s*.dj E M^1 where Njl' and N is the signal dimension; we assume that the signal is divisible by 2s setting s° = m, these coefficients are obtained from the following sums and differences as the following:

. For; = 1 —J (19)

The original sequence m=s° is simply recovered by computing:

■ " ^ Forj 1...J (20)

Similarly, we denote by a} and fi'J the scaling and wavelet coefficients of the original signal p at a given scale j. Note by linearity of the wavelet transform we have <p( d^) = cj-' and ijj(ji^) = [j^. then we have thanks to the properties of unnormalized Haar DWT the following results: a) It's an orthogonal transform, so we can split the MSE

into subband-specific error terms: (21)

This implies that we can minimize the MSE for each subband independently, while ensuring a global signaldomain 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.

2) Second step:

In this section we will compute the estimate SJ and ; SJ depends on the associated wavelet coefficients and the scaling coefficients s1 too. This leads to the following functional:

fiJ = &J (dJ,sy). In practice and as usual the lowpass residual is not processed i.e. a* .s-1 , hence the

MSE=Ij=1^ Si-Si : and next

which will be minimized

independently for each wavelet subband. The authors of PURE-LET lead to prove that 8*(d*. sJ~)is an estimate of the noise-free wavelet coefficients S = and then the random variable s. given by the following formula:

dT(9~fjl.s) + 9+(.d,s))-sT(9~(.d,s) - 9+(.d,s)~))

Is an unbiased estimate of the MSE for the subband under consideration, i.e., ¿¡(V) = %(MSEJ') with: 8£(.d,s) = fi^Cd + Bn,s - en] and

The vectors(en]ll=UJVi are considered as the canonical basis

of E"v'. Finally, we lead at unbiased |i 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.

3) 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. BLE7 (d, s.a) = Sf=1 fj, jy) Where the^ 's are generic estimators that we can be chosen arbitrarily (in this study, we will specify more elaborate && by 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 ael', 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: a — M_1c where :

For 1 < k, ï < K,

cfc [ dT Cd, s) + (d, e0) + s7 (4 s) - 6 + Cd. «))]/2

PURE-LET

VARMPN-CT PURE-LET

PSNR (db) Image 1 34,2111 5,3127

Image 2 33,115 5,2997

Runtime(s) Image 1 71,49 16,76

Image 2 71,33 15,98

Artifacts Image 1 none star-type artifact

Image 2 none Star-type artifact

Mw = e1ICdjS]7e][:djS] (23)

With S+(.d, sj and &~ (d, s) are similar as above. Besides, we propose a linearly parameterized thresholding function with K=2 parameters (ax and a2), whose n-th component is defined by:

(24)

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 nfl \vc simply compute the difference between the two scaling coefficients that surround fpj of; = 5,rl_1 — sf!+1. 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 [15]. 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.

VIII. RESULTS

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 (Digital Imaging Communication in Medicine) slices. We have leaded to results which could be resumed by the table below (comparative analysis between VAMPN-CT and PURE-LET) and the PSNR, the runtime and the presence of artifacts are the comparison parameters.

The labeled images are organized as following in order: the original image, then, the restored image by PURE-LET and finally the restored image by VARMPN-CT. We met to present here two images among several we tested and it's always the same types of results:

Further, the labeled images below; we have chosen one original image restored by the approaches already described. Finally, we should note that the images are 512*512 of size.

Original Image

Figure 1. Observed Image1

oLiiR +■ a.2 [1 — Bxp

Table 1. Comparative analysis between VAMPN-CT and

Restored Image PURE-LET , PSNR = 5.3127

Restored Image PURE-LET , PSNR = 5.2937

Fig. 2. Restored Image1 PURE-LET

Fig. 5. Restored Image2 PURE-LET

Restored Image VAMPN-CT. PSNR=34.2111 ISNR= 26.303

Restored Image VAMPN-CT. PSNR=33.115 ISNR= 23 7364

Fig. 3. Restored Image1 VARMPN-CT

Fig. 6. Restored Image2 VARMPN-CT

Original Image

We should note that the images which already tested, are well pretreated so-called by the radiologist by another method integrated with The MYRIAN2 software, except the pre-treatment methods who performs upstream. If we took an image which pre-treated only by the radiologist, we reach the results as following:

Fig. 4. Observed Image2

2 Myrian Software: A French software, specializing in the processing of CT-images and detection of tumor lesions, in particular, liver tumors.

Observed Image

Fig. 7. Observed image3 in arterial phase

Restored Image VARMPN-CT, PSNR=1124723 Lambda=0.1

Fig. 8. Restored image3 by VARMPN-CT, PSNR=112, 472

Restored Image PURE-LET , PSNR = 3 3743

Fig. 9. Restored Image3 by PURE-LET, PSNR=6, 374 We can see clearly the poor quality of the image restored by

the PURE-LET method (noise and volume partial). So, if we would yield the images medically significant, we apply an adequate fenestration according to the parenchyma studied (here the hepatic parenchyma for subsequent treatment of any liver nodules), we will then:

Restored Image VARMPN-CT, PSNR=112 4723 Lambda=0.1

Fig. 10. Restored Image3 by VARMPN-CT [Center=38, Width=425]

Restored Image PURE-LET , PSNR = B.3743

Fig. 11. Restored Image3 by PURE-LET [Center=42, Width=402]

We should note that the images are coded on 16 signed bits per pixel. Further the scratches which are visible just to the edges of the image, are the traces of the scanner table which can be eliminated easily with a simple approved mask applied to the image.

IX. CONCLUSIONS

Without losing generality, if we compare the two methods, we can say that the variational approach VARMPN-CT excels on the approach of wavelet PURE-LET, as regards the calculation of PSNR and eventually other parameters such as ISNR, MSE which are proportional to first, even, the

contrast resolution factor that has importance especially in medical imaging; PURELET excels in run time but with similar image quality, the runtime is meaningless; I focus mainly on the presence of artifacts called in this case startype artifacts is an annoying factor for image analysis, diagnostics and sometimes the image is out of Service (illegible);

The PURE-LET runtime is a little high, but still we seek an optimum solution will be in spite of run time; for example, we can be satisfied to 5 iterations in the VARMPN-CT algorithm and gives very good results always higher than that of PURE-LET, further the run time decreases about 14 seconds leaving the VARMPN-CT time is about 55 seconds. For a variational method, the indicated runtimes is very well if we are certain that the ratio between it and the wavelets one can be circumvented by the current progress in computers performance and the radiologist carry favor this latter, especially that this method doesn't let escape artifacts (some of them are insidious and not appears for the first shot) unlike others, particularly PURE-LET. Afterwards, it's clear that the variational approach VAMPN-CT excels on the wavelets approach PURE-LET. We must focus on the presence of artifacts for the second approach (PURE-LET), which is very troublesome as explained before and unacceptable in most time by the Radiologists' community.

The VAMPN-CT leads to highest PSNR, at least, in this case study and the CT images, 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], the PURE-LET approach is very faster that the PLATELETS [16] (a wavelet approach on which we implement: a multiscale approach for recovering edges and surfaces in photon-limited medical imaging) but the PSNR of the two approaches are very close, sometimes the chicken PSNR of PURE-LET is the most optimal and sometimes that of platelets prevails but with very small differences in the range from 0.1 or a little less. What makes us think our approach compared to that of platelets (of course the both approaches are derived from two different families but we will compare the image quality restored especially in the presence of Poisson noise for many types of medical images). Once again, according to the study done in [1], we would be ambitious that our approach VARMPN-CT excels on PLATELETS in image quality and even in runtime but that it must be demonstrated experimentally and elsewhere tried to have the platelets code from their authors, they didn't accept the challenge; In addition, these results are approved and confirmed by the radiologist Mrs. Ines Moussa Marzouk, aggregate radiologist, in the radiology center of Mongi Slim hospital situated in Marsa city, Tunisia and among their specialty, the hepatic imaging. This certificate said radiologist is as follows:

Dr Marzouk Moussa Ines, Associate Professor in Diagnostic and , Interventional medical imaging at the university hospital of Mongi ri Marsa Tunisia teaching in the medical school of Tunis, university of Tunis Manar. Performed publications in early drug evaluation research with the team of S1T£P in the Gustave Roussy institute, Villejuif, Franco in 2008. Member of a research laboratoiy in the field of gastrointestinal diseases in the medical school of Tunis. Working on elastography in breast imaging since 2012.

I undersigned Ms Ines Moussa Marzouk., CONFIRM having tested the Restore method of CT images implemented by Mr (rhazouani Kates and confirm that the method is efficient for noise elimination affecting CT images obtained from a multiplies machine (Aquilion 64 TOSHIBA 2012), also i confirm this method tested by Mr Ghazouani gave better results than the means we have implemented on the same machine for the quality of the picture and so for the medical interpretation. We performed a short survey with radiology resident* and seniors working in The University hospital of Mongi Slim Marsa to appreciate these results,

¡ICH*

To finish, he just said that this last variational method is effective to restore an image affected by a multiplicative Poisson noise. To eliminate artifacts that have different origins, it will be convenient to use other image processing processes unless the artifact derives Poisson noise. Finally, we can apply this method to any other parenchyma issue CT-imaging.

References

[1] F. Luisier, C. Vonesch, T. Blu, M. Unser. ''Fast Interscale Wavelet Denoising of Poisson-corrupted Images''. Signal Processing, vol. 90, no.12, pp.415-427, February 2010.

[2] Régent, A.Lisbona, F. Masson, A. Noel, "Scanner and X-rays", ''Scanner et Rayons X'', Elsevier Masson, 30 october 2013, 184p. Medical Imaging, "Imagerie Médicale''.

[3] P.L. Combettes et H. J.Trussell. ''Models and Algorithms for restoring digital X-ray images'', ''Modèles Et Algorithmes en vue de la Restauration Numérique D'images Rayons X''. MARI (Machines, Réseaux Intelligents), Image Electronique. Mai 1987, COGNITIVA 87.

[4] D. Foata, A. Fuchs"Poisson Process,"Markov chain and martingales , '' Processus de Poisson, chaînes de Markov et martingales'' , 14 october 2004, 236p,Dunod, '' Processus Stochastiques''

[5] G. Grimmett, D. Stirzaker: Probability and Random Processes. Oxford University Press, 2009.

[6] S. Battiato, " Entropy and Gibbs Distribution in Image Processing: an Historical Perspective'' - Image Processing & Communications, an International Journal ", Vol.7, No. 3- 4, pp. 35-50, 2001.

[7] G. Grimmett et D. Welsh. Probability, '' An introduction. Manual of Remote Sensing''. Oxford 23 Science Publications, 1986.

[8] G. Aubert and P. Kornprobst. "Mathematical Problems in Image Processing, Partial Differential Equations and the Calculus of Variations". Second Edition. United States. Mathematics Subject Classification 2000, 373 p. Applied Mathematical Sciences Vol. 147

[9] J.F.Aujol, G. Aubert, ''A Variational Approach to Remove Multiplicative Noise,'' SIAM Journal on Applied Mathematics, 68(4), 925-946 January 2008.

[10] L. Rudin, P-L. Lions, and S. Osher. Multiplicative denoising and deblurring: Theory and algorithms. In S. Osher and N. Paragios, editors, Geometric Level Sets in Imaging, Vision, and Graphics, pages 103-119. Springer, 2003.

[11] L. Rudin, S. Osher, and E. Fatemi. Nonlinear total variation based noise removal algorithms. Physica D, 60:259-268, 1992.

[12] C-U.David, A.Fiorenza. 'Variable Lebesgue Spaces, Foundation and Harmonic Analysis'. Springer Basel, 13 February 2013, 312 p. Applied and Numeric Analysis.

[13] I. Stakgold, M. Holst, "Green's functions and Boundary value

Problems", Third Edition, United States, Pure and Applied Mathematics,875p. A WILEY SERIES OF TEXTS, MONOGRAPHS, AND TRACTS.

[14] G.Cristobal, P.Schelkens, H.Thienpont. "Optical and Digital Image Processing: Fundamentals and Applications". Electrical & Electronics Engineering. 29 April 2011, 879p. Imaging Systems and Technology.

[15] F. Luisier, T. Blu, M. Unser, A new SURE approach to image denoising: interscale orthonormal wavelet thresholding, IEEE Transactions on Image Processing 16 (3) (2007) 593-606.

[16] R.M. Willett, R.D. Nowak, Platelets: a multiscale approach for recovering edges and surfaces in photon-limited medical imaging, IEEE Transactions on Medical Imaging 22 (3) (2003) 332-350.

Ghazouani. Kaies, PhD student at who won the

Masters in 2008, in the process of achieving an expert diagnostic aid system in medical imaging works as a teacher of Computer Science seconded to higher education since 2001, he joined the LR-SITI research team in 2006 and performs some publications in image processing that are being especially restoring images and the related mathematical problems.

Ellouze. Noureddine received a Ph.D. degree in 1977 from l'Institut National Polytechnique at Paul Sabatier University (Toulouse-France), and Electronic Engineer Diploma from ENSEEIHT in 1968 at the same

University. In 1978, Dr. Ellouze joined the Department of Electrical Engineering at the National School of Engineer of Tunis (ENIT-Tunisia), as assistant professor in statistic, electronic, signal processing and computer architecture. In 1990, he became Professor in signal processing; digital signal processing and stochastic process. He has also served as director of electrical department at ENIT from 1978 to 1983. General Manager and President of the Research Institute on Informatics and Telecommunication IRSIT from 1987-1990, and President of the Institute in 1990-1994. He is now Director of Signal Processing Research Laboratory LSTS at ENIT, and is in charge of Control and Signal Processing Master degree at ENIT. Pr Ellouze is IEEE fellow since 1987, he directed multiple Masters and Thesis and published over 200 scientific papers both in journals and proceedings. He is chief editor of the scientific journal ''Annales Maghrébines de l'Ingénieur''. His research interest includes neural networks and fuzzy classification, pattern recognition, signal processing and image processing applied in biomedical, multimedia, and man machine communication.

Dr Marzouk Moussa.Ines, Associate Professor in Diagnostic and Interventional Medical imaging at the university hospital of Mongi Slim Marsa Tunisia teaching in the medical school of Tunis, university of Tunis Manar. Performed publications in early drug evaluation research with the team of SITEP in the Gustave Roussy institute, Villejuif, France in 2008. Member of a research laboratory in the field of gastrointestinal diseases in the medical school of Tunis. Working on elastography in breast imaging since 2012.