Jules Sadefo Kamdem
email: sadefo@lameta.univ-montp1.fr

Coecient of variation and Power Pen's parade computation

Keywords: Gini index, Income inequality, Ranks, Harmonic Number, Pen's Parade. JEL Classication: D63, D31, C15

Under the the assumption that income y is a power function of its rank among n individuals, we approximate the coecient of variation and gini index as functions of the power degree of the Pen's parade. Reciprocally, for a given coecient of variation or gini index, we propose the analytic expression of the degree of the power Pen's parade; we can then compute the Pen's parade.

Introduction

Interest in the link between income and its rank is known in the income distribution literature as Pen's parade following [START_REF] Pen | Income Distribution[END_REF][START_REF] Pen | A parade of dwarfs (and a few giants)[END_REF]. The precede has motivated research on the relationship between Pens parade and the Gini index that is a very important inequality measure. However, this research has so far focused on the linear Pens parade for which income increases by a constant amount as its rank increases by one unit (see [START_REF] Milanovic | A simple way to calculate the Gini coecient, and some implications[END_REF].

Furthermore, a linear pen's parade does not closely t many real world income distributions Pens parade of which is convex in the absence of negative incomes (i.e., income increases by a greater amount as its rank increases by each additional unit). [START_REF] Mussard | Quadratic Pen's parade and the computation of Gini Index[END_REF] has recently introduced the computation of Gini index with convex quadratic Pens parade (or second degree polynomial) parade for which income is a quadratic function of its rank (see also [START_REF] Ogwang | The gini index for a quadratic prn's parade[END_REF]).

This paper extends the computation of gini index by using a more general and empirically more realistic case for which Pens parade is a power function. Hence, the Gini indices for a linear Pen's parade and for quadratic Pen's parade becomes special cases of our thus introduced power Pens parade under some constraints on parameters inducing convexity. Under the assumption that income y is a power function of its rank among n individuals, we approximate the coecient of variation and gini index as functions of the power degree of the pen's parade. Reciprocally, for a given coecient of variation or gini index, we propose the analytic expression of the degree of the power pen's parade. It follows that, for a given coecient of variation or gini index we can compute the pen's parade. We use the [START_REF] Milanovic | A simple way to calculate the Gini coecient, and some implications[END_REF] data to illustrate the interest of our results.

The rest of the paper is organized as follows: In Section 2, the specication of an a power Pens parade is provided. In Section 3., the problem of tting a power Pens parade to real world datasets of [START_REF] Milanovic | A simple way to calculate the Gini coecient, and some implications[END_REF] is discussed.

The concluding remarks are made in Section 4.

High degree Power Pen's Parade

Suppose that positive incomes, expressed as a vector y, depend on individuals' ranks r y in any given income distribution of size n. Suppose that incomes are ranked by ascending order and let r y = 1 for the poorest individual and r y = n for the richest one. Hence, following [START_REF] Lerman | A note on the calculation and interpretation of the Gini index[END_REF], the Gini index may be rewritten as follows:

G = 2 cov(y, r y ) n ȳ . (1)
Here, cov(y, r y ) represents the covariance between incomes and ranks and ȳ the mean income. It is straightforward to rewrite (13) as:

G = 2 σ y σ ry ρ(y, r y ) nȳ , (2)
where ρ(y, r y ) is Pearson's correlation coecient between incomes y and individuals' ranks r y , where σ y is the standard deviation of y and where σ ry is the standard deviation of r y .

Following [START_REF] Lerman | A note on the calculation and interpretation of the Gini index[END_REF] and under the assumption of a linear Pen's parade (i.e. [START_REF] Milanovic | A simple way to calculate the Gini coecient, and some implications[END_REF] demonstrates that for a suciently large n,

y = a + b r y ),
the Gini index can be further expressed as:

G = σ y √ 3y ρ(y, r y ) . (3)
Milanovic's result is very interesting since it yields a simple way to compute the Gini index. However, as mentioned by Milanovic (1997, page 48) himself, "in almost all real world cases, Pen's parade is convex: incomes at rst rise very slowly, and then their absolute increase, and nally even the rate of increase, accelerates". Thus, ρ(y, r y ) which measures linear correlation will be less than

y = k-1 i=0 b i r α i + b k r α k . ( 4
)
with k ∈ N * , α 0 = 0 and α i ∈ R * + for i = 1, . . . , k, and

α k > max{α i , i = 1, . . . , k -1}.
The covariance between y and r α j for j ∈ N * is given by: cov(y,

r α j ) = k i=1 b i cov ( r α i , r α j ) (5) = k i=1 b i n r α i +α j +1 - b i n 2 n r=1 r α i n r=1 r α i = k i=1 b i n n r=1 r α i +α j +1 - b i n 2 n r=1 r α i n r=1 r α j
The mean income ȳ is then:

y = k i=0 b i r α i = 1 n n r=1 k i=0 b i r α i (6)

The coecient of variation and Gini computation

Since incomes y are positive, we use (13) by assuming that b k > 0 and b j for j = 1, . . . , k -1 are chosen such that y > 0. For instance, if k = 2 then we can use b j for j = 1, . . . , k such that b 2 1 -4b 2 b 0 < 0. We are now able to compute the coecient of variation of incomes as follows:

σ y y = k i=1 k j=1 b i b j cov(r α i , r α j ) k i=1 b i r α i (7) = k i=1 b i k j=1 b j cov(r α i , r α j ) k i=1 b i r α i (8) = k j=1 b j cov(y, r j y ) k i=1 b i r α i (9) = k j=1 b j k i=1 b i n n r=1 r α i +α j +1 y -b i n 2 n ry=1 r α i n r=1 r α j 1 n n r=1 k i=0 b i r α i .
where the variance of

r α k is cov(r α k , r α k ) = σ 2 r α k = 1 n n r=1 r 2α k - 1 n n r=1 r α k 2 , (10)
the covariance between r α j and r α i for 0 < i, j ≤ k, is:

cov (r α i , r α j ) = 1 n n r=1 r α i +α j - 1 n 2 n r=1 r α i n r=1 r α j . ( 11
)
After a double summation

cov k i=1 b i r α i , k j=1 b j r α j = k i=1 k i=1 b i b j 1 n n r=1 r α i +α j - 1 n 2 n r=1 r i n r=1 r α j . (12)
Lemma 2.1 When n → ∞, for q ∈ R * + and r ∈ N * , we have that n r=1 r q ≡ n q+1 q + 1 .

(13)

Based on the preceding lemma, the variance of y when n → ∞ is equivalent to:

cov k i=1 r α i , k j=1 r α j ≡ b 2 k n n α k +α + +1 α k + α k + 1 - b 2 k n 2 n α k +1 α k + 1 n α k +1 α k + 1 ≡ (n α k b k k) 2 (2α k + 1)(α k + 1) 2 .
(14) Therefore when n → ∞ the standard deviation of y is equivalent to

σ y ≡ (n α k b k k) 2 (2α k + 1)(α k + 1) 2 ≡ α k |b k | (α k + 1) √ 2α k + 1 n α k . ( 15
)
When n → ∞, the mean of y is equivalent to

y = 1 n n r=1 k i=0 b i r α i ≡ b k n n α k +1 α k + 1 ≡ b k n α k α k + 1 (16)
Thereby, as n → ∞ the coecient of variation is equivalent expressed as:

σ y y ≡ α k |b k | (α k +1) √ 2α k +1 n α k b k n α k α k +1
(17) then we have the following limit which depends on k and the sign of b k :

lim n→∞ σ y y = |b k | b k α k √ 2α k + 1 = sign(b k ) α k √ 2α k + 1 . ( 18
)
We have then proved the following theorem:

Theorem 2.1 Under the assumption of a power Pen's parade, i.e., y = k i=0 b i r i y , with b k = 0, when n → ∞, the coecient of variation of the revenue y has the following limit:

lim n→∞ σ y y = |b k | b k α k √ 2α k + 1 (19)
On the other hand, following [START_REF] Milanovic | A simple way to calculate the Gini coecient, and some implications[END_REF]:

lim n→∞ 2 σ ry n = lim n→∞ n 2 -1 3n 2 = 1 √ 3 . (20)
Remark 2.1 The following theorem provide an interesting result in practise.

It tells us that, for a given income data such as the coecient of variation, we can compute the degree of the power Pen's parade. Therefore, we can compute the Pen's parade.

Theorem 2.2 For b k > 0, If we know the coecient of variation CV emp by using incomes data, then we can nd the parameter s = α k as a solution of the following equation:

s √ 2 s + 1 = CV emp . (21)
We can then use multiple regression to nd b i and get the Pen's parade that corresponds to incomes data. After a straightforward calculus, the positive solution of the equation ( 25) is :

α k = CV emp CV emp + CV 2 emp + 1 . (22)
The product of (20), (18) and ρ(y, r y ) entails the following result:

Theorem 2.3 Under the assumption of a power Pen's parade, i.e., y = k i=0 b i r α i , with b k = 0, when n → ∞, the Gini index G k can be approximately computed as follows:

G α k 1 √ 3 |b k | b k α k √ 2α k + 1 ρ(y, r y ), (23)
where ρ(y, r y ) is correlation coecient between the incomes y and their ranks r y .

Corollary 2.1 For b k > 0, the Gini index is approximately equal to

G α k 1 √ 3 |b k | b k α k √ 2α k + 1 ρ(y, r y ), (24)
3 Pen's parade of dierent countries data and their Gini indexes

In this section, we use the [START_REF] Milanovic | A simple way to calculate the Gini coecient, and some implications[END_REF] countries incomes data for the computation of the Pen's parade and gini index of each country income data.

In relation with the paramater k and b k > 0, we have the following table:

Concluding Remarks

In this paper, we proposed the power pen's parade as an alternative to the linear and quadratic pen's parade. Following [START_REF] Milanovic | A simple way to calculate the Gini coecient, and some implications[END_REF], we have proposed another simple way to calculate the Gini coecient under the assumption of a general power Pen's parade of order k . By using the data of Table 2, we concluded that the computation of the Gini index of each income data country need to nd a specic α k ∈ R * + which is the order of a specic power Pen's parade. We have therefore compute the Pen's parade and the gini index for each of the six countries.

It appeared that for a given coecient of variation of a country incomes data, we found the associated degree of the power Pen's parade, hence enabling compute the pen's parade using multiple regression. It also provides a convenient and straightforward way to compute the Gini index of each country by using the coecient of correlation between incomes and ranks.

One immediate and practical implications result from this new Gini expression. Estimating the coecients b i for i = 1, . . . , k, e.g. with Yitzhaki's Gini regression analysis or a multiple regression, enables a parametric Gini index to be obtained that depends on parameters reecting the curvature of Pen's parade, which may be of interest when one compares the shape of two income distributions.

Remark 3 . 2

 32 In our power Pen's parade, if r y = 1, then y = y 1 = k i=0 b k . In practical applications, the parameters b k which are estimated from the observed data using multiple regression (OLS) are such that the revenue y > 0 (i.e. b k > 0). It very important to note that, the correlation coecient between y and its rank r y depend on b k . Ignoring this fact equivalent to arbitrarily restricting the parade to pass through the origin and may result in less accurate estimates of the Gini index.

Table 1 :

 1 The Positive solution s = α k (see (22)) of equation (25) in relation with the empirical coecient of variation CV emp of each country.

	CV emp 0.43 0.56 0.57 0.60 0.68 0.68 0.71 0.76 1.07 1.63
	α k	0.653 0.955 0.981 1.060 1.285 1.285 1.375 1.532 2.712 5.774

where α k is the solution of the following equation:

where CV emp is the coecient of variation obtained by using incomes data y.

Remark 2.2 Note that for α k = 1, we have the result of [START_REF] Milanovic | A simple way to calculate the Gini coecient, and some implications[END_REF],

i.e.

Correlation between the incomes and ranks

The covariance between income y and the rank r y is

where

is an harmonic number function for s > 1 and n ∈ N * . In the simple case where, y = b 0 + b 1 r α 1 , we have:

The variance of the income y is:

The standard deviation of income is:

The standard deviation deviation of individuals ranks is given by:

(30)

The correlation coecient between y and r is:

b k > 0), then the Gini index is approximately equal to

where

with CV emp being the coecient of variation obtained by using income data y.


M Battaglia Parodi

P Iacono

G Ravalico

MD Francesco Bandello
email: francesco.bandello@hsr.it

Subthreshold Laser Treatment for Retinal Arterial Macroaneurysm

Purpose: To assess the effects of subthreshold laser treatment(STLT) for retinal arterial macroaneurysms(RAM) associated with foveal exudative manifestations and visual acuity deterioration.

Methods: Patients with RAM associated with foveal exudative manifestations and best corrected visual acuity (BCVA) worse than 20/80 Snellen equivalent underwent a ophthalmologic examination, including ETDRS visual acuity, optical coherence tomography (OCT), and fluorescein angiography. The patients were prospectively observed for 4 months period, and in absence of spontaneous improvement, they underwent STLT using an infrared diode laser.

Primary outcome measures were reduction of mean central point thickness(CPT) and BCVA changes at the 12-month examination. Secondary outcomes included changes in mean total macular volume(TMV) and central subfield thickness(CST).

Results: Nine patients were enrolled and prospectively followed up. Mean baseline values of BCVA, CPT, TMV, and CST were 0.8±0.1 (logMAR±SD), 340±49µm, 7.14±0.05mm 3 , and 366±37μm, respectively. At the 4-month examination following STLT, mean BCVA improved to 0.6±0.2, whereas mean CPT, TMV, and CST decreased to 274±29µm, 6.87±0.11mm 3 , and 296±33μm. At the 12-month examination, mean BCVA was 0.36±0.2, mean CPT was 195±11µm, mean TMV was 6.55±0.19mm 3 , and mean CST was 239±14μm, respectively. No side-effects were noted. In particular, no sign of retinal thinning and underlying backscattering typical of conventional laser treatment could be detected at the site of the laser application on OCT.

Conclusion:

The current pilot investigation of STLT for the treatment of symptomatic RAM shows encouraging data. A randomized clinical trial is required to ascertain the real efficacy this technique and the most appropriate settings to be employed.

Introduction

Retinal Arterial Macroaneurysms (MAR) are acquired saccular or fusiform dilations of the large arterioles of the retina. MAR are strongly associated with systemic hypertension and atherosclerosis. [START_REF] Nadel | Macroaneurysms of the retinal arteries[END_REF][START_REF] Lewis | Acquired arterial macroaneurysms of the retina[END_REF][START_REF] Mn | Retinal macroaneurysm: natural history and guidelines for treatment[END_REF][START_REF] Lavin | Retinal arterial macroaneurysms: a retrospective study of 40 patients[END_REF][START_REF] Rabb | Retinal arterial macroaneurysms[END_REF][START_REF] Panton | Retinal arterial macroaneurysms: Risk factors and natural history[END_REF][START_REF] Brown | Retinal arterial macroaneurysm:long term visual outcome[END_REF][START_REF] Moosavi | Retinal artery macroaneurysms: clinical and fluorescein angiographic features in 34 patients[END_REF] Natural history may lead to a spontaneous obliteration of the lesion, but also to a number of vision threatening complications, including subretinal, preretinal, or vitreous hemorrhages, macular edema, serous macular detachment, macular deposition of hard exudates, macular hole, and branch retinal vein occlusion. [START_REF] Nadel | Macroaneurysms of the retinal arteries[END_REF][START_REF] Lewis | Acquired arterial macroaneurysms of the retina[END_REF][START_REF] Mn | Retinal macroaneurysm: natural history and guidelines for treatment[END_REF][START_REF] Lavin | Retinal arterial macroaneurysms: a retrospective study of 40 patients[END_REF][START_REF] Rabb | Retinal arterial macroaneurysms[END_REF][START_REF] Panton | Retinal arterial macroaneurysms: Risk factors and natural history[END_REF][START_REF] Brown | Retinal arterial macroaneurysm:long term visual outcome[END_REF][START_REF] Moosavi | Retinal artery macroaneurysms: clinical and fluorescein angiographic features in 34 patients[END_REF][START_REF] Tsujikawa | Retinal structural changes associated with retinal arterial macroaneurysm examined with optical coherence tomography[END_REF][START_REF] Russell | Branch retinal artery occlusion after dye yellow photocoagulation of an arterial macroaneurysm[END_REF][START_REF] Ciardella | Ruptured retinal arterial macroaneurysm associated with a stage IV macular hole[END_REF][START_REF] Colucciello | Macular hole following ruptured retinal arterial macroaneurysm[END_REF][START_REF] Sagara | Macular hole formation after macular hemorrhage associated with rupture of retinal arterial macroaneurysm[END_REF][START_REF] Parodi | Branch retinal vein occlusion after spontaneous obliteration of retinal arterial macroaneurysm[END_REF][START_REF] Takahashi | Serous macular detachment associated with retinal arterial macroaneurysm[END_REF] Treatment of MAR is still controversial. [START_REF] Nadel | Macroaneurysms of the retinal arteries[END_REF][START_REF] Lewis | Acquired arterial macroaneurysms of the retina[END_REF][START_REF] Mn | Retinal macroaneurysm: natural history and guidelines for treatment[END_REF][START_REF] Lavin | Retinal arterial macroaneurysms: a retrospective study of 40 patients[END_REF][START_REF] Rabb | Retinal arterial macroaneurysms[END_REF][START_REF] Panton | Retinal arterial macroaneurysms: Risk factors and natural history[END_REF][START_REF] Brown | Retinal arterial macroaneurysm:long term visual outcome[END_REF][START_REF] Moosavi | Retinal artery macroaneurysms: clinical and fluorescein angiographic features in 34 patients[END_REF][START_REF] Tsujikawa | Retinal structural changes associated with retinal arterial macroaneurysm examined with optical coherence tomography[END_REF] In general, treatment is reserved to MAR causing exudative manifestations involving the fovea, and it is carried out using threshold laser photocoagulation in an attempt to achieve the obliteration of the vascular lesion, which leads to visible scar. Nevertheless, several complications have been related to laser photocoagulation, including arteriolar obliteration and consequent retinal ischemia. [START_REF] Russell | Branch retinal artery occlusion after dye yellow photocoagulation of an arterial macroaneurysm[END_REF] Subthreshold laser treatment (STLT) has been proposed in an attempt to minimize the adverse aspects of conventional laser photocoagulation by reducing the duration of laser exposure and by using a subvisible clinical end point. [START_REF] Friberg | The treatment of macular disease using a micropulsed and continuous wave 810-nm diode laser[END_REF][START_REF] Roider | Retinal sparing by selective retinal pigment epithelial photocoagulation[END_REF][START_REF] Moorman | Clinical applications of the MicroPulse diode laser[END_REF][START_REF] Parodi | Subthreshold grid laser treatment of macular edema secondary to branch retinal vein occlusion with micropulse infrared (810 nanometer) diode laser[END_REF][START_REF] Parodi | Intravitreal triamcinolone acetonide combined with subthreshold grid laser treatment for macular ooedema in branch retinal vein occlusion: a pilot study[END_REF] Repeated laser pulses potentially localises the thermal effects to the retinal pigment epithelium (RPE) with limited protein denaturation. [START_REF] Desmettre | Micropulse and continuous wave diode retinal photocoagulation: Visible and subvisible lesion parameters[END_REF] The aim of the present pilot study is to evaluate the effect of subthreshold laser treatment (STLT) in reducing the exudative manifestations secondary to RAM. Inclusion criteria included the diagnosis of RAM, regarded as symptomatic whenever associated with visual acuity deterioration and distortion. More specifically, we selected eyes presenting RAM, defined as an acquired dilation of the retinal arteries on biomicroscopy and fluorescein angiography (FA), which was associated with exudative manifestations involving the fovea (including: subretinal/intraretinal fluid, and/or hard exudates, and/or hemorrhages), and presented best corrected visual acuity (BCVA) worse than 20/80 Snellen equivalent, corresponding to 0.6 LogMAR. Patients who had undergone any other treatment for RAM in the past, or were affected by any other ocular disease were excluded from the study.

Methods

Each patient underwent a complete ophthalmological examination, including BCVA, on standard ETDRS charts, FA and optical coherence tomography (OCT). The patients were regularly scheduled at each 4-month follow-up for up to 12 months after laser treatment. BCVA measurement and OCT were performed by masked staff members, who were unaware of the purpoe of the study.

OCT (Stratus, Carl Zeiss, Meditec Inc, Germany) was performed using the fast macular thickness map protocol, consisting of six 6-mm radial lines, oriented 30° apart, centred on the fovea. When the automatic segmentation algorithm of fast macular thickness map failed to identify inner and outer retinal boundaries correctly, the central point thickness (CPT) was manually measured for each radial scan and the average of the 6 values obtained was used for statistical analysis. The mean CPT was defined as the distance from the inner retinal surface to the inner border of the RPE. Total macular volume (TMV) was determined by the sum of the volumes of the 9 quadrants obtained by a central macular thickness map. Central subfield thickness (CST) was defined as the mean value of the thickness values obtained in the central subfield.

In accordance with the protocol, each patient was first observed for a 4 months period to allow any spontaneous improvement of the symptoms to occur. If no improvement in BCVA and CPF was detected after 4 months, the patient was considered not able to achieve a spontaneous recovery and thus, eligible for the treatment with STLT, in an attempt to limit the exudation from the RAM, and to improve the visual function. Our protocol also made provision for each patient to undergo conventional visible laser treatment if BCVA had no improved by at least 1 ETDRS line and/or CPT had not reduced at least by 10% after 4 months following STLT.

STLT was performed using an infrared diode laser (Iris Medical Oculight, SLx Photocoagulator, Iridex Corp, Mountain View, Ca). The laser used parameters were: 125μm laser spot diameter, 0.3 second exposure, and 15% duty cycle. The laser power was arbitrarily set at 1400mW. This setting was chosen on the basis of our previous experience with STLT in the treatment of other diseases [START_REF] Parodi | Subthreshold grid laser treatment of macular edema secondary to branch retinal vein occlusion with micropulse infrared (810 nanometer) diode laser[END_REF][START_REF] Parodi | Intravitreal triamcinolone acetonide combined with subthreshold grid laser treatment for macular ooedema in branch retinal vein occlusion: a pilot study[END_REF] and the specific treatment of RAM (Parodi MB unpublished data). In particular, from our previous data regarding the treatment of RAM, we found out that a lower laser power had no beneficial effect on the leakage from RAM. In each case, laser spots were delivered in a contiguous mode with no free space between each spot application, and covered both the entire RAM lesion and the surrounding area, in order to promote a response in the neighbouring RPE too.

Primary outcome measures were the reduction of CPT and the improvement of BCVA at the 12month examination. Secondary outcomes included also the changes in TMV and CST.

Results

Twelve consecutive patients affected by RAM with exudative manifestations involving the fovea were considered for the study. Nine patients who fulfilled the eligibility and exclusion criteria were enrolled (6 females and 3 males), with a mean age of 74.2 years (range: 56-81, SD: 11.3). Three patients were excluded because they had previously undergone to laser treatment.

Of the 9 patients enrolled, 7 were affected by hypertension and 2 by diabetes mellitus. At the initial examination, 6 eyes showed a macular hemorrhage centered on the fovea, comprising a subinternal limiting membrane hemorrhage associated with a retinal hemorrhage. The remaining 3 eyes revealed on OCT a retinal edema especially at the outer retinal layer, associated with serous retinal detachment involving the fovea. FA showed profuse leakage from RAM in all cases.

All patients performed all the planned controls and completed the 12-month follow-up.

The complete data of patients is listed in Table 1. At the enrolment, all the 9 patients had a BCVA ≤ 0.6 LogMAR (approximately 20/80 Snellen equivalent), with a mean BCVA was 0.8±0.1 SD LogMAR (approximately 20/125 Snellen equivalent) (range: 1-0.6 LogMAR), a mean CPT of 314±22μ SD, a mean TMV of 7.05 mm 3 ±0.06 SD, and a mean CST of 344±27μ SD.

In accordance with the protocol, the 9 patients enrolled were simply observed for 4 months to permit a possible improvement in BCVA to take place and to monitor the spontaneous development of the clinical picture. After the 4-month observation period no patient showed any BCVA improvement. On the contrary CPT was increased in all the 9 patients with a mean value of 340±49μ SD, TMV had showed a mean value of 7.14 mm 3 ±0.05 SD, and mean CST was 366±37μ SD. Thus, all the patient were considered eligible to undergo STLT.

At the 4-month examination following STLT, BCVA and CPT improved in all 9 patients. More specifically, mean BCVA was 0.6±0.2 SD LogMAR (approximately 20/80 Snellen equivalent) (range: 0.9-0.4 LogMAR), CPT reduced to 274±29μ SD, TMV to 6.87 mm 3 ±0.11 SD, and CST to 296±33μ SD.

Mean BCVA readings at the 8-and 12-month examinations were to 0.5±0.1 SD LogMAR (approximately 20/63 Snellen equivalent) (range: 0.3-0.7 LogMAR), and 0.36±0.2 SD LogMAR (approximately 20/50 Snellen equivalent) (range: 0.2-0.6 LogMAR), respectively (Fig. 1). Mean CPT reduced to 221±14μ SD and 195±11μ SD at the 8-and the 12-month control examinations, mean TMV passed to 6.68 mm 3 ±0.2 SD and 6.55 mm 3 ±0.19 SD, and mean CST was 257±20μ SD and 239±14μ SD, respectively (Figure 1 and2).

The frequency distribution of BCVA changes revealed that all the patients gained a mean of 4.4 lines at the end of the follow-up, with 8 out of 9 eyes gaining 3 or more lines. No patient required additional treatment over the follow-up.

Baseline OCT performed at the level of the RAM revealed a vascular structure with a reflective wall associated with retinal thickening in the surrounding area. OCT scans after STLT revealed a reduction of RAM size associated with disappearance of the intraretinal oedema, with no sign of the retinal thinning and underlying backscattering, typical of conventional laser treatment.

CPT was manually measured in six case, because the automatic segmentation algorithm failed to precisely identify inner and outer retinal boundaries.

All cases showed a complete resolution of exudative manifestation at the end of the follow-up and a cessation of the dye leakage on FA. Interestingly, no biomicroscopic, or angiographic evidence of the treatment turned out to be detectable at the laser application site (Fig. 3). No side-effect was noted over the follow-up.

Discussion

The treatment of RAM is controversial and at present there is no general consensus about the optimum approach. Spontaneous obliteration of RAM may occur over the follow-up, but the long persistence of blood or subretinal fluid may lead to progressive photoreceptor damage, with consequent irreversible visual impairment. [START_REF] Tsujikawa | Retinal structural changes associated with retinal arterial macroaneurysm examined with optical coherence tomography[END_REF] More specifically, the presence of subretinal hemorrhage beneath the fovea leads to a diminished visual function. [START_REF] Tsujikawa | Retinal structural changes associated with retinal arterial macroaneurysm examined with optical coherence tomography[END_REF] Treatment is generally indicated for RAM associated with exudative manifestations involving the fovea which do not tend to resolve spontaneously. The most commonly employed approach is direct or indirect laser photocoagulation with visible end-point. Nevertheless, several complications may be associated with conventional threshold laser for macular lesions, including enlargement of laser scar, choroidal neovascularization, and subretinal fibrosis. [START_REF] Schatz | Progressive enlargement of laser scars following grid laser photocoagulation for diffuse diabetic macular edema[END_REF][START_REF] Morgan | Atrophic creep of the retinal pigment epithelium after focal macular photocoagulation[END_REF][START_REF] Treatment | Focal photocoagulation treatment of diabetic macular edema. Relationship of treatment effect to fluorescein angiographic and other retinal characteristics at baseline: ETDRS report no 19[END_REF][START_REF] Lewen | Subretinal neovascularization complicating laser photocoagulation of diabetic maculopathy[END_REF][START_REF] Lewis | Choroidal neovascularization after laser photocoagulation for diabetic macular edema[END_REF][START_REF] Guyer | Subretinal fibrosis after laser photocoagulation for diabetic macular edema[END_REF][START_REF] Rutledge | Sub-pigment epithelial membranes after photocoagulation for diabetic macular edema[END_REF][START_REF] Fong | Subretinal fibrosis in diabetic macular edema. ETDRS report 23. Early Treatment Diabetic Retinopathy Study Research Group[END_REF][START_REF] Han | Submacular fibrosis after photocoagulation for diabetic macular edema[END_REF] In addition to these complications, arteriolar obliteration, increased retinal exudation and scarring, with possible retinal traction, have also been reported as possible sequelae of the laser photocoagulation of RAM. [START_REF] Lewis | Acquired arterial macroaneurysms of the retina[END_REF][START_REF] Mn | Retinal macroaneurysm: natural history and guidelines for treatment[END_REF][START_REF] Lavin | Retinal arterial macroaneurysms: a retrospective study of 40 patients[END_REF][START_REF] Rabb | Retinal arterial macroaneurysms[END_REF][START_REF] Panton | Retinal arterial macroaneurysms: Risk factors and natural history[END_REF][START_REF] Brown | Retinal arterial macroaneurysm:long term visual outcome[END_REF][START_REF] Moosavi | Retinal artery macroaneurysms: clinical and fluorescein angiographic features in 34 patients[END_REF][START_REF] Tsujikawa | Retinal structural changes associated with retinal arterial macroaneurysm examined with optical coherence tomography[END_REF][START_REF] Russell | Branch retinal artery occlusion after dye yellow photocoagulation of an arterial macroaneurysm[END_REF] The use of STLT represents a different approach to the treatment of RAM. STLT has been proposed in an attempt to minimize the negative aspects of conventional grid laser photocoagulation and has achieved encouraging results, especially in macular oedema secondary to diabetic retinopathy, proliferative diabetic retinopathy, and branch retinal vein occlusion. [START_REF] Parodi | Subthreshold grid laser treatment of macular edema secondary to branch retinal vein occlusion with micropulse infrared (810 nanometer) diode laser[END_REF][START_REF] Parodi | Intravitreal triamcinolone acetonide combined with subthreshold grid laser treatment for macular ooedema in branch retinal vein occlusion: a pilot study[END_REF][START_REF] Friberg | Infrared micropulsed laser treatment for diabetic macular edema--subthreshold versus threshold lesions[END_REF][START_REF] Luttrull | Subthreshold diode laser micropulse photocoagulation for the treatment of clinically significant diabetic macular oedema[END_REF][START_REF] Luttrull | Serial optical coherence tomography of subthreshold diode laser micropulse photocoagulation for diabetic macular edema[END_REF][START_REF] Figueira | Prospective randomized controller trial comparing subthreshold micropulse diode laser photocoagulation and conventional green laser for clinically significant diabetic macular oedema[END_REF][START_REF] Luttrull | Subthreshold diode micropulse panretinal photocoagulation for proliferative diabetic retinopathy[END_REF] The purpose of employing STLT is to attempt to decrease laser damage by reducing the duration of laser exposure and by using a subvisible clinical end point. A micropulse diode laser releases low energy micropulses in order to confine the energy to the RPE cells, with limited protein denaturation, [START_REF] Desmettre | Micropulse and continuous wave diode retinal photocoagulation: Visible and subvisible lesion parameters[END_REF] avoiding lateral thermal spreading.

The current study presents the result of a pilot investigation of STLT for the treatment of symptomatic RAM show encouraging data. All the patients achieved a cessation of leakage on FA, with improvements in BCVA, CPT, CST, and TMV. No biomicroscopic, angiographic, or OCT sign of laser treatment was detected over the follow-up. Moreover, no evidence of arterial obliteration, increased retinal exudation, or retinal traction was noted.

The underlying mechanisms by which STLT acts remains a matter of speculation, because no specific histopathologic study in humans has been carried out. However, bearing in mind the selective damage to the RPE cells induced by a subthreshold laser application, [START_REF] Chong | Selective RPE damage by micro-pulse diode laser photocoagulation[END_REF][START_REF] Chong | A retinal laser which damages only the RPE: Ultrastructural study[END_REF][START_REF] Morimura | Histological effect and protein expression in subthreshold transpupillary thermotherapy in rabbit eyes[END_REF] it can be hypothesized that the effects of STLT may lead to an improved balance in cytokines release, perhaps including an upregulation of basic fibroblast growth factor. [START_REF] Xiao | Effects of retinal laser photocoagulation on photoreceptor basic fibroblast growth factor and survival[END_REF] Thus, STLT may act on the complex interaction among the different molecules implicated in the vascular permeability, promoting a vascular repair together with an improved fluid reabsorption.

We are perfectly aware that ours is a pilot study designed only to find out if a positive effect can be obtained using STLT in the treatment of long-standing symptomatyc RAM. The major limitations of the present study are represented by the small number of patients involved and the absence of a control group. In an attempt to be more confident that the effect noted are indeed due to STLT and not to the spontaneous evolution of RAM as part of its natural history, we designed the protocol so as to turn to STLT only after a 4-month period of simple observation in order to allow for the possible spontaneous resolution of the exudative manifestations of RAM. Patients were then there treated, as it was consider improper to prolong the observation beyond this period in patients having low visual acuity. Moreover, we chose the laser setting on the basis of our personal experience and cannot rule out that a different setting might lead to a even better functional and anatomical outcome. Legend: Mean central point thickness (μm) over the follow-up. A statistically significant difference was found (p<0.01, Student's t-test) from the 4-month examination after subthreshold laser treatment on.

  Patients referred to Department of Ophthalmology of the University of Udine and Trieste between January 2007 and January 2008 with the diagnosis of RAM were prospectively enrolled in the study. The research adhered to the tenets of the Declaration of Helsinki, and investigational review board approval was obtained. Each patient was carefully informed about the purpose of the research, providing signed consent to all procedures.

Figure and Table Legends

  Figure and Table Legends

Figure 1 Mean

 1 Figure 1

Figure 2 Mean

 2 Figure 2

Figure 3 Case 8 .

 38 Figure 3

Fig. 1 .

 1 Fig. 1. Mean Best corrected visual acuity values (LogMAR) over the follow-up.

Figure 2 .

 2 Figure 2. Mean central point thickness (μm) over the follow-up.

Table legend BCVA

 legend indicates best corrected visual acuity; CPT indicates central point thickness; CST indicates central subfield thickness; TMV indicates total macula volume. Statistical analyses were performed using Wilcoxon test.

Table 1 .

 1 Best corrected visual acuity and OCT data of the nine patients treated with subthreshold laser treatment over the follow-up Legend: BCVA indicates best corrected visual acuity; CPT indicates central point thickness; CST indicates central subfield thickness; TMV indicates total macula volume. Statistical analyses were performed using Wilcoxon test.

				BCVA					CPT					CST					TMV
	CASE baseline observation 4 mos 8 mos 12 mos baseline observation 4 mos 8 mos 12 mos baseline observation 4 mos 8 mos 12 mos baseline observation 4 mos 8 mos	12 mos
	1	1,0	1,0	0,9	0,7	0,6	345	450	257	230	213	361	436	263	241	229	7,01	7,20	6,78	6,80	6,70
	2	0,8	0,7	0,7	0,6	0,5	289	299	222	220	189	297	308	253	251	210	6,99	7,08	6,80	6,80	6,59
	3	0,7	0,7	0,7	0,6	0,5	278	270	240	231	187	299	319	267	245	230	7,03	7,12	6,85	6,90	6,78
	4	0,9	0,9	0,5	0,4	0,3	310	327	289	201	198	346	367	301	236	249	7,07	7,20	7,01	6,90	6,80
	5	0,6	0,6	0,4	0,4	0,3	340	329	296	212	201	367	359	341	278	256	7,18	7,10	7,00	6,50	6,50
	6	0,8	0,9	0,4	0,4	0,2	324	350	302	231	212	356	373	349	289	248	7,01	7,13	6,80	6,40	6,30
	7	0,7	0,8	0,5	0,3	0,2	301	328	289	201	178	356	367	301	234	243	6,99	7,08	6,90	6,40	6,30
	8	0,9	0,9	0,7	0,5	0,4	325	355	301	243	189	367	378	302	267	237	7,12	7,20	7,00	6,70	6,42
	9	0,8	0,8	0,6	0,5	0,2	314	360	269	222	190	354	387	289	278	249	7,06	7,14	6,73	6,70	6,53
	media	0,8	0,8	0,6	0,5	0,4	314	341	274	221	195	345	366	296	258	239	7,05	7,14	6,87	6,68	6,55
	sd	0,1	0,1	0,2	0,1	0,2	22,3	49,8	28,8	14,4	11,8	27,3	37,3	33,1	20,6	14,2	0,06	0,05	0,11	0,20	0,19
	P
	value			0.016	0.003 0,003			0,003	0,003	0,003			0,003 0,003	0,003			0,003 0,003	0,003

Thus, this research is to be considered an initial investigation acting as the basis for a multicentre randomised clinical trial. In essence, STLT offers promising results in the treatment of symptomatic RAM without bringing about negative side-effects. A multicentre randomized clinical trial is required to ascertain the real efficacy of this technique, and the most appropriate settings to be employed.


Eduardo Cuesta

Mokhtar Kirane

Salman A Malik

Image structure preserving denoising using generalized fractional time integrals

Keywords: Image processing, Fractional integrals and derivatives, Volterra equations, Convolution quadrature methods. 2010 MSC: 44A35, 44K05, 45D05, 65R20, 68U10, 94A08

A generalization of the linear fractional integral equation u(t) = u 0 + ∂ -α Au(t), 1 < α < 2, which is written as a Volterra matrix-valued equation when applied as a pixel-by-pixel technique, has been proposed for image denoising (restoration, smoothing,...). Since the fractional integral equation interpolates a linear parabolic equation and a hyperbolic equation, the solution enjoys intermediate properties. The Volterra equation we propose is well-posed for all t > 0, and allows us to handle the diffusion by means of a viscosity parameter instead of introducing non linearities in the equation as in the Perona-Malik and alike approaches. Several experiments showing the improvements achieved by our approach are provided.

Introduction

Partial differential equations based methods for image processing (filtering, denoising, restorations, segmentation, edge enhancement/detection,...) have been largely studied in the literature (see [START_REF] Weickert | Anisotropic Diffusion in Image Processing[END_REF] and references therein).

In that framework the first, and most investigated equation is the linear heat equation with homogeneous Neumann boundary condition

             ∂ t u(t, x) = ∆u(t, x), (t, x) ∈ [0, T ] × Ω, u(0, x) = u 0 (x), x ∈ Ω, ∂u ∂η (t, x) = 0, (t, x) ∈ [0, T ] × ∂Ω, (1)
where ∂ t , and ∆ stand for the time derivative, and twodimensional Laplacian operator, respectively, Ω ⊂ R 2 is typically a square domain, ∂Ω represents the boundary of Ω, ∂/∂η stands for the outward normal derivative, and u 0 is the noisy image from which the objective is to restore the original (ideal) image. Let us notice that u(t, x) stands for the restored image at the time level t.

Email addresses: eduardo@mat.uva.es (Eduardo Cuesta), mokhtar.kirane@univ-lr.fr (Mokhtar Kirane), salman.malik@univ-lr.fr (Salman A. Malik)

The interest for this model is due to the fact that the solution of (1) can be written as a convolution

u(t, x) = R 2 G √ 2t (x -y)u 0 (y) dy,
where G is the two-dimensional Gaussian kernel

G σ (x) := 1 2πσ 2 e -|x| 2 /2σ 2 .
Since convolution with a positive kernel is the basic tool in linear filtering, computing the solution of ( 1) is equivalent to Gaussian filtering in a classical way.

However, in this equation the diffusion is isotropic which, in the context of image processing, means that smoothing applies uniformly in the whole image and is equivalent in the direction of both coordinates axis, therefore independently of the image structure itself. This yields that in most of cases edges and corners are severely blurred disabling this filter for practical applications.

In view of this, an edge-preserving regularization, i.e. a non uniform diffusion model, will require more sophisticated approaches. One of this was proposed by Perona and Malik in [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF] and reads as the non linear heat equation based problem

               ∂ t u(t, x) = div c(|∇u(t, x)| 2 )∇u(t, x) , u(0, x) = u 0 (x), ∂u ∂η (t, x) = 0, x ∈ ∂Ω, (2)
for (t, x) ∈ Q = [0, T ] × Ω. The diffusion coefficient c : [0, +∞) → [0, +∞), is chosen to be close to zero near edges and corners, that is, pixels where gradient is large. On the contrary, c should be large in pixels with low gradient variation. Functions satisfying these assumptions are commonly called edge stopping functions, and examples of such a function were firstly proposed by Perona and Malik, e.g,

c(s) = 1/(1 + s) or c(s) = e -s . (3)
Unfortunately, the edge stopping functions lead in general to backward-forward problems that turn out to be ill-posed. Numerical experiments carried out by some authors with these models (involving edge stopping functions) reported that no significant instabilities are observed; moreover, for large final times, restored images seem to preserve and enhance edges (see e.g. [START_REF] Hamza | Unifing probabilistic and variational estimation[END_REF]). The reason for not observing any instability in Perona-Malik model is that in many cases the numerical scheme used does not correspond to their equation but rather to a time-regularized one which is a well-posed equation as reported by H. Amman [START_REF] Amann | Time-delayed Perona-Malik type problems[END_REF]. Even in the case of the discretization corresponding to the equation (i.e. without regularizing) the only one artifact usually observed in the numerical experiments is the staircasing effect. This occurs when a standard spatial discretization of the equation is considered instead of the continuous one. In that case it has been proved in [START_REF] Solina | A semidicrete nonlinear scale-space theory and its relation to the Perona-Malik paradox[END_REF] that (2) becomes a well-posed system of nonlinear ordinary differential equations.

In the same way, some other approaches have been proposed as for example the ones based on the total variation of suitable functionals (see [START_REF] Rudin | Nonlinear total variation based noise removal algorithm[END_REF]). All these works have promoted the idea of replacing (2) by nearby equations keeping on the one hand the same image structure, and on the other hand, lying in a reasonable functionals space setting where the well-posedness can be guaranteed as well as the bounded variation, and further analytical and numerical properties. The first perturbed model was proposed in [START_REF] Catte | Image selective smoothing and edge detection by nonlinear diffusion[END_REF] where, for a suitable extension of u over R 2 (e.g. by 0) denoted ũ, c(|∇u|) is replaced by c(|∇(G σ * ũ)|) (G σ as defined above). In that case, for u 0 ∈ L 2 (Ω), the regularized problem admits a unique solution in C([0, T ], L 2 ) ∩ L 2 ((0, T ), H 1 ) where H 1 is the Sobolev space defined as

H 1 (R 2 ) := { f ∈ L 2 (R 2 ), ▽ f ∈ L 2 (R 2 )}.
Despite of some features of this model, images become uniformly grey (for grey-scale images) in the long run, thus the information gets lost (see [START_REF] Amann | Time-delayed Perona-Malik type problems[END_REF]). Variants of this approach have been studied by many authors (see e.g. [START_REF] Kacur | Solution of nonlinear diffusion appearing in image smoothing and edge detection[END_REF][START_REF] Kacur | Slow and fast diffusion effects in image processing[END_REF] and references therein).

Lately, further approaches have been proposed, e.g. by adding regularizing terms, like ǫ∆∂ t u, to the diffusion equation in [START_REF] Amann | Time-delayed Perona-Malik type problems[END_REF] (see [START_REF] Bellettini | A regularized Perona-Malik functional: some aspects of the gradient dynamics[END_REF]), or for example higher order partial differential equations based regularization (see [START_REF] Didas | Properties of higher order Figure 7: Denoising of textured images: (a) Original image of wood, (b) Noisy image perturbed by gaussian noise, (c) (PM), (d) (VE) (e) Original naive image, (f) Noisy image perturbed by gaussian noise, (g) (PM), (h) (VE). nonlinear diffusion filtering[END_REF]), and in particular fourth order partial differential equations (see [START_REF] Lysaker | Noise removal using fourth-order partial differential equation with applications to medical magnetic resonance images in space and time[END_REF][START_REF] You | Fourth-order partial differential equations for noise removal[END_REF]). However, despite of the practical results which seem to be quite good in most of the cases, some of them have not been closely studied from the analytical and the numerical view point.

Many image denoising models based on local average (also known as neighborhood filters) have been proposed for controlling the diffusion process. Let us mention some of the neighboring filters: Yaroslavsky neighborhood filter [START_REF] Yaroslavsky | Digital Picture Processing. An Introduction[END_REF], sigma filter [START_REF] Lee | Digital image smoothing and the sigma filter[END_REF], bilateral filtering [START_REF] Tomasi | Bilateral filtering for grey and color images[END_REF] are few examples. The faster diffusion is obtained in these filters by assigning appropriate weights where the neighborhood pixels have gray scale value close to one another (uniform region) and slower diffusion is obtained across the boundaries of the region, resulting in the preservation of edges in the image. All the local nonlinear filters (neighborhood filters) create artifact boundaries (staircase effect) in the restored images. In [START_REF] Buades | The staircasing effect in neighborhood filters and its solution[END_REF], the authors justify the phenomena of artifact boundaries by showing that the Yaroslavsky neighborhood filter has exactly the same qualitative behavior as Perona-Malik model [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF]. The similar behavior of staircase for the bilateral filter was observed (see [START_REF] Barash | A fundamental relationship between bilateral filtering, adaptive smoothing and the nonlinear diffusion equation[END_REF]).

Finally let us comment that one can face up to diffusion filtering in a more general framework as the one of anisotropic models. In fact most of models we mentioned above turn out to be isotropic since the diffusion is governed by means of a scalar-valued function c which allows to reflect the structure of the underlying image, but the diffusion turns out to be the same in the two orthogonal directions of coordinate axis. A generalization of such approaches consist in replacing c by a tensor which allows on the one hand to rotate the flux and so the diffusion orthogonal directions toward a suitable orientation, and on the other hand a suitable choice of the eigenvalues of the tensor can lead to different diffusions on each single direction (see e.g. [START_REF] Aja-Fernández | Tensors in Image Processing and Computer Vision[END_REF][START_REF] Bartels | Stable discretizations of scalar and contrained vectorial Perona-Malik equation[END_REF][START_REF] Weickert | A review of nonlinear diffusion filtering[END_REF][START_REF] Weickert | Coherence enhancing diffusion filtering[END_REF][START_REF] Weickert | Diffusion and regularization of vector-and matrix-valued images[END_REF] for more details). In the present work we focus on the scalar approach, i.e. no tensor diffusion are considered, keeping in mind that a tensorial approach in the framework of fractional calculus fits and will be the topic of our future works.

In this work, we present a new approach based on fractional calculus which allows us to handle the diffusion, i.e. the smoothing in the image terminology, by means of a single parameter which plays the role of "viscosity" parameter in a linear partial differential equation. Let us mention the new features of our proposed model although satisfactory practical results are obtained (see section 6). First, the linear model we are proposing here is well-posed; since the final objective is the practical implementation, very efficient numerical discretizations can be applied for the proposed model because a large variety of them have been closely studied by many authors and they are at our disposal for the experiments in Section 6. Moreover, the model we propose enjoys the well-posedness for all positive time levels, i.e. for t > 0, on the contrary to what typically happens when nonlinear equations are involved. In fact, despite of the Perona-Malik model turn out to be a well-posed systems of nonlinear equations after discretizing in space as we described earlier, the existence of a unique solution can be locally guaranteed i.e. for a finite time interval (normally small). Therefore, when using the semidiscrete Perona-Malik model in practical applications some estimation of the final time should be done. Second, we have to highlight the simplicity of our model for controlling the diffusion process, which has been achieved by choosing a scalar function which allows a suitable choice of the viscosity parameter as we describe in Section 5. Moreover, due to the numerical scheme chosen for our experiments, some additional properties are guaranteed, in fact the positivity of the image seen in the continuous setting is also guaranteed in the discrete one (Th. 1, in Section 4).

This paper is organized as follows. In section 2, we recall some facts concerning fractional calculus, and the first approaches to image processing by using fractional calculus. Section 3 focusses on generalized fractional integrals based approach to image processing which is the main novelty of this work. Numerical discretizations are presented in Section 4. The discussion on the implementation, and practical experiments, are shown in Sections 5 and 6 respectively. Finally, we end with conclusions in Section 7.

Fractional calculus

Image filtering by means of fractional calculus is first considered in [START_REF] Cuesta | Image processing by means of a linear integro-differential equation[END_REF]. In that work, a generalization of the heat equation ( 1) is proposed which reads as

                   ∂ α t u(t, x) = ∆u(t, x), (t, x) ∈ [0, T ] × Ω u(0, x) = u 0 (x), x ∈ Ω, ∂ t u(0, x) = 0, x ∈ Ω, ∂u ∂η (t, x) = 0, (t, x) ∈ [0, T ] × ∂Ω, (4)
where ∂ α t stands for the fractional time derivative of order 1 < α < 2 in the sense of Riemann-Liouville. Integrating in both sides, the problem (4) can be expressed as

u(t, x) = u 0 (x) + 1 Γ(α) t 0 (t -s) α-1 ∆u(s, x) ds, also with homogeneous Neumann boundary condition, or in a compact format u(t, x) = u 0 (x) + ∂ -α ∆u(t, x), (5)
where ∂ -β , for β > 0, stands for the fractional integral of order β ∈ R + , in the sense of Riemann-Liouville. Let us recall that, for g : [0, +∞) → R, g ∈ AC[0, +∞) (the space of absolutely continuous functions) the integral of order β ∈ R + in the sense of Riemann-Liouville is defined as the convolution integral

∂ -β g(t) := t 0 k β (t -s)g(s) ds, t ≥ 0, (6)
where k β (t) := t β-1 /Γ(β), for t > 0 (see [START_REF] Kilbas | Theory and Applications of Fractional Differential Equations[END_REF]). Now, the definition of the fractional derivative of order β ≥ 0 is

∂ β g(t) := d m dt m ∂ β-m g(t), t ≥ 0, where m ∈ N, m -1 < β ≤ m.
The interest of our model in the framework of image processing is due to the fact that, for 1 < α < 2, the problem (4) interpolates the linear (parabolic) heat equation ( 1) corresponding to α = 1 (with no need of ∂ t u(0, x) = 0), and the linear (hyperbolic) wave equation

                   ∂ 2 t u(t, x) = ∆u(t, x), (t, x) ∈ [0, T ] × Ω u(0, x) = u 0 (x), x ∈ Ω, ∂ t u(0, x) = 0, x ∈ Ω, ∂u ∂η (t, x) = 0, (t, x) ∈ [0, T ] × ∂Ω, (7)
corresponding to α = 2 (with zero initial velocity). Therefore, some properties of the solution of ( 4) are intermediate between the ones of ( 1) and [START_REF] Brunner | Collocation Methods for Volterra Integral and Related Functional Equations[END_REF] (see e.g. [START_REF] Fujita | Integro-differential equation which interpolates the heat equation and the wave equation[END_REF][START_REF] Fujita | Integro-differential equation which interpolates the heat equation and the wave equation (ii)[END_REF]). In particular, it can be easily proved that for the scalar equation (i.e., by replacing ∆ in (4) with a complex λ with non-positive real part) the solution decays as 1/(1 + t α |Re(λ)|), i.e., the solution is o(t -α ), as t → +∞. In that case, the solution decays slower than the one of the scalar heat equation (with the same λ) for which the solution decays exponentially, and faster than the one of the scalar wave equation whose solution does not decays (but oscillates). The diffusion is now handled by the parameter α. We show below that the solution of the model we propose is well behaved.

At the same time, fractional calculus was proposed for edge detection in [START_REF] Mathieu | Fractional differentiation for edge detection[END_REF], and later for image denoising in [START_REF] Bai | Fractional anisotropic diffusion for image denoising[END_REF]. In these papers, the authors proposed anisotropic equations, where the anisotropy is handled by means of spatial fractional derivatives; however, the papers do neither include the study of the well-posedness of the problem nor the study of the behavior of the fully numerical discretization.

Finally, let us mention [START_REF] Didas | Denoising and enhancement of digital images[END_REF][START_REF] Didas | Regularity and scale-space properties of fractional high order linear filtering[END_REF][START_REF] Didas | Discrete multiscale wavelet shrinkage and integrodifferential equations[END_REF][START_REF] Didas | Integrodifferential equations for continuous multiscale wavelet shrinkage: The discrete case[END_REF][START_REF] Didas | Integrodifferential equations for continuous multiscale wavelet shrinkage[END_REF] where the fractional calculus, also applied to image processing, is understood as the fractional powers of the twodimensional Laplacian, i.e. (-∆) β , for β > 0, therefore in a different framework and with some different features.

Volterra equations

Despite of the fact that the approach (5) in Section 2 seems to be very promising, the diffusion (smoothing) is still uniform all over the whole image as in [START_REF] Aja-Fernández | Tensors in Image Processing and Computer Vision[END_REF].

The above idea leaded us to a refinement which consists in splitting the whole image into sub-images, and apply (4) with different values of α for each sub-image (see [START_REF] Cuesta | Some advances on image processing by means of fractional calculus[END_REF]). However the approach we propose in this paper is more than a slight extension of the one in [START_REF] Cuesta | Some advances on image processing by means of fractional calculus[END_REF], from a technical viewpoint. Roughly speaking, in [START_REF] Cuesta | Some advances on image processing by means of fractional calculus[END_REF] the choice of each α is carried out by setting values close to 1 for the sub-images with lower mean gradient variation, and close to 2 for the sub-images with higher mean gradient variation. This approach does not provide satisfactory results as extended in the framework of image denoising despite of the fact that this approach has provided good results in some practical situations such as in satellite image classification (see [START_REF] Quintano | Improving satellite image classification by using fractional type convolution filtering[END_REF]). It is mainly due to the fact that in most of cases the borders of each single sub-image keep clearly sharpened which is undesirable in image denoising.

However this idea suggested us a finer approach which is the main contribution of this work, and which intends to be the limit of the above situation, i.e., the application of the fractional equation ( 4) with a different value of α (i.e. order of derivative) for each single pixel. The values of α are chosen according to the gradient variation at each single pixel as we discuss in Section 5.

Hereafter, we will consider gray-scale images since for colored images processing becomes a more difficult task. A naive way for colored images is to perform a similar methodology based on the one we present here but separately for each of the three layers (one per color) e.g. in the case of RGB format.

Let us start by taking a spatial discretization of the Laplacian in (4) based on a second order central difference scheme with mesh length h > 0. In such a way, ∆ transforms into a M 2 × M 2 pentadiagonal matrix ∆ h (see Fig. 1), and in the same fashion, u(t, x) is transformed into a M 2 × 1 vector-valued function u(t) which stands for the vector-arranged image pixels at time level t where M is the length and width of the image. As we commented above, the novelty of our approach consists in replacing the (only one) order derivative α of equation ( 4) with a different value of α for each single pixel of the image. This approach reads now as the linear Volterra matrix-valued equation

u(t) = u 0 + t 0 K(t -s)u(s) ds, 0 ≤ t ≤ T, ( 8
)
where u 0 is the vector-arranged initial data (noisy image), and the convolution kernel K is defined as

K(t) = I(t) • ∆ h with I(t) =                               t α 1 Γ(α 1 + 1) 0 . . . 0 0 t α 2 Γ(α 2 + 1) . . . 0 . . . . . . . . . . . . 0 . . . 0 t α M 2 Γ(α M 2 + 1)                               and, 1 < α j < 2, for j = 1, 2, . . . , M 2 .
Let us notice that K is a pentadiagonal matrix valued function, and since the Laplace transform of K exists, the well-posedness of ( 8) is then guaranteed all over the positive real line (see e.g. [START_REF] Pruss | Evolutionary Integral Equations and Applications[END_REF]).

Time discretizations

Background

Time discretizations of Volterra equations as ( 8) have been largely studied in literature; let us mention, e.g., the convolution quadrature based methods (see [START_REF] Lubich | Convolution quadrature and discretized operational calculus I[END_REF][START_REF] Lubich | Convolution quadrature and discretized operational calculus II[END_REF][START_REF] Lubich | Convolution quadrature revisited[END_REF]). In particular Runge-Kutta convolution quadrature methods (the convolution quadrature is based on classical Runge-Kutta methods) of that equations provide high order numerical methods jointly with good stability properties. In [START_REF] Lubich | Runge-Kutta methods for parabolic equations and convolution quadrature[END_REF], these discretizations have been studied in the abstract setting of sectorial operators, i.e. for convolution kernels whose Laplace transform is of sectorial type. Let us recall that a complex function G is of sectorial type if there exist 0 < θ < π/2, c ∈ R, and µ, M > 0 such that G is analytic in the sector

S θ := {λ ∈ C : | arg(λ -c)| < π -θ}, and |G(λ)| ≤ M |λ| µ , λ ∈ S θ .
Under these assumptions, the inverse Laplace transform of G can be written by means of the Bromwich formula as

g(t) = 1 2πi γ e λt G(λ) dλ,
where γ is a complex path connecting -i∞ and +i∞ parallel to the boundary of S θ with increasing imaginary part.

The convergence of these methods has been recently extended to analytic semigroups (see [START_REF] Calvo | Runge-Kutta convolution quadrature methods for well-posed equations with memory[END_REF]) where the only one requirement on the kernel is the existence of the Laplace transform (weaker assumption than the assumption of sectorial type). Since the Laplace transform of each function t α j /Γ(α j +1) in ( 8) exists, and therefore the Laplace transform of K exists as a matrix-valued function, these methods turn out to be appropriate for our purposes. In this case, if K denotes the Laplace transform of K, then the inverse Laplace transform can be written as

K(t) = 1 2πi γ e λt K(λ) dλ,
where γ(r) = a + ri, -∞ < r < +∞, for a ∈ R + . Let us notice that now the whole path γ lies on the right-hand complex plane Re(λ) ≥ a.

In order to set a suitable Runge-Kutta convolution quadrature method, we must take into account that the time regularity of the solution of ( 8) is constrained by the nature of the convolution kernel, in particular the continuity of their derivatives is guaranteed only up to the first order. Therefore, in this work, we will focus on the backward Euler convolution quadrature method, i.e., avoiding higher order schemes whose convergence will require more regularity on the solutions. This method will be sufficient to show the new features of our approach.

Convolution quadratures

For the sake of the readers convenience, we first recall the definition of the backward Euler convolution quadratures, and for the sake of the simplicity of the explanation we refer the readers to [START_REF] Calvo | Runge-Kutta convolution quadrature methods for well-posed equations with memory[END_REF][START_REF] Lubich | Convolution quadrature revisited[END_REF] for the definition of Runge-Kutta convolution quadratures in the general case and further references.

Let τ > 0 be the time step of the discretization. The convolution integral in [START_REF] Buades | The staircasing effect in neighborhood filters and its solution[END_REF]

reads t 0 K(t -s)u(s) ds = t 0 1 2πi γ e λ(t-s) K(λ) dλ u(s) ds = 1 2πi γ K(λ)Y(λ, t) dλ,
where Y(λ, t) stands for the solution of the ordinary differential equation

y ′ (t) = λy(t) + u(t), 0 ≤ t ≤ T, with y(0) = 0. ( 9
)
The backward Euler convolution quadrature is obtained as

t n 0 K(t n -s)u(s) ds ≈ 1 2πi γ K(λ)Y n (λ) dλ,
where t n = nτ, and Y n (λ) stands for the approximation of Y(λ, t n ) reached by the backward Euler method applied to [START_REF] Calvo | Runge-Kutta convolution quadrature methods for well-posed equations with memory[END_REF]. Therefore, the convolution quadrature reads

t n 0 K(t n -s)u(s) ds ≈ n j=0 Q (α) n-j u(t j ),
where the weights Q (α) j 's are defined in terms of the backward Euler characteristic polynomials quotient ρ(z)/δ(z) = z/(z -1) evaluated in the variable ξ = 1/z. In fact, such weights turn out to be the coefficients of

K 1 -ξ τ = +∞ j=0 Q (α) j ξ j .
which explicitly written read

Q (α) j = τ α                                  α 1 j 0 . . . 0 0 α 2 j . . . 0 . . . . . . . . . . . . 0 . . . 0 α M 2 j                                  • ∆ h , (10)
for j = 0, 1, 2, . . ., as M 2 × M 2 matrices. (see [START_REF] Calvo | Runge-Kutta convolution quadrature methods for well-posed equations with memory[END_REF][START_REF] Cuesta | A numerical method for an integrodifferential equation with memory in Banach spaces[END_REF][START_REF] Lubich | Convolution quadrature revisited[END_REF] for more details).

Numerical method. Convergence

Let u n be the approximation of u(t n ), for n ≥ 0. Then the time discretization of ( 8) by means of the backward Euler convolution quadrature method reads

u n = u 0 + n j=1 Q (α) n-j u j , n ≥ 1,
and keeping in mind the practical implementation, since the matrix ∆ h is not singular, the unique n-th approximation is reached by solving the linear system

I -Q (α) 0 u n = u 0 + n-1 j=1 Q (α) n-j u j , n ≥ 1. ( 11
)
In the abstract setting considered in [START_REF] Calvo | Runge-Kutta convolution quadrature methods for well-posed equations with memory[END_REF], optimal error bounds are obtained. In particular, for the backward Euler based method, the first order is reached by assuming the existence and boundedness of the second derivative of the solution. However, since the second derivative of the solution of ( 8) is merely integrable but not continuous, these results cannot be directly applied. In our case, if one takes into account the nature of the convolution kernel, then the stability proven in [START_REF] Cuesta | A numerical method for an integrodifferential equation with memory in Banach spaces[END_REF] jointly with Theorem 3.1 in [START_REF] Lubich | Convolution quadrature and discretized operational calculus I[END_REF] allow us to guarantee that the method is convergent of first order.

Besides, in [START_REF] Calvo | Runge-Kutta convolution quadrature methods for well-posed equations with memory[END_REF] an interesting result is also proven which becomes even more interesting when applications fit into the framework of image processing. For the readers convenience we recall this result in the case of the backward Euler convolution quadrature method we apply in this paper. Theorem 4.1. If u is the solution of ( 8), and u n , for n ≥ 0, is the numerical solution yielded by [START_REF] Catte | Image selective smoothing and edge detection by nonlinear diffusion[END_REF], then there exists a probability density ρ n,τ : [0, +∞) → R such that

u n = +∞ 0 u(s)ρ n,τ (s) ds, n ≥ 1. ( 12
)
The interest of this theorem is that, since ρ n,τ is a probability density, i.e., positive and satisfying +∞ 0 ρ n,τ (s) ds = 1, if u is for example positive, then the representation [START_REF] Cuesta | Some advances on image processing by means of fractional calculus[END_REF] guarantees the positivity of u n . In other words, the numerical solution [START_REF] Catte | Image selective smoothing and edge detection by nonlinear diffusion[END_REF] preserves, among other properties, the positivity of the continuous solution.

Moreover, let us point out that in [START_REF] Calvo | Runge-Kutta convolution quadrature methods for well-posed equations with memory[END_REF] it is proven that the density of probability in [START_REF] Cuesta | Some advances on image processing by means of fractional calculus[END_REF] does not depend on the equation itself but only on the numerical method. In particular, for the method considered here (see [START_REF] Cuesta | A numerical method for an integrodifferential equation with memory in Banach spaces[END_REF][START_REF] Calvo | Backward Euler method as a positivity preserving method for abstract integral equations of convolution type[END_REF]), there is an explicit expression for such a density

ρ n,τ (t) = 1 τ(n -1)! t τ n-1 e -t/τ , n ≥ 1.
Finally, we can mention some other very efficient methods to discretize (8) as for example the ones based on the discretization of the inverse Laplace transform (see [START_REF] López-Fernández | On the numerical inversion of the Laplace transform of certain holomorphic mappings[END_REF]), the ones based on the adaptive fast and oblivious convolutions (see [START_REF] López-Fernández | Adaptive, fast,and oblivious convolution in evolution equations with memory[END_REF]), or the collocation methods (see [START_REF] Brunner | Collocation Methods for Volterra Integral and Related Functional Equations[END_REF] and references therein).

Implementation

In this section we will discuss some facts concerning the implementation of (11) itself.

First of all, the choice of α's should be done according to the idea of preserving edges and corners and removing noise. Therefore, in view of the discussion in Section 2, pixels where the gradient is large should be associated with values of α close to 2. On the contrary, pixels with lower gradients should be associated with values of α close to 1. Let us notice that the practical computation of the gradient variation turns out to be very simple in a discrete setting as the one we consider in the spatial variables.

However some remarks have to be taken into account:

• On the one hand, avoiding singular situations can be yielded in both-sides values of α, i.e. α = 1 and α = 2, at least from the numerical point of view, we will not reach these values when setting α's.

In particular, in Section 6, for j = 1, 2, . . . , M 2 ,

α j ∈ [1 + ǫ, 2 -ǫ] with ǫ = 10 -3 .
• On the other hand, since extreme situations appear such as isolated noisy pixels (see Figure 3 where a gray-scale image is shown as a three-dimensional surface), a particular choice of α's is expected; in particular for those pixels, values close to 1 will be associated. On the contrary, near edges and corners (very high gradient variation) no smoothing should be required, therefore α's close to 2 will be set for that pixels.

According to this criteria, in this work, the setting of α values follows a profile distribution as in Figure 2.

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Notice that this simple and naive procedure allows one to establish different settings of distributions depending, e.g., on the characteristics of each image.

• From a computational point of view, the number of different values of α's should be limited, otherwise if one admits a number of α's as large as the number of pixels, the implementation becomes unavailable in practical cases. In fact, in this work, we set α ∈ {1 + j/N, 1 ≤ j ≤ N}, for a fixed integer N which in Section 6 turns out to be N = 100.

Let us also mention that in [START_REF] Calvo | Backward Euler method as a positivity preserving method for abstract integral equations of convolution type[END_REF] a fixed number of weights Q (α) j are computed once for all for each α j . Moreover, the practical computation of that weights has been carried out by means of the Fast Fourier Transform as in [START_REF] Calvo | Runge-Kutta convolution quadrature methods for well-posed equations with memory[END_REF], therefore saving a noticeable run-time.

Another fact of interest concerns the measure of goodness of an implementation or procedures in image filtering, restoration, and in general, to measure the quality of a processed image. Despite of in many cases a visual analysis turn out to be sufficient to determine the goodness of a methodology when applied to image processing, in this work, we consider two numerical criteria, S NR and PS NR, which have been largely used in literature (see e.g [START_REF] Dabov | Color image denoising via sparse 3d collaborative filtering with grouping constraint in luminance-chrominance space[END_REF][START_REF] Emmanuel | The curvelet transform for image denoising[END_REF][START_REF] Portilla | Image denoising via adjustment of wavelet coefficient magnitude correlation[END_REF]), and which are commonly applied to determine the quality of a processed image in the sense commented above (filtering, restoration,...). In fact, S NR and PS NR stand for the Signal to Noise Ratio and Peak Signal to Noise Ratio, respectively, and the unit for both of these ratios are decibels (dB). To be more precise, S NR of a restored image R compared to an ideal image I is defined as

S NR = 10 • log 10 var(I) var(I -R) ,
where var(x) stands for the variance of the vector x.

Here I and R are considered as vector arranged grayscale images, as in Section 3, with 256 gray levels. In the same way, PS NR is defined as

PS NR = 10 • log 10 i, j 255 2 i, j (I i, j -R i, j ) 2 ,
where I i, j , R i, j are the pixel values of I and R respectively. Notice that in restoration problems, we have a corrupted image, and try to restore the ideal image which in general is not available, but for the calculation of S NR and PS NR (from above formulas) ideal image is required. For the experiments in Section 6, we take an image (ideal image) and we perturb that image by adding up some noise. This image is then used for restoration; therefore S NR and PS NR can be computed for the restored images.

Practical results

In this section, we show the improvements in the restorations provided by our approach. To this end, we perform some experiments where a noisy image is evolved by using the Perona-Malik model ( 2) with c(s) = e -s as commented in (3) (PM), the contrast parameter (λ) is chosen such that 1.5 * σ < λ < 2 * σ as suggested in [START_REF] Gerig | Nonlinear anisotropic filtering of mri data[END_REF] (where σ is the noise variation) and the model [START_REF] Buades | The staircasing effect in neighborhood filters and its solution[END_REF] we propose (VE).

First experiment:

In the first experiment the six images shown in (Figure 4) are considered for the experimental validation of the proposed method. Each original image from the six images (Figure 4) has been perturbed by additive Gaussian noise of variance ranging from 10 to 30. In Table 1, we show the results (in terms of SNR) yielded by the restoration carried out with the mentioned procedures.

In fact, in view of Table 1, it can be observed that similar results are reached with (PM) and (VE). However we should recall that (VE) stands for a linear model whose well-posedness is guaranteed on the contrary to what happens with (PM).

For further close observation of the procedures applied to the images of Lena (Figure 4.(a)) and Elaine (Figure 4.(c)), with σ = 25 a part of such image is considered for all denoising models; in particular, a 200 × 200 part of the image has been considered (see Figure 5). The restored images by Perona-Malik model for both Lena and Elaine images, preserves edges/corners (Figure 5.(c) and Figure 5.(g)) but smooths very strongly the flat areas which causes lost of information regarding the texture of the image. Also, strong artifacts at the edges and corners has been observed as reported in [START_REF] Buades | The staircasing effect in neighborhood filters and its solution[END_REF]. The (VE) denoises the image (Figure 5.(j) and Figure 5.(h)), also by preserving edges/corners of the image; but here smoothing in flat areas is not so strong as with (PM) and it preserves the structure of the image without creating any artifacts.

In Table 2 PS NR analysis for the six images of Figure 4 with the variance of noise from 10 to 30 is given. The values of the table shows that the proposed method performs more or less in the same manner. In fact for the images from efficiency of our approach seems to be based on no more than an optical evidence. However, we will show that the efficiency of our approach is more than optical, and to this end, we will consider images where the texture plays a crucial role, and where the restoration procedures can be stressed. The images we consider in these experiments are a 250 × 250 size image of a wood (Figure 7. (a)), and a 512 ×512 size naive image (Figure 7. (e)), both strongly perturbed with gaussian noise (ratios S NR = 0.014 and PS NR = 12, see Figure 7. (b), S NR = 0.054 and PS NR = 16.2, see Figure 7. (f), respectively). In these kind of images, the texture turns out to be more important than edges preservation, and among the numerical results in Table 3, a simple overview shows the goodness of our approach vs. (PM).

However, to be more precise in our analysis, Table 3 shows a numerical evidence of the efficiency of our method; in fact it must be highlighted that S NR and PS NR are improved by using (VE) in comparison with (PM).

Conclusions and outlook

In the present work, we propose a partial differential equation based approach to image processing (filtering, denoising, enhancing,...) whose main novelty is that it fits into the framework of fractional calculus (derivatives and integrals) hence Volterra equations.

The interest of our work is twofold: On the one hand, the model we propose allows us to handle the smoothing by means of certain "viscosity" parameters which define the matrix-valued linear Volterra equation we propose. In other words, the smoothing is now handled by means of a linear partial equation, i.e. without introducing tricky nonlinear terms in the equation as many authors propose.

On the other hand, the model we propose fits into a closed mathematical setting, both from the analytical and the numerical point of view. To be more precise, the well-posedness of the Volterra equation we propose is as we said guaranteed all over the whole positive real line on the contrary what typically happens with nonlinear models, and numerical methods for its discretization have been largely studied in the literature.

As an additional interesting property of our proposal is that one can change the profile of the filter merely by changing the "viscosity" parameters setting, i.e. the distribution function we use. This allows the user to adapt easily the filter to each single image according to its characteristics.

Figure 1 :

 1 Figure 1: Sparsity pattern of the discretized Laplacian

Figure 2 :

 2 Figure 2: Profile distribution of α's.

Figure 3 :

 3 Figure 3: Three dimensional representation of (expected) isolated noisy pixels for a gray-scale image.

Figure 4 :

 4 Figure 4: Images for the first experiment: (a) Lena, (b) Boats, (c) Elaine, (d) Baboon, (e) Lady, (f) Zebra.

Figure 4 .

 4 (a) to Figure 4.(e) the values of PS NR for the purposed method are slightly larger then the (PM) model. For the image of zebra (PM) model has values of PS NR larger then (VE) model. Consider the image of Lady (Figure 4.(e)) perturbed by an additive Gaussian noise of variance σ = 25 for further investigation. Notice that the size of the image is 256 × 256 and in Figure 6, we provide original image (Figure 6.(a)), noisy image (Figure 6.(b)), restored images by (PM) and (VE) in (Figure 6.(c)-Figure 6.(d)), respectively. Figure 6.(e) and Figure 6.(f) are the residual to (PM) model and (VE) model with the noisy image, i.e, Residual = Noisy image -restored image. The residual figure 6.(e) contain more structure part of the image as compare to the figure 6.(f), which also hints that the denoising by proposed (VE) method is structure preserving. Second experiment: As in the previous discussion of S NR and PS NR for (PM) and (VE), keep very close one to each other, the

6 Fig. 7 .

 67 (h)

Figure 5 :

 5 Figure 5: Zoomed 200 × 200 part of the Lena and Elaine images: (a) Original image of Lena, (b) Noisy image perturbed by Gaussian noise (σ = 25), (c) Restoration by (PM), (d) Restoration by (VE), (e) Original image of Elaine, (f) Noisy image perturbed by Gaussian noise (σ = 25), (g) Restoration by (PM), (h) Restoration by (VE).

Figure 6 :

 6 Figure 6: Analysis of Lady image: (a) Original image of Lady, (b) Noisy image perturbed by Gaussian noise (σ = 25), (c) Restoration by (PM), (d) Restoration by (VE), (e) Residual to (PM) restoration, (f) Residual to (VE) restoration.

Table 1 :

 1 First experiment: SNR analysis

	σ	10	15	20	25	30	10	15	20	25	30
	Input SNR	13	10	7	5.6	4	13.3	9.88	7.4	5.5	4
	Method (PM)	Lena (512 × 512) 16.81 14.79 13.82 13.32 11.27	Boats (512 × 512) 15.07 13.69 12.26 11.32	9.72
	(VE)	17.31 15.00 14.46 13.73 12.78	15.74 14.10 12.65 11.57 10.32
	Input SNR	13.2	9.7	7.3	5.3	3.8	12.5	9	6.5	4.5	3
	Method (PM)	Elaine (512 × 512) 15.60 14.56 13.52 13.01 11.22	Baboon (512 × 512) 13.42 10.09 8.34	7.27	6.15
	(VE)	16.32 14.74 13.71 13.52 11.96	12.88 10.29	8.51	7.49	6.89
	Input SNR	12	8.73	6.34	4.5	3.1	13.79 10.25	7.81	5.88	4
	Method (PM)	Lady (256 × 256) 16.06 14.23 12.50 11.90	9.94	Zebra (256 × 256) 15.18 12.38 10.19	9.13	7.49
	(VE)	16.49 14.54 13.19 12.00 11.07	14.33 11.44	9.63	8.14	7.60

Table 3 :

 3 Second experiment (textured images).

Table 2 :

 2 First experiment: PSNR analysis .29 28.46 28.00 27.04 29.64 26.83 24.63 22.05 19.95 (VE) 32.14 30.78 29.35 27.85 25.76 28.00 25.45 22.57 21.45 19.00

	σ	10	15	20	25	30	10	15	20	25	30
	Input PSNR	26.66	24.5	22	20.2	18	28	24.6	22	20.27	18.7
	Method (PM)	Lena (512 × 512) 28.14 27.40 26.73 26.00 25.15	Boats (512 × 512) 29.45 28.00 27.12 25.64 24.15
	(VE)	30.00 29.33 28.37 27.08 25.94	30.49 28.42 27.25 26.32 25.02
	Input PSNR	28	24.6	22.16 20.25	18.7	28	24.6	22	20.19	18.6
	Method (PM)	Elaine (512 × 512) 29.11 28.50 26.73 25.94 23.61	Baboon (512 × 512) 28.44 24.91 22.62 22.57 21.75
	(VE)	30.98 30.19 28.41 28.00 26.08	29.02 25.92 24.07 22.87 22.34
	Input PSNR	28.24 24.88 22.49 20.69 19.23	28.24	24.7	22.27	20.3	18.78
	Method (PM)	32.06 30	Lady (256 × 256)		Zebra (256 × 256)

Acknowledgements

The research of the first author has been supported by DGI-MCI under project MTM2010-19510 (subprogram MTM) cofinanced by FEDER funds.