B Issenmann

E Falcon

Gravity wave turbulence revealed by horizontal vibrations of the container

Keywords: numbers: 47.35.-i, 05.45.-a, 47.52.+j, 47.27.-i

We experimentally study the role of the forcing on gravity-capillary wave turbulence. Previous laboratory experiments using spatially localized forcing (vibrating blades) have shown that the frequency power-law exponent of the gravity wave spectrum depends on the forcing parameters. By horizontally vibrating the whole container, we observe a spectrum exponent that does not depend on the forcing parameters for both gravity and capillary regimes. This spatially extended forcing leads to a gravity spectrum exponent in better agreement with the theory than by using a spatially localized forcing. The role of the vessel shape has been also studied. Finally, the wave spectrum is found to scale linearly with the injected power for both regimes whatever the forcing type used.

When waves of large enough amplitudes propagate within a dispersive medium, the nonlinear interactions generate waves at different scales. This energy transfer from the large scales (where the energy is injected) to the small scales (where it is dissipated) is called wave turbulence. It occurs in various domains: optical waves, surface or internal waves in oceanography, astrophysical plasma waves, Rossby waves in geophysics, elastic or spin waves in solids (for recent reviews see [1][2][3]). Since the end of the 1960s, weak turbulence theory describes the wave turbulence regimes in almost all fields involving waves [START_REF] Zakharov | Kolmogorov Spectra of Turbulence I[END_REF]. It assumes strong hypotheses such as those addressing weakly nonlinear, isotropic and homogeneous random waves in an infinite size system with scale separation between injection and dissipation of energy. It notably predicts analytical solutions for the spectrum of a weakly nonlinear wave field at equilibrium or in a stationary out-of equilibrium regime.

While homogeneity and isotropy are two premises of the theory, laboratory experiments generally use spatially localized forcing to generate wave turbulence (e.g., elastic waves on plates, or surface waves on a fluid). The use of a spatially homogeneous forcing is thus of primary interest to probe the validity domain of this theory. Previous experiments were performed by vertically vibrating a vessel filled with a fluid using the Faraday instability to homogeneously generate capillary wave turbulence [START_REF] Henry | [END_REF][6][7]. However, this forcing generates localized structures and discrete resonance peaks in the wave spectrum.

In oceans, the gravity wave spectrum depends on numerous forcing parameters (wind, fetch, sea severity, etc). Consequently, in situ spectra are usually fitted with many parameters [START_REF] Ochi | Ocean waves[END_REF], and some data are quantitatively in rough agreement with the weak turbulence predictions [START_REF] Toba | [END_REF]. However, recent well controlled laboratory experiments show deviations from the predictions for the scaling of the gravity wave spectrum when waves are locally generated (using vibrating blades at the surface of a fluid) [10,11]. In this case, the exponent of the frequency power-law spectrum of gravity waves depends on the forcing parameters (amplitude and frequency bandwidth) [10,11] instead of being independent, as expected theoretically. The origin of this discrepancy remains an open problem. It has been suggested to be due to finite size effect [10] or to the presence of strong nonlinear waves [12].

In this Letter, we study gravity-capillary wave turbulence subjected to horizontal random vibrations of the whole container. The frequency power-laws of the wave spectrum are found independent of the forcing parameters in both gravity and capillary regimes, with a rough agreement with weak turbulence theory. The probability distribution of the wave height and the scaling of the wave spectrum with the forcing amplitude is also measured. Note that horizontal vibrations of a container at a single frequency have been used to study liquid sloshing motions [13], the shape of steep capillary-gravity waves arisen through Kelvin-Helmholtz instability of two immiscible liquids [14], and in other dissipative systems driven far from equilibrium such as granular materials [15]. To our knowledge, no experiment using horizontal random vibrations of the container has been performed so far to study hydrodynamic wave turbulence.

The experimental setup is shown in Fig. 1. A circular vessel, 22 cm in diameter, is mounted on 4 ball bearing wheels and is horizontally vibrated using an electromagnetic shaker. The container is filled with water up to a depth h = 3 cm leading to an almost deep water limit (λ 2πh for our range of wavelengths λ). The shaker (LDS V406/PA 100E) is driven by a random noise forcing low-pass filtered within a frequency bandwidth between 1 Hz and f p (f p being from 5 to 7 Hz). A force sensor (FGP Instr. NTC) is fixed to the shaker axis to measure the instantaneous force F (t) applied by the shaker to the container. The instantaneous velocity V (t) of the container is measured using a home made coil placed on the shaker axis [16]. A magnet links the container to the shaker axis (and the axis of both sensors) to impose a force on the container in the direction of the shaker axis. The surface wave height η(t) is measured by a home Typeset by REVT E X made capacitive sensor [10]. This sensor moves together with the vessel, η(t) thus being measured in the container framework. Typical wave mean steepness s ranges from 0.01 to 0.10, estimated as s = k * /σ η with σ η the rms value of η(t), and k * the wavenumber of the first normal mode of the vessel (roughly corresponding to the maximum frequency of the spectrum). F (t) and V (t) are recorded for 5 min to compute the mean power P injected to the system (see below). η(t) is recorded for 5 min and 30 min, respectively, to compute its power spectrum and its probability distribution. The location of the capacitive sensor has no influence on the spectrum. We are far from conditions of resonance sloshing generating waves strongly coupled with the bulk flow such as swirling waves [13]. Note also that the maximum forcing amplitude is less than the onset of the water drop ejection or wave breaking.

A typical temporal recording of η(t) is shown in the inset of Fig. 2. η(t) displays erratic motion with η = 0. Its power density spectrum is shown in Fig. 2. Note that two peaks are visible (at 3.4 and 4.5 ± 0.2 Hz) that correspond to theoretical vessel eigenvalue modes [START_REF] Lamb | Hydrodynamics[END_REF]. Here, we are interested in the part of the spectrum not directly excited by the forcing (f > 6 Hz). At low forcing amplitude, no power law is observed and no wave turbulence regime occurs. At high enough forcing, two frequencypower laws are observed in the spectrum corresponding to the gravity and capillary wave turbulence regimes at low and high frequency, respectively. Similar results were obtained with a vibrating blade forcing [10]. The transition between both regimes occurs at a crossover frequency close to 20 Hz corresponding to λ ≈ 1 cm. The spectrum strongly decreases at high frequency ( 100 Hz) due to dissipation. When the forcing amplitude is increased, the frequency-power law fits are roughly parallel for each regime. The exponents of the frequency-power laws are shown in Fig. 3. Both gravity and capillary exponents are found independent of the forcing parameters for our range of injected power, taking values of -4.5±0.2 and -2.4 ± 0.3, respectively. These differ from results of previous studies with localized vibrating blades [10,11] where the gravity spectrum exponent was strongly dependent on the forcing parameters, taking values between -7 to -4 for the same range of injected power [10]. The exponents obtained here are however slightly different from the theory. Indeed, the gravity exponent is between -4 and -5. These values correspond respectively to the weak turbulence spectrum S grav η (f ) ∝ ǫ 1/3 gf -4 [19], and the Phillips spectrum ∝ ǫ 0 gf -5 for sharp crested waves [21], ǫ being the energy flux, f the frequency, and g the acceleration of gravity. A possible explanation for this deviation is that most of the waves are strongly nonlinear with the wave crest propagating with a preserved shape [22]. The spectrum of wave crest ridges having a fractal dimension in the range 0 ≤ D < 2 is predicted to scale as ∝ ǫ (2-D)/3 g 1+D f -3-D [11]. The gravity spectrum found experimentally in f -4.5 thus corresponds to D = 1.5. The capillary exponent is also found slightly shifted (see Fig. 3) with respect to the weak turbulence prediction [20], γ and ρ being the surface tension and the density of the fluid.

S cap η (f ) ∝ ǫ 1/2 (γ/ρ) 1/6 f -17/6
The crossover frequency f c between both regimes is measured on the spectrum in Fig. 2 as the intersection of both power laws. f c is shown in the inset of Fig. 3 for different forcing parameters. For such a horizontal forcing of the whole container, f c is roughly independent of the forcing parameters. This result differs from previous studies with a vibrating blade forcing [10] where f c depended on the forcing parameters in a range from 15 to 35 Hz for the same range of injected power and the same container size.

Consequently, horizontally vibrating the whole container is better than using vibrating blades or parametric forcing to reach a continuum wave turbulence regime independent on the forcing parameters. The main reason is that this forcing is expected to be more spatially homogeneous, and thus better approaches the corresponding theoretical hypothesis even if other assumptions are still not met, like weak nonlinearity and infinite vessel size.

To test the role of the vessel shape and of the type of forcing on gravity-capillary wave turbulence, experiments on two vessels were performed: the circular vessel (22 cm in diameter) and a rectangular one (15 × 19 cm 2 ). Two types of forcing were tested for each vessel: a localized vibrating blade and a horizontal forcing of the container as shown on Fig. 4. To avoid the predominance of the eigenvalue frequencies and sloshing modes of the rectangular container, its diagonal is chosen in the same direction as the shaker axis (see Fig. 4). We find that the frequency power-law exponent of the gravity spectrum depends on the forcing parameters: i) with the vibrating blade forcing regardless of the vessel shape, ii) with the rectangular vessel whatever the forcing type. The gravity spectrum exponent is found independent of the forcing parameters only when horizontally vibrating the circular container. Although the direction of the forcing is favored in any case, a circular vessel is more isotropic than a rectangular one due to the various wave reflection directions generated by the curved boundary. Thus, beyond homogeneity of the forcing, isotropy is also necessary to reach a gravity spectrum exponent independent of the forcing parameters.

Let us now focus on the scaling of the wave height spectrum with the mean injected power. The power injected by the shaker to the system corrected by its inertia is P(t) = (F -mdV /dt)V . m = 3.1 kg is the moving system mass (including the fluid). The mean injected power, P ≡ P , linearly increases with the variance of the shaker velocity σ 2 V ≡ V 2 (see inset of Fig. 5). • denotes the temporal average.

The height spectrum is found to scale as P 1±0.1 for both regimes over almost one order of magnitude in P (see Fig. 5). This scaling does not depend on the vessel geometry used. A similar spectrum scaling ∼ P 1 has been observed for both regimes with a vibrating blade forcing [10] for the same range of P , for the capillary regime with a parametric forcing [7], and for the inverse cascade of gravity wave turbulence [23]. This linear scaling is in disagreement with the weak turbulence theory that predicts a spectrum ∼ ǫ 1/3 in the gravity regime and ∼ ǫ 1/2 in the capillary regime (see above). Experimentally, the mean energy flux ǫ is usually estimated by the measurement of P/(ρS), S being the immersed moving surface. It is likely that a part of the power is directly provided to the bulk flow and dissipated by viscosity without cascading through the wave field. Although this mechanism is certainly present, it is unlikely to be the dominant one. Indeed, the scaling law of the spectrum with P is the same when the forcing is parametric, or by using wave makers, or by horizontally vibrating the container, while those forcings generate very different bulk flows. We rather think that the part of the injected power directly generating large scale waves only transfers a small amount of energy flux to higher harmonics compared to the direct dissipation of large scale waves by viscosity. This speculation is strengthened by recent experiments of decaying wave turbulence on the surface of a fluid that have shown that only a small part of the initial power injected into the waves feeds the capillary cascade, whereas the major part is dissipated at large scales [START_REF] Deike | [END_REF]. This unknown dissipated fraction of injected power could explain the discrepancy with weak turbulence theory for the scaling of the spectrum with P . Other possible origins of this discrepancy might be due to finite size effects (by inhibiting the energy transfers among large scale waves) [10], or the presence of strong fluctuations of the injected power [16].

Finally, the probability density function (PDF) of the wave height normalized by its rms value, η/σ η , is shown in Fig. 6. At low forcing, it is symmetric and is roughly fitted by a Gaussian function with zero mean and unit standard deviation. At high enough amplitude, it becomes asymmetric suggesting that large crests are more probable than deep troughs as usual in laboratory experiments with a vibrating blade forcing [10,24,25] or in oceanography [START_REF] Ochi | Ocean waves[END_REF]26,27]. At high forcing, the PDF tends towards a Tayfun distribution (the first non linear correction to the Gaussian) that reads p

[η] = ∞ 0 exp -x 2 -(1 -c) 2 /( 2s
2 ) /(πsc)dx where c = 1 + 2sη + x 2 , η = η/σ η , and s the mean wave steepness [25,28]. No adjustable parameter is used here. The shape of PDF(η/σ η ) thus is similar to the one obtained with a vibrating blade forcing. We also find that σ 2 η = aP for both forcing types with different proportionality constants a (see inset of Fig. 6). For a vibrating blade forcing, P was shown to be proportional to S [10]. Here, we have checked that a ∼ 1/S for both methods of forcing. Indeed, the ratio of slopes in the inset of Fig. 6 is equal (with a 4% accuracy) to the inverse of the ratio of the immersed surfaces of the blade and of the container boundary. Finally, the experimental results, P ∼ σ 2 η and S η (f ) ∼ P 1 , are coherent since by definition ∞ 0 S η (f )df = σ 2 η /(2π). In conclusion, we have introduced a new type of forcing to study gravity-capillary wave turbulence. With this spatially extended forcing, the frequency power-laws of the height spectrum are found independent of the forcing parameters for both gravity and capillary regimes. This contrasts with results of previous experiments using a spatially localized forcing where the gravity spectrum exponent depended on the forcing parameters [10,12]. Our study suggests that the dependence should be related to the inhomogeneity and the anisotropy of the localized forcing. The gravity spectrum exponent found here is slightly different from the one predicted by weak turbulence theory due to the presence of strong nonlinear waves. Finally, an explanation for the discrepancy observed with the theory for the spectrum scaling with P is also given, and applies regardless of the forcing used.

FIG. 1 .

 1 FIG. 1. Experimental setup.

FIG. 2 .FIG. 3 .

 23 FIG. 2. (Color online) Spectra of the wave height for different injected powers: P = 1.2, 14.6, 17.2, 23.6 and 28.5 mW (bottom to top). Frequency bandwidth of the forcing: 1-6 Hz (colored area). Curves are vertically shifted for clarity. Dashed (red) lines: Power-law fits of the gravity spectra. Dash-dotted (blue) lines: Power-law fits of the capillary spectra. Inset: Typical temporal evolution of the wave height.

FIG. 4 .

 4 FIG. 4. (Color online) Top view of the tested setups. The shaker horizontally vibrates a blade on the surface of water (left) or the whole container (right) as in Fig. 1. Both experiments are performed with a circular and a rectangular vessel. The light gray (red) color denotes motion parts, and black color denotes fixed parts.

FIG. 6 .

 6 FIG. 6. (Color online) Probability density function of the wave height for P = 1.2 (red), 2.7 (blue), 14.8 mW (magenta) (arrows indicate increasing power) corresponding to s =0.018, 0.037 and 0.075. Solid black line: Gaussian with zero mean and unit standard deviation. Dashed black line: Tayfun distribution with s = 0.075. Forcing: 1 -6 Hz. Inset: σ 2 η vs. P for a forcing of the whole container (•) or with a vibrating blade ( ). Slopes are respectively 8.1 and 53 cm 2 /W. Forcing: 1 -6 Hz.

We thank M. Berhanu for fruitful discussions, A. Lantheaume, C. Laroche and J. Servais for technical assistance. B. I. thanks CNRS for funding his postdoctoral research fellow. This work has been supported by ANR Turbulon 12-BS04-0005.


Julie Delon

Agn Ès Desolneux

A PATCH-BASED APPROACH FOR REMOVING IMPULSE OR MIXED GAUSSIAN-IMPULSE NOISE *

Keywords: restoration, denoising, Gaussian noise, impulse noise, non local methods, image patch, maximum likelihood estimator, order statistics

In this paper, we address the problem of the restoration of images which have been affected by impulse noise or by a mixture of Gaussian and impulse noise. We rely on a patch-based approach, which requires careful choices for both the distance between patches and for the statistical estimator of the original patch. Experiments are led in the case of pure impulse noise and in the case of a mixture. The method proves to be particularly powerful, especially for the restoration of textured regions, and compares favorably to recent restoration methods.

1. Introduction. Impulse noise in images is generally due to data loss happening during the acquisition stage (for instance in the presence of faulty pixels), storage (due to defective memory locations) or image transmission. This loss usually occurs after, or possibly at the same time as the first stages of image acquisition. During these stages, other sources of noise affect images: mainly, the photon shot noise, which is Poisson distributed, and the readout noise, which is Gaussian distributed [START_REF] Aguerrebere | Study of the digital camera acquisition process and statistical modeling of the sensor raw data[END_REF]. When no data is lost, variance stabilization techniques [START_REF] Foi | Practical Poissonian-Gaussian noise modeling and fitting for single-image raw-data[END_REF] can efficiently transform Poisson-Gaussian noise mixtures into an additive white Gaussian noise, which explains the huge literature solely dedicated to Gaussian denoising. In practice, a user may have to deal with an image that is damaged by a mixture of Poissonian, Gaussian and impulse noise. Removing this kind of mixture while preserving image details and textures is of great importance before most image analysis tasks (edge detection, segmentation, etc). In this paper, we study a simplified version of this noise model, and concentrate on the mixture of Gaussian and impulse noise.

Two models of impulse noise are generally used in the literature. In the first one, called saltand-pepper noise, each gray level is replaced with a given probability by 0 or M , where [0, M ] is the range of the original image (M = 255, in general). In this paper, we will focus on the second model of impulse noise, called random-valued impulse noise, where each gray level value is replaced with probability p, called noise ratio, by a random value in the set {0, 1, . . . , M }. Observe that detecting and removing random-valued impulse noise is much more difficult in practice than removing salt-and-pepper noise.

There is a broad literature on pure impulse noise, but the mixture of Gaussian and impulse noise is generally less studied, despite the fact that it is a more realistic noise model. The traditional approaches for impulse noise removal act locally and non linearly on images. Among them, let us mention the median and its extensions [START_REF] Pratt | Median filtering[END_REF][START_REF] Ko | Center weighted median filters and their applications to image enhancement[END_REF]. These approaches modify all pixels indifferently, while impulse noise affects only a portion of the pixels. In order to avoid this shortcoming, the trend for nearly twenty years has been to propose different impulse noise detectors and to restrict the restoration to pixels detected as corrupted. For instance, this idea underlies the switching median filter [START_REF] Sun | Detail-preserving median based filters in image processing[END_REF], the adaptive center weighted median filter (ACWMF) [START_REF] Chen | Adaptive impulse detection using center-weighted median filters[END_REF] or the pixel-wise median absolute deviation [START_REF] Crnojevic | Advanced impulse detection based on pixel-wise MAD[END_REF]. Unfortunately, median-based methods tend to destroy details and textures in images when the noise ratio is large. A successful alternative consists in combining a well chosen impulse detector, generally relying on local order statistics of gray level differences, with a global or at least semi-local restoration approach. This is the case of [START_REF] Garnett | A universal noise removal algorithm with an impulse detector[END_REF], which combines the Rank Order Absolute Differences (ROAD) for detection with a trilateral filter for restoration. Moreover, this method is used in both cases: pure impulse noise removal or mixture of Gaussian and impulse noise removal. This detection/restoration scheme is also followed by the authors of [START_REF] Chan | An iterative procedure for removing random-valued impulse noise[END_REF][START_REF] Dong | A detection statistic for random-valued impulse noise[END_REF], who rely on ACWMF, or ROLD (Rank Ordered Logarithmic Difference) for detection, before applying a variational approach to restore corrupted pixels. Other approaches for removing a mixture of Gaussian and impulse noise generally start by estimating or detecting the impulse noise and then adapt (by reducing the influence of the pixels affected by impulse noise) a Gaussian noise removal method. This is for instance the case of [START_REF] Xiao | Restoration of images corrupted by mixed Gaussian-impulse noise via l 1 -l 0 minimization[END_REF], where outliers are first detected by a median type filter and then a K-SVD dictionary learning is performed on impulse-free pixels to finally solve a l 1 -l 0 minimization problem. This is also the case of [START_REF] Cai | Fast two-phase image deblurring under impulse noise[END_REF][START_REF] Cai | Two-phase approach for deblurring images corrupted by impulse plus Gaussian noise[END_REF][START_REF] Huang | Fast Image Restoration Methods for Impulse and Gaussian Noises Removal[END_REF][START_REF] Li | Framelet Algorithms for De-Blurring Images Corrupted by Impulse Plus Gaussian Noise[END_REF] where a median type filter (ACWMF in general) is used to detect the impulse noise and then a variational framework (TV norm, or l 1 norm of framelet coefficients for instance) is used to restore the whole image. Let us mention that some of these papers (namely [START_REF] Cai | Fast two-phase image deblurring under impulse noise[END_REF][START_REF] Cai | Two-phase approach for deblurring images corrupted by impulse plus Gaussian noise[END_REF][START_REF] Li | Framelet Algorithms for De-Blurring Images Corrupted by Impulse Plus Gaussian Noise[END_REF]) tackle at the same time the problem of blur in images. The framework is thus a bit different, being at the same time more difficult (there is an additional deterioration of the image because of the blur), and an easier one (impulse noise is more efficiently detected on a blurry and therefore smoother image). Detecting the impulse noise is also the first step of [START_REF] Ghita | A new GVF-based image enhancement formulation for use in the presence of mixed noise[END_REF] where an impulse noise estimator is used to modify an anisotropic diffusion model, or of [START_REF] Yang | Mixed Gaussian and uniform impulse noise analysis using robust estimation for digital images[END_REF] where impulse noise is first removed by an adaptive cascade of median filters before removing Gaussian noise using BM3D [START_REF] Dabov | Image restoration by sparse 3d transform-domain collaborative filtering[END_REF]. The method of [START_REF] López-Rubio | Restoration of images corrupted by Gaussian and uniform impulsive noise[END_REF] takes advantage of the full probabilistic mixture model: an external method is used to predict original pixel values, followed by an EM algorithm which permits to learn the full noise model, and finally some kernel regression to get the restored image. The authors of [START_REF] Li | A new method for removing mixed noises[END_REF] use the ROAD detector of [START_REF] Garnett | A universal noise removal algorithm with an impulse detector[END_REF] to obtain a modified ("impulse controlled") distance between image patches that they use to compute the weights in the denoising scheme of Non Local means (NL-means) [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF]. The approach presented in this paper for removing mixed Gaussian-impulse noise will also be patchbased, which means that it will rely on the patch redundancy inside images. In the last fifteen years, a great deal of image processing techniques have been developed in order to take advantage of selfsimilarities of images [START_REF] Efros | Texture synthesis by non-parametric sampling[END_REF][START_REF] Awate | Unsupervised, information-theoretic, adaptive image filtering for image restoration[END_REF][START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF][START_REF] Delon | Stabilization of flicker-like effects in image sequences through local contrast correction[END_REF]. Since the introduction of the NL-means by Buades et al. [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF] to tackle Gaussian noise, the idea of relying on patch redundancy to reduce noise variance has proved to be particularly powerful. The NL-means have since been extended with success to other noise models [START_REF] Goossens | An improved non-local denoising algorithm[END_REF][START_REF] Deledalle | Iterative weighted maximum likelihood denoising with probabilistic patch-based weights[END_REF]. The idea of these restoration methods is both simple and nice, since it relies on the assumption that in a natural image, a given patch can be found almost identically in different places. The patch can then be restored in two steps : first, by finding all of its corresponding patches in the image, and second by estimating the real underlying patch behind these different damaged versions. The mathematical framework adapted to deal with this redundancy is the one of statistical estimation, as underlined in the recent contributions [START_REF] Kervrann | Local adaptivity to variable smoothness for exemplar-based image regularization and representation[END_REF][START_REF] Van De Ville | SURE-Based Non-Local Means[END_REF][START_REF] Salmon | Patch reprojections for Non Local methods[END_REF][START_REF] Duval | On the parameter choice for the Non-Local Means[END_REF]. The goal of this paper is to model properly the two steps of this estimation scheme in the case of images suffering from a mixture of impulse and Gaussian noise.

The paper is organized as follows. We start in Section 2 by giving some details about the important steps of patch-based approaches for denoising. In Section 3, we compare different measures of similarity between patches, designed to be robust to impulse and Gaussian noise. The estimation step is then tackled in Section 4, in which we also analyse the statistical properties of several possible estimators. In Section 5 we present in details our full denoising scheme. For practical reasons, this scheme is named PARIGI 1 in the paper. Experiments and comparisons with recent approaches are displayed in Section 6. For the sake of completeness, let us mention that a conference proceeding version of this work has appeared in [START_REF] Delon | A patch-based approach for random-valued impulse noise removal[END_REF].

Patch-based approaches for denoising.

2.1. Noise model. In the rest of the paper, the discrete damaged image is denoted by u and the original image is denoted by u 0 . These images are defined on a discrete domain Ω, assumed to be a bounded rectangle of Z 2 . We consider in this paper that the damaged image u is a realization of a random image U which can be written

U = (1 -T ).(u 0 + W ) + T.V, (2.1)
where T , V and W are independent random images such that • the values T (x), x ∈ Ω are i.i.d. random variables with a Bernoulli distribution of parameter p (in the sense that P[T (x) = 1] = p); p is called the noise ratio; • the values V (x), x ∈ Ω are i.i.d. random variables with a discrete uniform distribution on the range [0, M ]; • the values W (x), x ∈ Ω are independent centered Gaussian random variables with variance σ 2 : W (x) ∼ N (0, σ 2 ). In other words, the value u 0 (x) of a pixel x is first modified by adding a centered Gaussian noise and then, with probability p, replaced by a uniform value on [0, M ]. Pure impulse noise can be considered as a particular case of the previous model, obtained when σ = 0.

The fundamentals of patch-based denoising.

The idea of patch-based approaches is to take advantage of the redundancy of patches (also sometimes called the self-similarity) in images. The central hypothesis of these approaches is that for each patch P 1 in u, there exist other patches P 2 , . . . , P n in u such that all the P i are realizations of the same random patch, obtained from an underlying non-noisy patch P 0 of the original image u 0 . The first step of the whole restoration process consists in finding, for each patch P 1 , this set of corresponding replicas. In a second step, a denoised version of P 1 can be obtained by relying on any estimator P of P 0 , computed from its different damaged versions P 1 , . . . , P n .

The first step generally boils down to the choice of a similarity measure between patches. This measure should constitute a good trade-off between discriminative power and robustness to noise, as presented in Section 3. In the classical version of NL-means [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF], the similarity between patches relies on the Euclidean distance, chosen in particular for its robustness to Gaussian noise. This distance is used to compute dissimilarity weights, and these weights then balance the influence of the different patches in the final estimator P. Alternatively, the authors of [START_REF] Louchet | Total variation as a local filter[END_REF] propose to threshold the similarity measure in such a way that "a patch and one of its replicas have a probability of 0.99 to be considered as similar". The estimator P keeps only patches whose distance to P 1 is below this threshold. Other valid solutions only retain the n-nearest neighbors of P 1 [START_REF] Brox | Iterated nonlocal means for texture restoration[END_REF], or make use of clustering [START_REF] Brox | Efficient nonlocal means for denoising of textural patterns[END_REF][START_REF] Chatterjee | Clustering-based denoising with locally learned dictionaries[END_REF] or of PCA [START_REF] Deledalle | Image denoising with patch based pca: local versus global[END_REF] strategies in order to determine the set of replicas for P 1 . For the sake of simplicity, the denoising scheme PARIGI described in this paper relies on the n-nearest neighbors solution, the main difficulty to be tackled being the design of a robust similarity measure in presence of impulse noise.

The second step of the restoration amounts to a very classical estimation problem : given U 1 , . . . , U n n i.i.d. random vectors following a distribution P(θ) with θ a vector parameter, find a good estimator θ of θ. When the noise is purely Gaussian, P(θ) = N (θ, σ 2 ), with σ known, the best estimator in terms of Mean Squared Error (or quadratic risk) is the mean, which explains the usual averaging formulation of the NL-means [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF]. Recently, many authors have adopted the point of view of quadratic risk minimization as a way to optimize the parameters of the NL-means [START_REF] Kervrann | Local adaptivity to variable smoothness for exemplar-based image regularization and representation[END_REF][START_REF] Van De Ville | SURE-Based Non-Local Means[END_REF][START_REF] Duval | On the parameter choice for the Non-Local Means[END_REF] or to propose further improvements [START_REF] Salmon | Patch reprojections for Non Local methods[END_REF]. In Section 4.2, we will also make use of the quadratic risk in order to compare the efficiency of different estimators in presence of impulse and Gaussian noise.

These two steps permit to restore any patch in the noisy image. Now, since a given pixel x belongs to several patches, multiple choices are possible to obtain a restored value û(x) at x. This last step, called reprojection in the work of Salmon and Strozecki [START_REF] Salmon | Patch reprojections for Non Local methods[END_REF], corresponds to the general framework of estimators aggregation in statistics. Let us denote by P x a patch centered at x. In this paper, we assume for the sake of simplicity that patches are (2f + 1) × (2f + 1) squares (but more general shapes could be considered, as proposed in [START_REF] Deledalle | Non-local Methods with Shape-Adaptive Patches (NLM-SAP)[END_REF]), and we will denote by

Ω f = [-f, f ] × [-f, f ]
the domain of definition of the patches. A patch P x in an image u is then an image defined on Ω f by: ∀δ ∈ Ω f , P x (δ) = u(x + δ). Moreover we define the support of P x as the set x + Ω f . As a consequence, a pixel x belongs to the support of all the patches P x+δ , δ ∈ Ω f . In order to restore the value u 0 (x), we can thus rely on all the values P x+δ (-δ), δ ∈ Ω f . A common choice for the estimator of u 0 (x) is to take only into account the value P x (0) of the restored patch centered at x, as proposed in the original NL-means [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF]. Another natural choice, studied for instance in [START_REF] Salmon | Patch reprojections for Non Local methods[END_REF] is to take the mean

δ∈Ω f P x+δ (-δ) (2f + 1) 2 . ( 2.2)
The interest of this simple aggregation is that it divides the variance of the estimator of u 0 (x) by (2f + 1) 2 if the different estimations P x+δ (-δ), δ ∈ Ω f , are i.i.d. This assumption is not true in practice since the patches are overlapping, but the quality of the results is nevertheless visibly improved by this procedure, at least in the case of Gaussian noise. In this paper, we propose to merge this aggregation step with the estimation of the denoised patches. More precisely, the final estimator u(x) of u 0 (x), described in Section 5, will rely on all the values u(y) for which there exists δ ∈ Ω f such that P x+δ and P y+δ are similar.

3. Robust distance between patches. This section focuses on the search of similar patches in an image affected by a mixture of impulse and Gaussian noise. As underlined before, the success of any patch-based denoising procedure depends greatly on the ability to find the replicas of a patch in the noisy image u. More precisely, for a given patch P in u, we aim at discovering all patches Q such that the unknown original patches P 0 and Q 0 in u 0 are equal or at least similar. In presence of impulse noise, the Euclidean ( 2 ) distance between P and Q contains outliers and cannot be trusted. In the following, we compare several alternative measures designed to be robust to this kind of noise.

Generalized likelihood ratio.

A first possibility, suggested by the recent work [START_REF] Deledalle | How to Compare Noisy Patches? Patch Similarity Beyond Gaussian Noise[END_REF], is to make use of the generalized maximum likelihood ratio

GLR(P, Q) = sup T P[P, Q|P 0 = Q 0 = T] sup T P[P|P 0 = T] sup T P[Q|Q 0 = T] . (3.1)
In the case of pure random-valued impulse noise, this ratio becomes

GLR(P, Q) = 1 + (M + 1)(1 -p) p 2(n sim -|P|) , (3.2)
where n sim is the number of pixels in Ω f that have exactly the same grey level value in P and Q, and |P| = (2f + 1)2 is the size of the definition domain of P. Obviously, this measure is too rigid to be used in practice. Indeed, two patches can be very similar and at the same time be such that n sim = 0.

In the mixed Gaussian-impulse case, the ratio becomes (with f θ the Gaussian distribution with mean θ and variance σ 2 ):

GLR(P, Q) = sup T x∈Ω f p M + 1 + (1 -p)f T(x) (P(x)) p M + 1 + (1 -p)f T(x) (Q(x)) x∈Ω f p M + 1 + (1 -p)f P(x) (P(x)) p M + 1 + (1 -p)f Q(x) (Q(x)) . (3.3)
The denominator of this formula can be simplified into the constant ( p M +1 + 1-p σ √ 2π ) 2|P| , independent of P and Q. The numerator is more complex. It consists of a product of different terms num(x)

when x ∈ Ω f : num(x) equals ( p M +1 + 1-p σ √
2π e -(P(x)-Q(x)) 2 /8σ 2 ) 2 when (P(x) -Q(x)) 2 /σ 2 is small (because its value is a sup that is achieved for T(x) = (P(x)+Q(x))/2) and num(x) approximately equals

( p M +1 + 1-p σ √ 2π )( p M +1 + 1-p σ √ 2π e -(P(x)-Q(x)) 2 /2σ 2 ) when (P(x)-Q(x)) 2 /σ 2 is large (it is achieved for T(x) P(x) or T(x) Q(x)
). As a consequence, the log of GLR(P, Q) is a sum over x of a function of (P(x) -Q(x)) 2 , the function depending on whether (P(x) -Q(x)) 2 /2σ 2 is small or not. The main difficulty when using this generalized likelihood ratio (GLR) is that we do not know σ in practice, but its main advantage is that it automatically adapts to large or small values of (P(x) -Q(x)) 2 . It is thus related to robust similarity measures, that we explain in the following section.

Weighted Euclidean distances.

A second possibility is to rely on a robust similarity measure, inspired by order statistics [START_REF] David | Order statistics. Wiley series in probability and mathematical statistics[END_REF]. It is known that, for non-noisy patches, the 2 distance is a good distance to measure the similarity between the patches. But when impulse noise is present, many large values of (P -Q) 2 i , 1 ≤ i ≤ n, are observed that are due to the noise and the 2 distance then contains many outliers. Our distance is built to be similar to the Euclidean distance but robust to large values created by the noise. Let |P -

Q| (1) ≤ |P -Q| (2) ≤ • • • ≤ |P -Q| (2f +1) 2 be
the values obtained by ordering the (2f + 1) 2 values of the differences |P(y) -Q(y)|. We propose to rely on a distance of the form

D(P, Q) = (2f +1) 2 k=1 w k |P -Q| 2 (k) , (3.4)
where w 1 , . . . , w n are positive weights. There is a wide choice for the weights w k . For instance, if they are all equal, we find again the 2 distance. If we set w k = 1 for k below a number T and w k equal to 0 above, we obtain a trimmed sum of order statistics : the distance relies only on the smallest distances between gray levels in the support of the patch. We denote this distance D trimmed and we choose for T the value T = sup{k; B((2f + 1) 2 , k, (1 -p) 2 ) > 0.99}, where B denotes the tail of the binomial distribution 2 . This choice comes from the following property: if P and Q are two independent random patches such that P 0 = Q 0 , the probability that the k th difference |P -Q| (k) stems from two untouched pixels is B((2f + 1) 2 , k, (1 -p) 2 ) (with the approximation that the smallest distances correspond to untouched pixels).

In cases of low signal-to-noise ratio, this trimmed distance becomes difficult to use, since the number T of pixels taken into account rapidly tends toward zero. In order to avoid this problem, we propose to keep all the weight w k positive and equal to the probability that the value |P -Q| (k) stems from untouched pixels. We denote this distance D weighted :

D weighted (P, Q) = (2f +1) 2 k=1 B((2f + 1) 2 , k, (1 -p) 2 )|P -Q| 2 (k) . (3.5)
Notice that when there is no impulse noise (p = 0), then the two distances D weighted and D trimmed are both the usual 2 distance between patches.

The choice of the distance D weighted is reinforced by modeling the problem the following simple way: assume the two patches P 0 and Q 0 are both constant patches, and denote by a and b their respective constant gray level value. Then, when the two patches are damaged independently by some impulse noise of intensity p, the random variables 2 and dF is a probability distribution on R + (that we don't need to compute). In this framework, we have the following theorem that shows that choosing w k = B((2f +1) 2 , k, (1-p) 2 ) for the weighted distance in (3.4) is, in some sense, optimal to discriminate the two constant patches. Theorem 3.1. Let dF be a probability distribution on R + and let

(P -Q) 2 i , 1 ≤ i ≤ n, are i.i.d. following a distribution of the form (1-p) 2 δ (b-a) 2 +(1-(1-p) 2 )dF , where δ (b-a) 2 denotes the Dirac distribution at (b -a)
X (α) 1 , . . . , X (α) n be n i.i.d.
random variables, distributed according to the mixture distribution qδ α + (1 -q)dF , where q ∈ [0, 1], and δ α denotes the Dirac distribution at α ∈ [0, +∞). Let w = (w 1 , . . . , w n ) be a vector of weights in [0, +∞), and let

D w (α) = E n k=1 w k X (α) (k) , where X (α) (k) denotes the k-th order statistic of X (α) 1 , . . . , X (α)
n . Assume that the probability distribution dF satisfies p α := P(Y ≤ α) → 0 when α → 0 where Y is a random variable distributed following the law dF . Then when α is small, we have

D w (α) = D w (0) + α n k=1 w k B(n, k, q) + o(α),
and notice that this first-order Taylor expansion is independent of the probability distribution dF . Now, a direct consequence of this theorem is that among all the distances such that w 2 k is fixed, the one that is the most able to discriminate, on the average, two constant equal patches (this corresponds to α = 0) from two constant but non-equal patches (α > 0), is the one such that n k=1 w k B(n, k, q) is as large as possible. This is achieved when w k is proportional to B(n, k, q) where q = (1 -p) 2 . We recognize the weighted distance defined in (3.5).

Proof. Let us first introduce some notations. We will denote m + k,l (α) and m - k,l (α) the following expected values:

m + k,l (α) = E(Y (k) | Y 1 ≥ α, . . . , Y l ≥ α) and m - k,l (α) = E(Y (k) | Y 1 ≤ α, . . . , Y l ≤ α),
where Y 1 , . . . , Y l are i.i.d. random variables following the law dF . Let also N α be the random variables that counts the number of X (α)

i , 1 ≤ i ≤ n, that are equal to α. Then N α follows the binomial distribution of parameters n and q, and we will denote, for 0

≤ k ≤ n, b k = P(N α = k) = b(n, k, q) = n k q k (1 -q) n-k . Now for α ≥ 0 we can write D w (α) = n k 0 =0 n k=1 w k E X (α) (k) | N α = k 0 P(N α = k 0 ) = n k 0 =0 b k 0 n-k 0 l=0 b(n -k 0 , l, p α ) l k=1 w k m - k,l (α) + (w l+1 + . . . + w l+k 0 )α + n k=l+k 0 +1 w k m + k-l-k 0 ,n-l-k 0 (α)   .
But notice that we also have, for any α ≥ 0,

D w (0) = n k 0 =0 b k 0 n-k 0 l=0 b(n -k 0 , l, p α )   l k=1 w k+k 0 m - k,l (α) + n k=l+k 0 +1 w k m + k-l-k 0 ,n-l-k 0 (α)   .
Consequently, we get

D w (α) = D w (0) + n k 0 =0 b k 0 n-k 0 l=0 b(n -k 0 , l, p α ) l k=1 (w k -w k+k 0 )m - k,l (α) + (w l+1 + . . . + w l+k 0 )α .
Finally, since 0 ≤ m - k,l (α) ≤ α, and p α goes to 0 when α goes to 0, we have that the first-order Taylor expansion is given by the term l = 0 in the sum above. More precisely, we get

D w (α) = D w (0) + n k 0 =0 b k 0 (w 1 + . . . + w k 0 )α + o(α) = D w (0) + α n k=1 w k B(n, k, q) + o(α),
which is the announced result.

Comparison of the different distances.

When an image is not affected by noise, the 2 distance is a good way to measure the similarity between patches. A nice way to compare different robust distances in presence of noise is to measure their ability to find similar patches at the same locations as those found by the 2 distance in the non-noisy image. A robust distance should be able to preserve as much as possible the ordering of the patches provided by the 2 distance in the non-noisy image. The comparison procedure works as follows. Given an original patch P 0 in an original image u 0 , and a neighborhood V (that is a square window of half-size t centered at the center of the patch), we will denote by P 0

x the patch of u 0 centered at x ∈ V . Let the previous patches in the noisy image u be respectively denoted by P and P x . For a distance D between patches, and for n ≥ 1 integer, let x 1 , . . . , x n denote the centers of the patches that achieve the n smallest values of the distance D to the noisy patch P: 0 = D(P, P x 1 ) ≤ D(P, P x 2 ) ≤ . . . D(P, P xn ). From left to right: the original patch (marked by a white window of size 7 × 7 pixels and its neighborhood V ), a noisy version of the image (p = 0.2 and σ = 10), and the mean number of selected patches for the different distances when σ = 10, p = 0.2 (third column) and p = 0.5 (fourth column). The 2 distance corresponds to blue curves, the trimmed distance D trimmed to green curves, the weighted distance D weighted to red curves and the "impulse-controlled" (IC) distance used by [START_REF] Li | A new method for removing mixed noises[END_REF] to indigo curves. See the text for some comments on these graphics. Now, we count how many of these patches would have been also chosen among the n-nearest patches for the 2 distance in the original image. More precisely, if we denote by x 0 1 , . . . , x 0 n the centers of the patches that achieve the n smallest 2 distances to the patch P 0 in u 0 , we then compute

N D (n) = #{1 ≤ i ≤ n such that x i ∈ {x 0 1 , . . . , x 0 n }}.
We make n vary from 1 to 40 and we plot the number N D (n) of "well-selected" patches as a function of n. We ran this experiment for four different distances between patches: the 2 distance (blue curves), the trimmed distance D trimmed (green curves), the weighted distance D weighted defined in Equation (3.5) (red curves) and the "impulse-controlled" (IC) distance used by [START_REF] Li | A new method for removing mixed noises[END_REF] (indigo curves). This last distance is also a weighted 2 distance but where the weights are now pixel-dependent and related to the noise detector ROAD [START_REF] Garnett | A universal noise removal algorithm with an impulse detector[END_REF] (see Section 5.1 for a more detailed description of ROAD).

The results are shown on Figures 3.1 and 3.2, for different patches: an edge patch (its boundary is marked by the white window in the top left image), a T-junction patch, and three different texture patches; and for different noise parameters: σ = 10 and p = 0.2; or σ = 10 and p = 0.5. The obtained curves are averaged from 200 noise samples for each set of parameters. The black diagonal (that recalls the value of n) is the value that N D (n) should try to reach in order to be a "good" distance between the noisy patches. As a consequence, the higher the curve n → N D (n) is, the better the distance is.

From these experiments we can make the following remarks: -The 2 distance (blue curves) between the noisy patches always achieve the worst performances.

-The trimmed distance (green curves) becomes unable to correctly select the patches when the impulse noise is too important (the parameter T used to trimmed the distance is too small -for Fig. 3.2. Comparison of different distances on noisy texture patches. From top to bottom: three different texture patches (the scarf of Barbara, the tablecloth in Barbara, the fur of the baboon). From left to right: the original patch (marked by a white window of size 7 × 7 pixels and its neighborhood V ), a noisy version of the image (p = 0.2 and σ = 10), and the mean number of selected patches for the different distances when σ = 10, p = 0.2 (third column) and p = 0.5 (fourth column). The 2 distance corresponds to blue curves, the trimmed distance D trimmed to green curves, the weighted distance D weighted to red curves and the "impulse-controlled" (IC) distance used by [START_REF] Li | A new method for removing mixed noises[END_REF] to indigo curves. See the text for some comments on these graphics.

instance its value is only T = 6 when p = 0.5 for patches of size 7 × 7 pixels).

-The IC distance (indigo curves) and the weighted distance (red curves) have similar performances on the edge patch and on the second and third texture patches.

-On the T-junction patch, the IC distance performs better, while on the first texture patch (a part of the scarf of Barbara), the weighted distance seems to be a more relevant choice.

-On the third texture, all distances have bad performances (and even very bad ones when p = 0.5). This can be explained by the fact that this texture looks like "noise", and in that case, the nearest patches selected by the distance are quite randomly chosen.

Our conclusion of these experiments is the following: since the only case where the IC distance performs clearly better than the weighted distance is on the T-junction patch and since such a patch is rather usually a "rare" patch in images, we decide in all the following to use the weighted distance between noisy patches. Indeed, the IC distance of [START_REF] Li | A new method for removing mixed noises[END_REF] is rather empirically defined, and relies on several parameters whose influence is not easy to interpret in practice. Conversely, the weighted distance has a simple definition and only relies on the value of p. We will see in Section 5.1 how the noise ratio p can be efficiently estimated.

4. Choice of the estimator P. In this section, we aim at defining a good estimator P of the underlying patch P 0 behind different damaged versions P 1 ,... P n . For the sake of simplicity, we assume that these patches are independent realizations of the same random patch, following the noise model (2.1).

Let δ ∈ Ω f be a fixed point in the domain of definition of the patches (for instance δ = 0 corresponds to the center of the patches) and let us denote by X 1 , . . . , X n the random variables corresponding to the realizations P 1 (δ), . . . , P n (δ). The X i are then real-valued random variables that are i.i.d. and follow the mixture distribution

X i ∼ (1 -p) N (µ, σ 2 ) + p U [0,M ] , (4.1)
where U [0,M ] denotes the uniform law on the discrete set [0, . . . , M ] and N (µ, σ 2 ) denotes the Gaussian distribution of mean µ and variance σ 2 . We wish to estimate µ (and possibly σ). Let us start with the simpler problem of pure impulse noise : in this case, the Gaussian part of the mixture degenerates into a Dirac distribution at µ.

4.1. Pure impulse noise. In this pure impulse noise case, the mean and variance of the X i 's can be written

E[X i ] = (1 -p)µ + p M 2 and (4.2) Var[X i ] = µp(1 -p)(µ -M ) + p 2M 2 + M 6 -p 2 M 2 2 . (4.3)
In the following, we study the relevance of different estimators of µ.

4.1.1. Mean. The mean X n = 1 n n i=1 X i is known to be the maximum likelihood estimator (MLE) of µ for additive Gaussian noise. For pure impulse noise, the mean of X n equals

E[X n ] = E[X i ] = µ + p M 2 -µ , (4.4)
We can thus derive a first unbiased estimator of µ, as

X n = 1 1-p X n -p M 2 .
The quadratic risk (i.e the mean squared error or MSE) of this unbiased estimator equals

R[ X n ] = E[| X n -µ| 2 ] = Var[ X n ] = 1 1 -p 2 Var[X i ] n . (4.5)
As a consequence, this risk makes the estimator useless for practical purposes. Indeed, for example, when M = 255, p = 0.5 and µ = 100, then Var[ X n ] 108 √ n , which means that we need more than 468 samples to make the square root of the risk go below 5 gray levels.

4.1.2. Median. The spatial median filter is one of the most famous filters used to remove impulse noise. Let us study its properties as an estimator of µ. Assume, for the sake of simplicity, that n is odd: the median of X 1 , . . . X n is X ( n+1 2 ) , where X (k) is the k th order statistic (or k th smallest value) of X 1 , . . . X n . The distribution function of the median can thus be easily computed as:

P[X ( n+1 2 ) ≤ k] = B n, n + 1 2 , F (k) := n j≥ n+1 2 n j F (k) j (1 -F (k)) n-j , (4.6)
where

F (k) := P[X i ≤ k] = p k+1 M +1 + (1 -p)1 k≥µ is the distribution function of the X i 's. After some computations, it follows that E[X ( n+1 2 ) ] = M - M -1 k=0 B n, n + 1 2 , F (k) . (4.7)
Numerically, this bias is negligible as long as p ≤ 0.4. Nonetheless, it increases dangerously when p becomes larger than 0.5. The variance of the estimator can be computed in the same way. At the end, the quadratic risk of the estimator equals: This huge difference in behavior is due to the bias of the median estimator which attains very large values as soon as p ≥ 0.5 and for gray levels µ far from the middle of [0, M ]. As a conclusion, the median estimator is clearly more interesting than the one built upon the mean, but in practice, its quadratic risk is not controlled, except for small values of p. where the histogram attains its maximum. Indeed,

R[X ( n+1 2 ) ] = (M -µ) 2 - M -1 k=0 (2(k -µ) + 1)B n, n + 1 2 , F (k) . ( 4
X n = arg max µ log P[X 1 , . . . , X n |µ ] = arg max µ n i=1 log P[X i |µ ] = arg max µ n i=1 log (1 -p)δ X i =µ + p M + 1 = arg max (h).
If the histogram has several maxima, one of them is chosen randomly, with equal probabilities.

Consequently, if we denote by X n the random variable corresponding to this MLE, the law of X n is given by

X n ∼ (1 -q n,p ) δ µ + q n,p U [0,M ] , (4.9)
where q n,p ∈ [0, 1] is a probability that depends on p, n (and M ). The probability for the estimator to estimate correctly µ is then P[ X n = µ] = 1 -q n,p + qn,p M +1 . The mean and the variance of X n are easily obtained by replacing p by q n,p in Equations (4.3) and (4.2). It follows that the quadratic risk of X n as an estimator of µ can be written

R[ X n ] = E[| X n -µ| 2 ] = (Bias[ X n ]) 2 + Var( X n ) = q n,p µ(µ -M ) + 2M 2 + M 6 . (4.10)
and is proportional to q n,p . In order to keep R[ X n ] smaller than 5 gray levels for M = 255 and any value of µ in [0, M ], the probability q n,p must be smaller than 10 -3 . Now, it can be shown that the value of q n,p is quickly decreasing when n increases. Indeed, let K := #{i, X i = µ} be the random variable corresponding to the number of samples equal to µ. Among the n -K remaining samples with values different from µ, let C K be the number of K-uples with equal values. Observe that the event X n = µ is realized when C K = 0. Of course, C K = 0 whenever K ≥ n+1 2 . As a consequence, if we assume for sake of simplicity that n is odd, we get 1 -q n,p + q n,p

M + 1 = P[ X n = µ] ≥ P[K ≥ n + 1 2 ] + P[2 ≤ K ≤ n -1 2 and C K = 0] = n k= n+1 2 n k p n-k (1 -p) k + n-1 2 k=2 n k p n-k (1 -p) k (1 -α k ), where p = P[X 1 = µ] = p -p (M +1
) and where

α k := P[C k ≥ 1] ≤ E[C k ] = n-k k 1 M k-1 . It follows that 1 -q n,p + q n,p M + 1 = P[ X n = µ] ≥ B(n, n + 1 2 , 1 -p) + n-1 2 k=2 n k p n-k (1 -p) k 1 - n -k k 1 M k-1 = 1 -p n -n(1 -p) p n-1 - n-1 2 k=2 n k n -k k 1 M k-1 (1 -p) k p n-k .
For reasonable values of p and n (p ≤ 0.8 and n ≤ 50), we checked empirically that the inequality in the above equation is in fact almost an equality, and that the last sum can be approximated by its first term (k = 2), which is broadly dominant. This allows us to estimate the parameter q n,p of the distribution (4.9) very easily in function of p. The left part of Figure 4.2 shows the values of q n,p as a function of p and n, when M = 255. The white line on the figure corresponds to the level line q n,p = 10 -3 . This line yields the number n of samples necessary to keep the square root of the risk always below 5 gray levels (whatever the value of µ), as a function of p. This number of samples remains quite reasonable in comparison to previous estimators. Moreover, the risk is controlled for all values of µ, unlike the median estimator. The risk of X n as a function of p and n is illustrated on the right part of Figure 4.2, for µ = 100.

Mixture case.

The previous study tends to prove that in the case of pure impulse noise, the Maximum Likelihood Estimator is preferable to the median estimator for denoising.

In the full mixture case, we would like to estimate both the mean µ and the variance σ 2 of the Gaussian part. Observe that in practice, this full mixture case is a realistic model for the estimation problem even in presence of pure impulse noise. Indeed, similar patches in an image are never exactly equal and it is a sound hypothesis to assume that their inner variability can be modeled with a standard deviation σ (that is pixel-dependent).

The expectation of the X i 's is unchanged (see Equation (4.2)) and their variance becomes

Var[X i ] = µp(1 -p)(µ -M ) + p 2M 2 + M 6 -p 2 M 2 2 + (1 -p)σ 2 . (4.11)
The only difference with the pure impulse case consists in the additional term +(1 -p)σ 2 .

Mean and Median.

The properties of the mean and median as estimators of µ are quite similar to the pure impulse case. The variance of the unbiased estimator X n is still obtained from Equation (4.5). The law of the median is easily computed by replacing the distribution function

F in Equation (4.6) by F (k) = p k+1 M +1 + (1 -p)G µ,σ (k)
, where G µ,σ is the distribution function of the N (µ, σ 2 ) distribution. Formulas (4.7) and (4.8) for the bias and risk of this estimator are generalized in the same way. Again, the risk of this estimator is uncontrolled as soon as p increases.

Maximum likelihood.

In the mixture case, the Maximum Likelihood Estimator (MLE) of (µ, σ), that we will denote by ( X n , σ), is defined as

( X n , σ) = arg max µ ,σ log P[X 1 , . . . , X n |µ , σ ] = arg max µ ,σ n i=1 log p M + 1 + (1 -p)g µ ,σ (X i ) ,
where g µ ,σ is the Gaussian probability density of mean µ and variance σ 2 . This can be rewritten by using the empirical distribution h of the values {X 1 , . . . , X n } on [0, M ], which yields

( X n , σ) = arg max µ ,σ h * f σ (µ )
where

f σ : m → log p M +1 + 1-p σ √ 2π e -m 2 σ
. Observe that the MLE of pure impulse noise can be recovered from this formula by taking σ = 0. In the same way, the MLE in the case of pure Gaussian noise is obtained by taking p = 0 in the previous formula.

The study of the bias and variance of this maximum likelihood estimator is far more complex than in the pure impulse case. In particular, there is no obvious close formula for the distribution of X n . In order to evaluate the quality of this estimator, we replace this study by an empirical estimation of the quadratic risk R[

X n ] = E[( X n -µ) 2 ]
as a function of p and n for different values of σ and for µ = 100 (as observed in the pure impulse case, this MSE does not depend too much on the value of µ). For each σ and for each value of p, we draw a white cross at the minimum number n that should be chosen in order to ensure that the square root of the risk or MSE is smaller than 5. This number obviously increases with σ and p. It remains nevertheless reasonable for p ≤ 0.5 and σ ≤ 15.

Choosing n in practice. The idea behind the previous empirical computation is to evaluate the appropriate number n of patches necessary to compute X n with a controlled risk. Assuming that we have an accurate estimation of p and σ, this number of trusted patches can be computed empirically and used to estimate µ in a reliable way. Now, recall that X 1 , . . . , X n correspond in practice to realizations P 1 (δ), . . . , P n (δ) for a set of similar patches P 1 , . . . , P n . In natural images, even the most similar patches are seldom equal, and present some differences. This situation has two consequences. First, the value of σ in the model (4.1) should account both for the Gaussian part of the noise and for the inner variability between patches. Even in the case of pure impulse noise, choosing σ > 0 in order to estimate the right number of trusted patches n might improve the results. Second, the number of trusted patches n should always result from a compromise between the control of the theoretical risk and the redundancy of the image to be denoised. Indeed, the inner variability between patches considered as similar might increase with n, depending on the image regularity.

In practice, we observed that choosing σ = 10 to estimate n in the case of images suffering from pure impulse noise usually yields a good compromise between image redundancy, patch variability and risk control. Once p is estimated, the values of n used in all our experiments on pure impulse noise (see Section 6) correspond to the white crosses of Figure 4.3 (b). In the case of mixture noise, the situation is more complex, since the values of σ are unknown in practice. In order to reduce the number of parameters, we resolved to choose the values of n obtained for σ = 15 in all our experiments involving a mixed noise (white crosses of Figure 4.3 (c)). For the sake of completeness, 5. Denoising scheme. We are now in a position to fully describe PARIGI, a denoising scheme built upon the maximum likelihood estimator presented in Section 4.2.2 and the weighted distance between patches introduced in Section 3. This denoising scheme requires to estimate the noise ratio p 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 n 13 17 20 27 37 56 86 101 p beforehand. We describe in Section 5.1 different approaches devoted to estimate p globally on the damaged image. Section 5.2 then details the successive steps of our patch-based denoising scheme.

5.1. Estimation of p. Different methods are possible to estimate p globally on the damaged image. In this paper, we propose to rely on impulse noise detectors.

Many impulse noise detectors are proposed in the literature, such as the detectors ROAD [START_REF] Garnett | A universal noise removal algorithm with an impulse detector[END_REF], ACWMF [START_REF] Chen | Adaptive impulse detection using center-weighted median filters[END_REF], ROLD [START_REF] Dong | A detection statistic for random-valued impulse noise[END_REF] or the pixel-wise MAD [START_REF] Crnojevic | Advanced impulse detection based on pixel-wise MAD[END_REF], to cite only a few. The goal of these schemes is to yield a map of noisy (impulse affected) pixels in an image damaged by impulse noise, or by a mixture of impulse and Gaussian noise. An estimation of p can naturally be derived from this estimation by computing the ratio between the number of noisy pixels and the image size. For the sake of simplicity, we restrict ourselves in this paper to the detectors ROAD and ACWMF, which are quite complementary.

The detector ROAD (for "Rank Ordered Absolute Differences"), proposed in [START_REF] Garnett | A universal noise removal algorithm with an impulse detector[END_REF], can be described as follows: for each pixel x, the absolute differences between u(x) and u(y) are computed for all y = x in a centered 3×3 patch around x. These differences are ordered. The value ROAD(x) is obtained by computing the sum of the 4 smallest differences. This value measures how close u(x) is from its neighbors. When ROAD(x) is above a given threshold τ , set as 70 in our experiments, x is considered as noisy.

The detector ACWMF [START_REF] Chen | Adaptive impulse detection using center-weighted median filters[END_REF] works as follows. For a given pixel x, a weighted median of order k is defined as

med k (x) = median   {u(y), y ∈ V x } ∪ {u(x), . . . , u(x) 2k times }   , (5.1)
where V x is a square window centered at x. This weighted median boils down to the usual median filter when k = 0 and to the identity when k is large enough. Now, let d k = | med k (x) -u(x)| and let t k be a decreasing sequence of well chosen thresholds. The pixel x is presumed to be noisy if there exists one k such that d k is above the threshold t k . In practice, the authors of [START_REF] Chen | Adaptive impulse detection using center-weighted median filters[END_REF] recommend to use 3 × 3 patches and to compute the four thresholds t k , k = 1, . . . , 4 as

t k = s. median {|u(y) -med 0 (x)|, y ∈ V x } + δ k , (5.2)
with [δ 0 , δ 1 , δ 2 , δ 3 ] = [START_REF] Wang | Image quality assessment: From error visibility to structural similarity[END_REF][START_REF] Garnett | A universal noise removal algorithm with an impulse detector[END_REF][START_REF] Chatterjee | Is denoising dead?[END_REF][START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF]] and 0 ≤ s ≤ 0.6. Figures 5.1 (a,b) show the quality of both estimators for different values of p on the 512 × 512 image Lena in the case of pure impulse noise (on the left) and in the case of a mixture of impulse noise and Gaussian noise with σ = 10 (on the right). For each value of p between 0.1 and 0.9 with a step of 0.1, we show the boxplot of the estimation of p for 100 different noise samples. The ROAD-based estimation of p, shown on Figure 5.1 (a), is quite accurate in both cases and fast to compute. As already underlined in [START_REF] Chen | Adaptive impulse detection using center-weighted median filters[END_REF], the detector ACWMF is efficient when the noise ratio is not too high, typically p ≤ 0.25. Above this value, the number of pixels considered as touched by noise is highly underestimated and the estimation of p is not usable in practice. As a consequence, we decided in all our experiments to use the estimation of p based on the detector ROAD.

Implementation details.

The estimation of p is the first step of PARIGI. The second step consists in estimating n, the number of nearest neighbors taken into account in the space of patches to compute the estimator P . Assuming that the estimation of p is accurate, this number of trusted patches is computed empirically as explained in Section 4.2.2.

The algorithm continues as follows. For each point x in Ω, we seek the n nearest neighbors P y of P x for the distance D weighted introduced in Section 3. We restrict this investigation to the points y spanning a (2t + 1) × (2t + 1) square V x centered at x (see Figure 5.2)

V x = {x + δ; δ ∈ [-t, t] × [-t, t]}. (5.3)
This permits to define V n x , the subset of V x corresponding to the n nearest neighbors of P x . The MLE ( u(x), σ(x)) is then computed at each x from the n-tuple (u(y), y ∈ V n x ). The image u constitutes a first denoised version of u. In practice, this restored image is sometimes a little bit too smooth. At the same time, some impulse pixels can remain. In order to recover the grain of the original image and to eliminate these last impulses, we take into account the estimated standard deviation σ(x) at each point x for defining a map of noisy pixels:

M = {x ∈ Ω; | u(x) -u(x)| > σ(x)}.
We then compute the mixture of u and u: u 2 (x) = u(x).1 x∈M + u(x).(1 -1 x∈M ), and eventually apply the previous denoising steps to u 2 in order to obtain u 2 . These steps can be repeated a few times to improve the final result. The interest of these iterations is shown on the last line of Figure 6.7. In the rest of the experiments, we repeat this step twice and the output of the algorithm is given by u 3 . In all our experiments, the half-size of the research neighborhood V x is set to t = 7. The half-size of the patches is set to f = 3 for pure impulse noise and to f = 8 in the case of a mixture of Gaussian and impulse noise.

A refined version of the algorithm can be obtained by following the approach introduced in [START_REF] Salmon | Patch reprojections for Non Local methods[END_REF]. Notice that a point x belongs to all patches P x+δ , δ ∈ Ω f . The idea of the refinement is to take into account in the estimation all the information from these patches. In this version, the MLE is computed at each x from the set of values

u(y -δ); δ ∈ Ω f and y ∈ V n x+δ .
This refined version is the one used in the experimental section. The whole refined algorithm is described in Algorithm 1.

6. Experiments and discussion. This section is devoted to the experimental analysis of the denoising scheme introduced in the previous section. We confront PARIGI with the recent state of the art approaches [START_REF] Garnett | A universal noise removal algorithm with an impulse detector[END_REF][START_REF] Dong | A detection statistic for random-valued impulse noise[END_REF] on pure impulse noise and with [START_REF] Garnett | A universal noise removal algorithm with an impulse detector[END_REF][START_REF] Xiao | Restoration of images corrupted by mixed Gaussian-impulse noise via l 1 -l 0 minimization[END_REF] on the mixture of impulse and Gaussian noise. Comparison results are provided both under the form of PSNR tables and of visual experiments. In the case of pure impulse noise, we also show MAE (Maximum Absolute Error) and SSIM (Structural SIMilarity index) tables. The last part of the section is devoted to a short discussion on the link between image regularity models and denoising approaches.

6.1. Pure impulse noise. In this paragraph, we present some experiments on pure impulse noise. Tables 6.1, 6.2 and 6.3 give the PSNR, MAE and SSIM obtained with the three methods ( [START_REF] Garnett | A universal noise removal algorithm with an impulse detector[END_REF][START_REF] Dong | A detection statistic for random-valued impulse noise[END_REF] and our approach) on the 512 × 512 classical images Lena, Bridge, Baboon, Barbara, Cameraman, Boat, Peppers and Goldhill. Let us mention at this point that most of these test images are available in more than one version on the internet, the main differences concerning their

= u.1 M + u.(1 -1 M ) where M = {x ∈ Ω; | u(x) -u(x)| > σ(x)}. end Return u N it +1
size, contrast, compression level and the way they are cropped 3 . In consequence, it is sometimes difficult to compare performance results provided by different papers, which may have used different image versions: on two different 512 × 512 versions of Lena, we observed PSNR differences of more than 3dB for the same algorithm. In this paper, we always run the different algorithms on the same images and noise samples (when the code is available), instead of using the performance results provided in the literature. We also provide the different images used in the paper on our website http://perso.telecom-paristech.fr/ ~delon/Demos/Impulse/.

Let us recall that the PSNR is a common way to measure the quality of a restored image v in comparison to the undamaged one u 0 . It is given by the formula

PSNR(u 0 , v) = 10 log 10 255 2 |Ω| x∈Ω (u 0 (x) -v(x)) 2 ,
where |Ω| is the size of the support of u 0 . An alternative to this quality measure is given by the Structural Similarity, an index (SSIM) measuring the similarity between images [START_REF] Wang | Image quality assessment: From error visibility to structural similarity[END_REF], designed to be more consistent with human perception than the usual measures PSNR or MSE. The code used here to compute SSIM is the one kindly provided by the authors of [START_REF] Wang | Image quality assessment: From error visibility to structural similarity[END_REF] on http://www.cns.nyu. edu/ ~lcv/ssim/, with the default parameters. For the sake of completeness, we also provide the Mean Absolute Error for all results:

MAE(u 0 , v) = x∈Ω |u 0 (x) -v(x)| |Ω| ,
The results of Tables 6.1, 6.2 and 6.3 are all obtained on the same noise samples, using the codes kindly provided by the authors of the papers, both for the trilateral filter 4 and for ROLD-EPR. The results for ROAD-trilateral [START_REF] Garnett | A universal noise removal algorithm with an impulse detector[END_REF] are obtained with the parameters set as recommended by the authors, i.e. σ s = 1, σ I = 40, σ J = 30 and 4 iterations of the algorithm. We observed that the parameter σ R of the trilateral filter, which controls its radiometric influence, has a great impact on the performances of the code and that its optimal value strongly depends on the noise level. For this reason, we keep for each experiment the parameter σ R which yields the best PSNR (the same value of σ R is used for the MAE and SSIM tables). For ROLD-EPR, we iterate the algorithm 7 times and keep all parameters as proposed in [START_REF] Dong | A detection statistic for random-valued impulse noise[END_REF]. In the PSNR table, we also put into brackets the best PSNR of the algorithm along the iterations. Table 6.1 shows quite similar performances between PARIGI and ROLD-EPR on most images (the SSIM results are particularly close), with a slight advantage of the patch-based approach for small noise ratios, and the opposite for large noise ratios. This proves that a patch-based approach is well founded for random-valued impulse noise removal. Much larger differences in PSNR, MAE and SSIM can be observed on the image Barbara. As we will see in the following, these remarkable differences in performance can be mostly explained by the way the three denoising approaches are handling regular textures in images. Let us insist on the fact that the parameters of our approach (patch size, research neighborhood size, number of iterations, ROAD threshold, number of trusted patches in function of p) are fixed for all these experiments, as explained in Section 5.2. Better results can be obtained by optimizing these values for each image.

In order to evaluate the precision of the PSNR results, for each value of p and each method, we compute the PSNR for 10 different noise realizations. In practice, the standard deviation of the PSNR for all methods is always close to 0.1: for instance, for the scheme PARIGI on the image Lena and p = 0.4, the empirical standard deviation is 0.096. This means that if we observe, between two methods, PSNR differences smaller than 0.2, these differences can be considered as meaningless in practice.

Figures 6.1, 6.3, 6.2 and 6.4 provide a visual comparison of the different approaches on the different images when p = 20% and p = 40% (images should be seen at full resolution on the electronic version of the paper). When the impulse noise ratio is low, all methods are efficient and results are visually quite similar in smooth regions. Now, fine geometrical structures (for instance the ropes on the Boat image) are clearly better handled by our patch-based scheme. Conversely, ROLD-EPR seems to preserve better the image grain than other methods.

When the noise ratio increases, one can observe that the result of PARIGI is clearly smoother than the ones of ROLD-EPR [START_REF] Dong | A detection statistic for random-valued impulse noise[END_REF] and ROAD-Trilateral [START_REF] Garnett | A universal noise removal algorithm with an impulse detector[END_REF], in which some clues of impulse noise remain. On Barbara, while all schemes yield reasonable results on constant regions, our approach is the only one to handle properly the regular stripes of the clothes, while the trilateral filter and ROLD-EPR replace them by mottled textures. Table 6.1 PSNR results of different restoration filters for the 512×512 images Lena, Bridge, Baboon, Barbara, Cameraman, Boat, Peppers and Goldhill for pure impulse noise. For the results of the trilateral filter, the parameter σR is optimized between 10 and 50 for each experiment. As for ROLD-EPR, we show into brackets the best possible PSNR obtained along the iterations of the algorithm. Other parameters are set as explained in the text.

ROAD-Trilateral

6.2. Gaussian and impulse noise mixture. We now investigate the performance of PARIGI for the denoising of images suffering from a mixture of Gaussian and impulse noise, by confronting it with the recent approaches [START_REF] Garnett | A universal noise removal algorithm with an impulse detector[END_REF][START_REF] Xiao | Restoration of images corrupted by mixed Gaussian-impulse noise via l 1 -l 0 minimization[END_REF]. To the best of our knowledge, the paper of Xiao et al. [START_REF] Xiao | Restoration of images corrupted by mixed Gaussian-impulse noise via l 1 -l 0 minimization[END_REF] on Gaussian-impulse noise can be considered as the state of the art on the subject, and [START_REF] Garnett | A universal noise removal algorithm with an impulse detector[END_REF] is one of the first paper proposing to handle such a noise mixture. To be completely fair, we do not include the results of ROLD-EPR in the comparison, mainly because it was not developed to handle such a noise mixture and does not work well as soon as the Gaussian standard deviation increases. Nevertheless, we think that a similar variational approach could be adapted to this form of degradation. Table 6.4 shows the PSNR results of the three methods for different combinations of Gaussian and impulse noise. The results for the l 1 -l 0 minimization approach of Xiao et al. are taken from Tables 8 and9 in [START_REF] Xiao | Restoration of images corrupted by mixed Gaussian-impulse noise via l 1 -l 0 minimization[END_REF], and thus come from different noise samples. The results of the trilateral filter are computed on the same noise samples as our approach. Again, for the trilateral filter, we choose for each mixture experiment the parameter σ R between 10 and 50 which yields the best PSNR, the optimal σ R begin strongly dependent on the noise ratio. Again, the PSNR standard deviation ROAD-Trilateral [START_REF] Garnett | A universal noise removal algorithm with an impulse detector[END_REF] for one experiment is generally around 0.1, which means that PSNR differences of this order of magnitude should not be considered as meaningful. The performances of the different schemes are quite similar, with the notable exception of Barbara. As is the case for the pure impulse noise experiments, the performance differences on this image can be explained by the presence of several regular textures, better preserved by our patch-based approach. These results are also illustrated on Figures 6.5 and6.6. 6.3. Discussion. It is remarkable that the patch-based method presented in this paper provides such good performances on Barbara, when compared to state of the art approaches. Such a difference cannot be observed on other images. The visual experiments suggest that this difference is mostly due to the way the different methods are handling the regular stripes on Barbara's clothes. This is confirmed by the experiments shown on Figure 6.7. In this example, we compare several denoising procedures on a synthetic 256 × 256 image composed of perfect vertical stripes and suffering from pure impulse noise, with p = 50%. Thanks to the huge patch redundancy of this particular image, PARIGI is able to recover the vertical stripes almost perfectly after a few iterations, providing a result very close the original image (because of the image regularity, the result is even better by choosing larger patches). The variational approach of [START_REF] Dong | A detection statistic for random-valued impulse noise[END_REF] and the trilateral ROAD-Trilateral [START_REF] Garnett | A universal noise removal algorithm with an impulse detector[END_REF] filter [START_REF] Garnett | A universal noise removal algorithm with an impulse detector[END_REF] do not fully take advantage of this redundancy, and this explains their quite poor results on this kind of image. On the contrary, we can observe on Figure 6.1 that the edge preserving regularization of [START_REF] Dong | A detection statistic for random-valued impulse noise[END_REF] seems to preserve the fur of the mandrill image slightly better than our patchbased approach. An interesting direction of research, which is beyond the scope of this paper, would consist in identifying the image classes (or, at least, the texture classes) optimally restored by each kind of denoising approach.

Conclusion.

In this paper, it is shown that a patch-based approach can be an efficient tool to remove mixtures of Gaussian and impulse noises. This result was known (and has been widely studied) in the particular case of Gaussian noise, but its extension to impulse degradations necessitates a careful choice of both the similarity measure between patches and the statistical estimator of the original patches. It is shown on several experiments that this patch-based approach permits to attain state of the art denoising performances on classical images. An important performance gain is demonstrated on geometrically regular textures. This work opens several perspectives. First, as explained in the previous section, we intend to explore the connections between image or texture regularity models and the denoising performances of different kinds approaches (local filters, variational methods, patch-based approaches). In particular, it would be of great interest to determine

1

  PARIGI stands for Patch based Approach for the Restoration of Images affected by Gaussian and Impulse noise.

Fig. 3 . 1 .

 31 Fig. 3.1. Comparison of different distances on noisy patches. From top to bottom: an edge patch and a T-junction patch.From left to right: the original patch (marked by a white window of size 7 × 7 pixels and its neighborhood V ), a noisy version of the image (p = 0.2 and σ = 10), and the mean number of selected patches for the different distances when σ = 10, p = 0.2 (third column) and p = 0.5 (fourth column). The 2 distance corresponds to blue curves, the trimmed distance D trimmed to green curves, the weighted distance D weighted to red curves and the "impulse-controlled" (IC) distance used by[START_REF] Li | A new method for removing mixed noises[END_REF] to indigo curves. See the text for some comments on these graphics.

. 8 )Figure 4 .

 84 Figure 4.1 shows the values of R[X ( n+12 ) ] for different values of p and n, with M = 255 and µ = 20 or µ = 100. For p = 0.5 and µ = 100, n = 18 samples are necessary to make the square root of the MSE decrease below 5, which is quite reasonable. However, for µ = 20, this number of necessary samples becomes larger than 100, as can be observed on the left part of Figure4.1. This huge difference in behavior is due to the bias of the median estimator which attains very large values as soon as p ≥ 0.5 and for gray levels µ far from the middle of [0, M ]. As a conclusion, the median estimator is clearly more interesting than the one built upon the mean, but in practice, its quadratic risk is not controlled, except for small values of p.

Fig. 4 . 1 .

 41 Fig. 4.1. Square root of the risk (or MSE) R[X ( n+1 2 ) ] as a function of p (between 0.1 and 0.8 on the vertical axis) and n (between 1 and 100 on the horizontal axis), with M = 255, µ = 20 on the left and µ = 100 on the right.

Fig. 4 . 2 .

 42 Fig. 4.2. Left: parameter qn,p of the distribution (4.9) in function of p (vertical axis) and n (horizontal axis). The white curve corresponds to the level line qn,p = 10 -3 . This line yields the number n of samples necessary to keep the square root of the risk always below 5 gray levels (whatever the value of µ). Right: Square root of the risk E[| Xn -µ| 2 ] as a function of p (vertical axis) and n (horizontal axis), for µ = 100.

Figure 4 .

 4 3 shows the empirical values of R[ X n ] for p ∈ {0.1, . . . , 0.8} and n ∈ {5, . . . , 100}. The four subfigures correspond respectively to σ = 5, 10, 15 and 20.

= 20 Fig. 4 . 3 .

 2043 Fig. 4.3. Square root of the risk R[ Xn] for the MLE in the mixture case, as a function of p (vertical axis) and n (horizontal axis), for µ = 100. White crosses correspond to the minimum values of n ensuring that the square root of the risk is smaller than 5.

  Fig. 5.1. Left: statistical results of the estimation of the impulse noise parameter p on the Lena image when there is no Gaussian noise. The results p are shown as boxplot graphics, where the horizontal axis represents the tested values of p: from 0.1 to 0.9. Boxes are the statistics obtained from 100 samples for each value of p. On the right, same experiment made on the Lena image in the case of a mixed noise: impulse and Gaussian with σ = 10.

Fig. 5 . 2 .

 52 Fig. 5.2. Neighborhood Vx and possible corresponding patches.

Fig. 6 . 1 .

 61 Fig.6.1. Comparative results on Lena,Bridge,Baboon and Barbara with p = 20% of random-valued impulse noise. For each column, from top to bottom: original image, noisy image, ROAD+trilateral[START_REF] Garnett | A universal noise removal algorithm with an impulse detector[END_REF], ROLD+EPR[START_REF] Dong | A detection statistic for random-valued impulse noise[END_REF], PARIGI. Images should be seen at full resolution on the electronic version of the paper.

Fig. 6 . 2 .

 62 Fig.6.2. Comparative results on Lena,Bridge,Baboon and Barbara with p = 40% of random-valued impulse noise. For each column, from top to bottom: original image, noisy image, ROAD+trilateral[START_REF] Garnett | A universal noise removal algorithm with an impulse detector[END_REF], ROLD+EPR[START_REF] Dong | A detection statistic for random-valued impulse noise[END_REF], PARIGI. Images should be seen at full resolution on the electronic version of the paper.

Fig. 6 . 3 .

 63 Fig. 6.3. Comparative results on Cameraman,Boat,Peppers and Goldhill with p = 20% of random-valued impulse noise. For each column, from top to bottom: original image, noisy image, ROAD+trilateral [25], ROLD+EPR [21], PARIGI. Images should be seen at full resolution on the electronic version of the paper.

Fig. 6 . 4 .

 64 Fig.6.4. Comparative results on Cameraman,Boat,Peppers and Goldhill with p = 40% of random-valued impulse noise. For each column, from top to bottom: original image, noisy image, ROAD+trilateral[START_REF] Garnett | A universal noise removal algorithm with an impulse detector[END_REF], ROLD+EPR[START_REF] Dong | A detection statistic for random-valued impulse noise[END_REF], PARIGI. Images should be seen at full resolution on the electronic version of the paper.

15 Fig. 6 . 5 .

 1565 Fig. 6.5. Comparative results on Lena and Barbara with p ∈ {10%, 30%} of random-valued impulse noise and σ ∈ {5, 15}. For each subfigure, from left to right: noisy image, result of the trilateral filter [25], PARIGI.

15 Fig. 6 . 6 .

 1566 Fig. 6.6. Comparative results on Cameraman and Boat with p ∈ {10%, 30%} of random-valued impulse noise and σ ∈ {5, 15}. For each subfigure, from left to right: noisy image, result of the trilateral filter [25], PARIGI.

74 Fig. 6 . 7 .

 7467 Fig. 6.7. Comparison of different denoising approaches on the image Stripes, with p = 50%.

Table 4 . 1

 41 Number n of patches used in the maximum-likelihood estimation for each value of p in the case of pure impulse noise.

Table 4 . 2

 42 Number n of patches used in the maximum-likelihood estimation for each value of p in the case of a mixture noise.

Algorithm 1 :

 1 Denoising algorithm PARIGI input : Image u, half patch size f , half size of the research zone t, number of iterations N it (N it = 2 by default) output: Denoised image u N it +1 Initialize u 1 = u; for i from 1 to N it do Estimate p on u i ; // See Section 5.1 , the subset of V x corresponding to the n nearest neighbors of P x for D weighted in u i ; for each δ ∈ Ω f do L x+δ = L x+δ ∪ {u i (y + δ), y ∈ V n x }; end end for each pixel x do Obtain ( u(x), σ(x)) as the MLE from the set L x ; // see Section 4.2.2.

	Compute n = N umberN earestN eighbors(p); // see Tables 4.1 and 4.2 in
	Section 4.2.2.
	for each pixel x do
	L x = ∅;
	end
	for each pixel x do
	Find V n x end
	Compute
	u i+1

Table 6 .

 6 2MAE results of different restoration filters for the 512 × 512 images Lena, Bridge, Baboon, Barbara, Cameraman, Boat, Peppers and Goldhill for pure impulse noise. For the results of the trilateral filter, the parameter σR is optimized between 10 and 50 for each experiment. Other parameters are set as explained in the text.

Table 6 .

 6 3SSIM results of different restoration filters for the 512 × 512 images Lena, Bridge, Baboon, Barbara, Cameraman, Boat, Peppers and Goldhill for pure impulse noise.

	ROLD-EPR [21] PARIGI

B(n, k, q) = n i=k n i q i (1 -q) n-i .

See for instance http://www.ece.rice.edu/ ~wakin/images/ for a discussion on the standard gray level version of Lena. This standard version is the one used in all our experiments.

The authors of[START_REF] Garnett | A universal noise removal algorithm with an impulse detector[END_REF] provide their code on their webpage www.ssc.wisc.edu/ ~thuegeri.

ROAD-Trilateral

 [START_REF] Garnett | A universal noise removal algorithm with an impulse detector[END_REF]Xiao [START_REF] Xiao | Restoration of images corrupted by mixed Gaussian-impulse noise via l 1 -l 0 minimization[END_REF] which image class is optimal for each denoising method. Another possible extension of this work includes the collaboration of different denoising approaches, depending on the image local regularity and redundancy. Finally, we intend to study the theoretical bounds of impulse noise removal, as pioneered by [START_REF] Chatterjee | Is denoising dead?[END_REF] and [START_REF] Levin | Natural image denoising: Optimality and inherent bounds[END_REF] in the case of Gaussian noise.