Charles-Alban Deledalle
email: charles-alban.deledalle@telecom-paristech.fr

Vincent Duval
email: vincent.duval@telecom-paristech.fr

Joseph Salmon
email: salmon@math.jussieu.fr

Non-Local Methods with Shape-Adaptive Patches (NLM-SAP)

Keywords: Image denoising, non-local means, spatial adaptivity, aggregation, risk estimation, SURE

Noise (AWGN) has drawn a lot of efforts from the image processing community. A wide variety of strategies were proposed during the last decades, from partial differential equations (PDE) to transform-domain methods such as wavelets, approximation theory or stochastic analysis.

In fact, a first major difficulty in image denoising is to handle efficiently regular parts while preventing edges from being blurred. To this end, diffusion methods [START_REF] Perona | Scale space and edge detection using anisotropic diffusion[END_REF] and variational methods [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] have shown their effectiveness. The latter approaches introduce a prior model of the noise-free image and search for a trade-off between data fitting and regularization. It is crucial for these techniques to define a suitable prior that guarantees both the smoothness of the denoised image and the preservation of its structure. Total variation [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] is an example of a prior that enforces smoothness while preserving edges. However, such priors tend to bias the denoised image [START_REF] Strong | Edge-preserving and scaledependent properties of total variation regularization[END_REF], especially in texture area or when high noise levels are considered. Donoho and Johnstone [START_REF] Donoho | Adapting to unknown smoothness via wavelet shrinkage[END_REF] have shown that spatial adaptive methods can be obtained by soft-thresholding wavelet decompositions. Wavelet shrinkage can be interpreted from a Bayesian perspective. It results from an ℓ 1 sparsity-promoting prior: the regularity of the signal is enforced by using sparsifying priors on the coefficients in the transform domain. Indeed, wavelet and cosine bases are well known to be able to capture most of a signal/image with very few coefficients, a good property for compression. For denoising applications, such orthogonal transforms lead to a separation of signal and noise. Noise can then be strongly suppressed by canceling the least significant coefficients. These types of methods culminated with the BLS-GSM algorithm [START_REF] Portilla | Image denoising using scale mixtures of gaussians in the wavelet domain[END_REF]. Other improvements of wavelet oriented methods have been proposed such as Shape Adaptive Domains [START_REF] Foi | Pointwise shapeadaptive DCT for high-quality denoising and deblocking of grayscale and color images[END_REF] and over-complete [START_REF] Starck | The curvelet transform for image denoising[END_REF] or geometric [START_REF] Pennec | Sparse geometric image representations with bandelets[END_REF] decompositions.

Spatial adaptivity can also be reached by considering adaptive neighbor filters, as the well known Sigma Filter [START_REF] Lee | Digital image smoothing and the sigma filter[END_REF], Yaroslavsky Filter [START_REF] Yaroslavsky | Digital picture processing[END_REF] and Bilateral Filter [49]. Moreover, applications of Lepski's method [START_REF] Lepski | Optimal spatial adaptation to inhomogeneous smoothness: an approach based on kernel estimates with variable bandwidth selectors[END_REF] (known in signal processing as the ICI rule, Intersection of Confidence Intervals) to kernel smoothing are also known to perform well for denoising. Their theoretical performance is as good as the best wavelet shrinking methods: they reach adaptive rate in the minimax sense [START_REF] Goldenshluger | On spatially adaptive estimation of nonparametric regression[END_REF] for regular functions. In practice, Lepski's rule has successfully been applied to image denoising with Adaptive Weights Smoothing [START_REF] Polzehl | Adaptive weights smoothing with applications to image restoration[END_REF] or Shape Adaptive Windows [START_REF] Katkovnik | Directional varying scale approximations for anisotropic signal processing[END_REF] for instance.

Though very efficient at dealing with edges, such methods cannot proceed efficiently in textured regions. To overcome this drawback, many authors have proposed to work with small sub-images, called patches, to take into account the redundancy in natural images, especially in textured parts. The interest of using patches lies in their robustness to noise. The NL-Means (Non-Local Means) [5] and UINTA [START_REF] Awate | Unsupervised, informationtheoretic, adaptive image filtering for image restoration[END_REF] algorithms are typical examples of this approach, as is their extension using Lepski's method [START_REF] Kervrann | Optimal spatial adaptation for patch-based image denoising[END_REF]. Those algorithms consist in averaging similar pixels by measuring their similarity through patches. Iterative versions have also been proposed [START_REF] Brox | Efficient nonlocal means for denoising of textural patterns[END_REF][START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF][START_REF] Deledalle | Iterative weighted maximum likelihood denoising with probabilistic patchbased weights[END_REF] as well as mixed methods using Total Variation [START_REF] Louchet | Total variation as a local filter[END_REF].

Another approach proposed by Aharon et al. [START_REF] Aharon | K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation[END_REF] and extended in [START_REF] Mairal | Learning multiscale sparse representations for image and video restoration[END_REF][START_REF] Mairal | Non-local sparse models for image restoration[END_REF] also uses patches to denoise the image. Though, for these methods, the key point is to get a good representation for each patch of the image by using dictionary learning techniques, well-suited for overcomplete representation (based for instance on ℓ 1 regularization or greedy algorithms). This leads to a state-of-the-art denoising method called Learned Simultaneous Sparse Coding (LSSC) [START_REF] Mairal | Non-local sparse models for image restoration[END_REF]. It has to be noted that this algorithm incorporates a prior on images by initializing a dictionary of patches on a large clean dataset. Another state-of-the-art method in denoising is BM3D [7]. It also relies on patches and consists of a smart combination of classical filtering techniques, such as wavelet denoising and Wiener's Filter, applied in the space of patches.

These methods are already quite efficient at dealing with smooth regions and textures. However, since they use patches with a fixed square shape and a fixed scale over the whole image, their performances may be limited for dealing with edges, mostly for high contrasted edges. Indeed, high contrasted edges present few redundancies and their denoising version suffer from a persistence of residual noise: this is named the noise halo. In order to overcome this drawback, more directional priors may be considered, using locally chosen scales and orientations of shapes. Few attempts have been made to use several patch sizes (see [START_REF] Mairal | Learning multiscale sparse representations for image and video restoration[END_REF] for learning with patches or [43] for the NL-Means).

As far as we know, the only work trying to handle variable shapes rather than simple square has recently been proposed by Dabov et al. [START_REF] Dabov | BM3D image denoising with shape-adaptive principal component analysis[END_REF] as a way to improve the BM3D algorithm. The authors propose to adapt the shapes used by the algorithm: they locally select a shape by applying Lepski's method, and then perform the same steps of the BM3D algorithm with these shapes rather than with common square blocks (i.e., patches).

Recently, spatial-adaptive methods for selecting the parameters of NL-means have been proposed in [16,[START_REF] Doré | Robust NL-Means Filter With Optimal Pixel-Wise Smoothing Parameter for Statisti-cal Image Denoising[END_REF]. Both methods propose to locally select the parameters which minimize a local estimate of the risk (i.e., the Mean Square Error, MSE) by considering respectively Stein's Unbiased Risk Estimate (SURE) [START_REF] Stein | Estimation of the mean of a multivariate distribution[END_REF][START_REF] Stein | Estimation of the mean of a multivariate normal distribution[END_REF] or the C p criterion [START_REF] Mallows | Some comments on c p[END_REF]. The use of SURE for NL-Means was originally proposed in order to select the best bandwidth parameter [START_REF] Van De Ville | SURE-based Non-Local Means[END_REF]. Applications of SURE emerged for choosing the smoothing parameter in families of linear estimates [START_REF] Li | From Stein's unbiased risk estimates to the method of generalized cross validation[END_REF] such as for model selection, ridge regression, smoothing splines, etc. It was then widely used in the wavelet community after the introduction of the SURE-Shrink algorithm [START_REF] Donoho | Adapting to unknown smoothness via wavelet shrinkage[END_REF]. Solo [START_REF] Solo | A sure-fired way to choose smoothing parameters in ill-conditioned inverse problems[END_REF] gave a general form of SURE for an estimator defined as a minimizer of a regular energy, especially for least square regression regularized by a Sobolev norm or the Total Variation. More recently, linear combinations of estimates based on SURE were considered [START_REF] Blu | The SURE-LET approach to image denoising[END_REF] instead of the selection of a single one. Moreover, Ramani et al. [START_REF] Ramani | Monte-Carlo SURE: a black-box optimization of regularization parameters for general denoising algorithms[END_REF] have proposed a Monte Carlo approach to evaluate SURE when a closed-form expression is not available or too computer-intensive. Our paper is in the same vein as [START_REF] Van De Ville | SURE-based Non-Local Means[END_REF]16] and proposes to locally select or aggregate the best shapes.

Our contributions-The aim of this paper is to investigate the potential benefit of replacing the simple square patches with more general shapes, in the classical NL-Means filter. To this end, we propose a fast algorithm, based on the fast Fourier transform, which allows to compute the solution of the NL-Means for arbitrary patch shapes. Then we explain how to combine the estimators associated with each shape in a suitable way. We select or combine locally the shape-based estimates by measuring the performance of their associated denoisers with SURE. We coin such type of algorithms Non-Local Means with Shape-Adaptive Patches (NLM-SAP).

The main advantage of using adaptive patch shapes in the context of NL-Means is to reduce the noise halo produced around high contrasted edges. Such edges can appear anytime in natural images and they appear commonly in high dynamic range images (HDR). It has been shown that high contrasts in synthetic aperture radar images interfere with Non-Local approaches since these kinds of images present very high contrasted features (known as bright scatterers) [START_REF] Deledalle | NL-InSAR: Nonlocal interferogram estimation[END_REF]. Our method is an improvement of the NL-Means taking into account the anisotropy of natural images. It is all the more relevant when the images to denoise present high contrasted edges for which the classical NL-Means fails.

Organization of the paper -We give in Section 2 a general overview of the NL-Means method. We introduce in Section 3 a more general framework using general shapes instead of square patches to measure the similarity between pixels. This leads in Section 4 to the natural problem of locally selecting or combining the best shapes in our NLM-SAP algorithm. In Section 5 we illustrate numerically, and above all visually, the gain in aggregating various shape-based estimates in a proper manner.

An overview of the NL-Means

In this section we define the model we use for the noisy image, and we recall the definition and basic facts about NL-Means. We also present some limitations and justify the need to introduce variable shapes instead of the standard fixed square patches.

The AWGN model

In this paper we are concerned with the problem of denoising. A noisy image Y is given and assumed to be a version of an unobserved deterministic image f corrupted by an Additive White Gaussian Noise (AWGN).

Let Ω ⊂ Z 2 be the (bounded) indexing set of the pixels. For any pixel x in the grid Ω, the model has the following formulation:

Y(x) = f (x) + ε(x) , (1)
where ε is a centered Gaussian noise with known variance σ 2 and the noise components ε(x) are independent. Even if we only focus on the case of AWGN, our methodology could be adapted to other noise models by considering NL-Means extensions such as [START_REF] Kervrann | Bayesian nonlocal means filter, image redundancy and adaptive dictionaries for noise removal[END_REF][START_REF] Deledalle | Iterative weighted maximum likelihood denoising with probabilistic patchbased weights[END_REF] and statistical results on unbiased risk estimation [START_REF] Hudson | A natural identity for exponential families with applications in multiparameter estimation[END_REF].

Definition and parameters

In this section we present the usual definition of the NL-Means procedure as introduced by Buades et al. [5]. It is an averaging filter, i.e., for each pixel the output of the procedure is a weighted average of the whole image. The weights used are selected using a "metric" which determines whether two pixels are similar or not. For instance, with the Bilateral Filter [49], the "metric" consists in comparing both pixels values and pixel locations. The core idea of the NL-Means is to create a metric governed by patches surrounding each pixel, regardless of their position, i.e., non-local in the image space. For any pixel of interest x, an estimate of the value f (x) is given by:

f (x) = x ′ ∈Ω ω(x, x ′ )Y(x ′ ) x ′ ∈Ω ω(x, x ′ ) , (2)
where the weights ω(x, x ′ ) depend on patches around x and x ′ . The denominator is a normalizing factor which ensures the weights sum to one. For a fixed (odd) width p, a patch P x is a subimage of width p, centered around the pixel of interest x:

P x = Y(x + τ ), τ ∈ - p -1 2 , p -1 2 2 . (3)
With this notation, the original weights in the NL-Means are of the following form:

ω(x, x ′ ) = ϕ P x -P x ′ 2 2,a 2h 2 , ( 4
)
where h > 0 is the bandwidth parameter, ϕ is the kernel used to measure similarity between patches, • 2,a is a weighted Euclidean norm using a Gaussian kernel, and a is the bandwidth that controls the concentration of the kernel around the central pixel. Let us now briefly recall the influence of each parameter (see [16] for a more extensive discussion on this subject).

The bandwidth h: this parameter has a smoothing effect and plays the same role as the bandwidth for kernel methods in statistics (cf. [START_REF] Wasserman | All of Nonparametric Statistics[END_REF] for more details). With our parametrization, the larger the bandwidth, the smoother the image becomes. Choosing this parameter is a difficult task and many solutions have been proposed in the literature. The simplest and most common one is to set a single h for the whole image, whose value is determined by cross validation on a small dataset of images. In [START_REF] Polzehl | Propagation-separation approach for local likelihood estimation[END_REF][START_REF] Kervrann | Optimal spatial adaptation for patch-based image denoising[END_REF], the authors set this parameter according to the quantile of a χ 2 distribution, due to the particular metric they consider to compare pixels (or patches). Van De Ville and Kocher [START_REF] Van De Ville | SURE-based Non-Local Means[END_REF] calculate an unbiased risk estimate of the NL-Means to globally select the bandwidth h. Pursuing this idea, Duval et al. [16] consider a method based on the same approach but to locally select the bandwidth.

The main interest of the last two approaches is that they give a closed-form expression of the SURE for NL-Means.

The search window size ℓ: the summation defined in Equation ( 2) is usually restricted to an ℓ × ℓ search window around the pixel of interest.

Defining

Ω ℓ = -ℓ-1 2 , ℓ-1 2 2
, the search window centered on each pixel x is then x + Ω ℓ . This was proposed in the seminal work [5] for computational acceleration. Though, some authors have also noticed that choosing locally the best search windows [START_REF] Kervrann | Optimal spatial adaptation for patch-based image denoising[END_REF] or restricting the average over small ones (see [16] or [START_REF] Salmon | On two parameters for denoising with Non-Local Means[END_REF] for more details) could benefit to the NL-Means procedure.

The kernel ϕ: function ϕ was chosen by Buades et al. [5] as t → exp(-t), but other choices may be considered, such as compactly supported smooth functions. It was noticed by several authors [START_REF] Goossens | An improved non-local denoising algorithm[END_REF] that weights with compact support yield better results. In this paper we have restricted our experiments to the kernel given in [5], but any other C 1 kernel could be handled by our framework.

The patch size p: this parameter is generally chosen equal to 5, 7 or 9. Using a width p = 1 would lead to a method close to the Bilateral Filter [49] or Yaroslavsky Filter [START_REF] Yaroslavsky | Digital picture processing[END_REF]. Moreover, in most papers the patch width is a global fixed parameter. Few works have tried to handle the difficult task of using several sizes of patches for a single image. To our knowledge, the first attempt was proposed in the context of learning patches by [START_REF] Mairal | Learning multiscale sparse representations for image and video restoration[END_REF], using Support Vector Machines (SVM).

Our approach extends the notion of patches to the more general concept of shapes. Our goal is to better handle anisotropy in natural images using oriented patches. Indeed, it is well known that edges are not satisfactorily denoised when using square patches with the NL-Means. This leads to a noise halo around edges (see Fig. 1). We propose in the next section an algorithm to compute efficiently the solution of NL-means for arbitrary shapes. We also consider in details several natural collections of shapes. Numerical results are given in Section 5.

3 From patches to shapes: beyond the rare patch effect

In practice, we have seen that the original algorithm suffers from a noise halo around edges, due to an abrupt lack of redundancy of the image. This phenomenon is sometimes referred to as the rare patch effect (Fig. 1). Statistically, it leads to an NL-Means estimator with large variance around edges. Several solutions have already been proposed to handle this drawback [43,16,[START_REF] Louchet | Total variation as a local filter[END_REF] (see Section 3.2 for more details). Here, we generalize these approaches by considering general shapes rather than simple square patches (cf. Fig. 2).

In order to deal with patches of arbitrary shapes, we reformulate the way the distance between two pixels is measured in terms of patches. The weighted Euclidean distance • 2,a used in Eq. ( 4) can be generalized using the following expression:

d 2 S (x, x ′ ) = τ ∈Ω S (τ ) (Y(x + τ ) -Y(x ′ + τ )) 2 , ( 5
)
where S encodes the shape we aim at. With this notation we can easily rewrite the original NL-Means with a simple S, cf. Eq. ( 8) and Eq. ( 9). One of our contributions is to provide an efficient algorithm, based on the Fast Fourier Transform (FFT), to compute the distances in Eq. ( 5). Our implementation is independent of the shape, and can thus be applied with different shapes (see Section 3.1 for details). As soon as we consider the use of anisotropic shapes, and not just squares or disks centered on the pixel of interest, two questions emerge. The first one is how to choose the collection of shapes to consider. The second issue is to propose a way to combine the estimators provided by each shape.

The collection of shapes should be composed of more than one shape to locally take into account the geometrical properties of natural images. Consider for instance the use of a single vertically elongated patch. This could be interesting to handle vertical features, but eventually we would not optimally deal with horizontal details with only one oriented shape. So, the collection should be diversified and numerous enough to identify directional features (see Fig. 8 for a visual illustration). At the same time, it should remain small enough so that the algorithm is not computationally intensive.

Fast algorithm to handle shapes

In this section, we present a fast way to compute the NL-Means weights for general shapes, based on the 2D-FFT. It is inspired from works initiated in [START_REF] Wang | Fast non-local algorithm for image denoising[END_REF] and [START_REF] Darbon | Fast nonlocal filtering applied to electron cryomicroscopy[END_REF] to speed up the NL-Means algorithm. However, contrary to these approaches, ours can deal with nonsquare and/or non-binary patches, i.e., with general shapes S. Like them, our method is independent of the shape size.

In [START_REF] Wang | Fast non-local algorithm for image denoising[END_REF] and [START_REF] Darbon | Fast nonlocal filtering applied to electron cryomicroscopy[END_REF], the authors propose to compute the Euclidean distances using "Summed Area Tables" [START_REF] Crow | Summed-area tables for texture mapping[END_REF] (also called "Integral Images" [START_REF] Viola | Rapid object detection using a boosted cascade of simple features[END_REF]). This allows them to reduce the computational cost of the NL-Means from

O(|W | • |Ω| • |P|) to O(|W | • |Ω|), where |W | = ℓ 2
is the number of pixels in the search window, |Ω| is the image domain size and |P| is the patch size (we refer to [START_REF] Darbon | Fast nonlocal filtering applied to electron cryomicroscopy[END_REF] for more details). To compute these integral images, the authors change the original algorithm by swapping the two "for" loops: instead of considering all the shifts for each pixel, they consider all the pixels for each shift.

We use basically the same swapping trick. Notice that Equation ( 5) can be reformulated for any translation parameter δ (i.e., taking x ′ = x + δ) as a discrete convolution:

d 2 S (x, x + δ) = τ ∈Ω S (τ ) (Y(x + τ ) -Y(x + δ + τ )) 2 = ( Š ⋆ ∆ δ )(x) , (6)
where (

S (τ ) = S (-τ ), ∆ δ (x) = (Y(x) -Y(x + δ)) 2
and ⋆ is the convolution operator. This term can be interpreted as the correlation between the shape S and Compute the 2D-FFT F (∆ δ ) Perform the convolution of ∆ δ by the shape Š

d 2 S (•, • + δ) ← " F -1 " F (S)F (∆ δ ) "" (•) ⊲ O(|Ω| • log |Ω|) operations using 2D-FFT
for all pixels x in Ω do Compute the weights

ω(x, x + δ) = ϕ " d 2 S (x, x + δ) 2h 2 « Update the accumulators A(x) ← A(x) + ω(x, x + δ)Y(x + δ) B(x) ← B(x) + ω(x, x + δ)
end for end for Final (normalized) estimator f (x) = A(x) B(x) for all pixel x Note: the central pixel (δ = 0) is treated as a special case ⊲ see Section 5 for details the square difference of the observe image and the δshifted version. The convolution Š ⋆ ∆ δ can be computed quickly thanks to following relation:

Š ⋆ ∆ δ = F -1 (F( Š)F(∆ δ )) = F -1 (F(S)F(∆ δ )) , ( 7
)
where F is the 2D discrete Fourier transform (2D-FFT) and F -1 is its inverse transform. According to Equation (7), and given a translation δ, we only need to perform one term by term multiplication in Fourier domain and two 2D-FFT (note that F(S) can be computed off-line). The repetition of this procedure for every translation δ covering the search window, leads to an algorithm (whose pseudo-code is detailed in Fig. 3) with a complexity of

O(|W | • |Ω| • log(|Ω|)).

Families of shapes

The main purpose of this paper is to show that the use of different shapes allows to reduce the rare patch effect. This point of view is a generalization of the NL-Means based on square patches with the reprojection studied in [43], since each translated patch can be regarded as a de-centered shape. Here, h is fixed and the challenge is to find shapes with enough similar candidates in the search window to reduce the noise. We now present several types of families we have considered. The first collections consist of classical squares and disks shapes. Then, we propose more directional shapes such as pie slices and bands displayed in Fig. 2.

Squares: To begin with, we apply our framework to the most commonly used shapes, i.e., the square shapes of odd length (so the squares have centers we can consider). For instance, choosing:

S (τ ) =    1, if τ ∞ ≤ p-1 2 , 0, otherwise, (8)
leads to the classical (simplified) NL-Means definition with square patches of size p × p and distance between patches measured by the Euclidean norm.

Gaussian: The original, but less common choice, is to set:

S (τ ) =    exp(-(τ 2 1 + τ 2 2 )/2a 2 ), if τ ∞ ≤ p-1 2 , 0, otherwise. (9)
Equation ( 9) means that the norm • 2,a is used to measure the distance between patches. This limits the influence of square patches corners and leads to a more isotropic comparison between patches.

Disks: Disk shapes are defined in the same way, using the Euclidean norm instead:

S (τ ) =    1, if τ 2 ≤ p-1 2 , 0, otherwise. (10)
A non-binary version may also be defined for pixels crossed by the boundary.

Pie slices:

We study a family of shapes, denoted as "pie", whose elements are defined with three parameters: two angles and a radius. These shapes represent a portion of a disk delimited by two lines and surrounding the discrete central pixel.

Bands: This family of shapes is simply composed of rectangles, potentially rotated and decentered with respect to the pixel of interest.

Connection with previous work

One of our main concern is to address the rare patch effect of the NL-Means algorithm. Different methods have been designed to limit this drawback or to improve the NL-Means in terms of quality or speed. The first attempt was proposed to speed-up the algorithm.

The idea is to denoise patch by patch rather than pixel by pixel. Taking into account patches overlaps, a fast implementation of the NL-Means is reached by using a sub-sampled grid of pixels [5,[START_REF] Kervrann | Bayesian nonlocal means filter, image redundancy and adaptive dictionaries for noise removal[END_REF]. Quality improvement can also be obtained by properly using overlapping patches. Indeed, we get |P| estimates for each pixel (where |P| is the number of elements in a patch). Some authors [5,[START_REF] Kervrann | Optimal spatial adaptation for patch-based image denoising[END_REF] simply propose to uniformly average those |P| estimates while a weighted average is performed in [43] (see Section 4.1 and 4.2). In our framework, these blockwise approaches are equivalent to use |P|, possibly decentered, square shapes (cf. Fig. 4).

Other methods have been introduced to reduce the noise halo. In [START_REF] Louchet | Total variation as a local filter[END_REF], the authors use a total variationbased pre-filtering of the image and set locally its parameter so that the NL-Means find enough similar patches. In [16], the authors aim to select locally the bandwidth parameter h using SURE (introduced in details in Section 4.3) to select enough patches according to a bias-variance trade-off.

Aggregation of shape-based estimates

In this section we investigate several ways to aggregate the NL-Means estimators based on different shapes of "patches". We have extended the standard square shape to other shapes such as disks, pies or bands (see Section 3.2). Thus, the new goal in this context is to determine how to locally take the most of each proposed denoiser.

Assume that for any pixel x in the image, we have built a collection of K pixel estimators f1 (x), • • • , fK (x) based on different shapes, and that estimates of their corresponding performances are available. With this information at hand, we can address different aggregation/optimization problems. Indeed, the ways of combining the estimators at hand may depend on the theoretical aggregation problem we aim to solve (as described by [START_REF] Nemirovski | Topics in non-parametric statistics[END_REF] and [START_REF] Tsybakov | Optimal rates of aggregation[END_REF]):

Uniformly weighted aggregation (UWA)

This is the simplest form of aggregation that we consider. The idea is to give the same weight to any shapebased estimator. With few shapes it is already an improvement in practice (both in term of PSNR and SSIM, see Table 3), but as the number of shapes increases, we can take into account positions that are completely irrelevant. It is in essence, the one proposed in [5] and referred as blockwise NL-means in [START_REF] Kervrann | Optimal spatial adaptation for patch-based image denoising[END_REF], with shapes being decentered square patches. Such an aggregation procedure still suffers from the rare patch effect (see Fig. 1, third line, for a visual illustration). We can express for any pixel x the uniform estimate fUWA (x) as the mean of the shape-based estimates:

fUWA (x) := 1 K K k=1 fk (x) . (11)
In our experiments, this method can be considered as the benchmark to be improved when using more complex aggregating strategies.

4.2 Variance-based decision, Weighted Average (WAV)

A possible way to limit the halo effect is to adapt WAV-reprojection [43] to general shapes. The idea, also proposed by Dabov et al. [7] in a different context, is to perform a weighted average of the estimates f1 (x), • • • , fK (x). To limit the noise halo, each weight should be chosen inversely proportional to the (approximate) variance of the corresponding estimator. In the context of NL-means, this approximate variance can be obtained in closed-form in the same way as in [START_REF] Kervrann | Optimal spatial adaptation for patch-based image denoising[END_REF], assuming that the coefficients ω(x, x ′ ) in Equation ( 2) can be treated as deterministic. Measuring the performances of the estimators in term of variance is well justified since the halo effect results to the high variance of our estimators around the edges (see [43]). However, it tends to over-smooth the edges and the thin details since it does not consider the bias of each estimator.

SURE-based decisions

A way to take the bias into account is to consider the risk estimate rather than the variance to locally attribute more weight to the estimators with small risks. In [START_REF] Van De Ville | SURE-based Non-Local Means[END_REF], Van De Ville and Kocher give a closedform expression of Stein's Unbiased Estimator of the Risk (SURE) for the NL-Means. They aim at selecting globally the best bandwidth for a given image. Here, our approach is different, despite the use of the same tool. Indeed, our choice of shape is done locally (i.e., for each pixel), since it is very unlikely that a single shape should be optimal for a whole natural image. Our method is closer to the one proposed by Duval et al. [16]. The authors rely on SURE to locally determine the parameters: the bandwidth h and the patch size p.

They have shown that a local choice of h reduces the visual artifacts, especially the rare patch effect. Let us now rephrase Stein's Lemma [START_REF] Stein | Estimation of the mean of a multivariate distribution[END_REF][START_REF] Stein | Estimation of the mean of a multivariate normal distribution[END_REF] in our NLM-SAP framework.

Proposition 1 (Stein) Let f (x) ∈ R, ε(x) ∼ N (0, σ 2 ), and Y(x) = f (x) + ε(x). Denote by γ(Y(x))
an estimate of f (x) and assume γ : R → R is absolutely continuous, and

i) lim |z|→∞ γ(f (x) + z)e -z 2 2σ 2 = 0, ii) E(γ(f (x) + ε(x))) 2 < +∞, and iii) E|γ ′ (f (x) + ε(x))| < +∞.
Then the following relation holds:

E|γ(Y(x)) -f (x)| 2 = E (γ(Y(x)) -Y(x)) 2 + 2σ 2 ∂γ(Y(x)) ∂ε(x) -σ 2 . ( 12
)
Proposition 1 provides an unbiased estimate of the risk for the K shape-based denoised values fk (x) at each pixel x. Indeed, if γ is one of the shape-based estimators and if the kernel ϕ is such that γ is absolutely continuous, then

r k (x) = ( fk (x) -Y(x)) 2 + 2σ ∂ fk (x) ∂ε(x) -σ 2 , (13)
is an unbiased estimate of the risk at pixel x for the k-th shape-based estimate, i.e., E(r

k (x)) = E| fk (x) -f (x)| 2 .
The main contribution of [START_REF] Van De Ville | SURE-based Non-Local Means[END_REF]16] is that they give a closed-form expression of ∂ fk (x) ∂ε(x) for NL-Means. Indeed, thanks to Eq. ( 2) and Eq. ( 4), its expression can be recast in the following form:

∂ fk (x) ∂ε(x) = ϕ(0) C x + 1 C x x ′ Y(x ′ ) ∂ω(x, x ′ ) ∂ε(x ′ ) - 1 C x x ′ Y(x ′ )ω(x, x ′ )
x ′′ ∂ω(x, x ′′ ) ∂ε(x) .

(

) 14
where

C x = x ′ ω(x, x ′ ) is a normalization constant.
In our NLM-SAP framework, our shape-based norm defined in Eq. ( 5) leads to the following expression of the derivative of the weights ω(x, x ′ ):

∂ω(x, x ′ ) ∂ε(x ′ ) = 1 h 2 S (x) Y(x) -Y(x ′ ) +S (2x -x ′ ) Y(x) -Y(2x -x ′ ) . ( 15
)
where S encodes the shape of our k-th shape-based estimator. Finally, combining equations ( 13), ( 14) and ( 15) leads to an unbiased estimate of the risk of our NLM-SAP denoiser. Now, since we are using K shapes, we get r 1 (x), • • • , r K (x), unbiased risk estimates respectively for the shape-based estimators f1 (x), • • • , fK (x).

Minimizing the risk of linear combinations

Suppose we want to solve the linear problem of aggregation. A natural way to find a good linear combination is to solve the following problem:

fLIN (x) := K k=1 α * k fk (x) , (16)
where α * := arg min

α∈R K E K k=1 α k fk (x) -f (x) 2 .
Here, the linearity of the combination allows to use Stein's Lemma (Proposition 1). It provides an unbiased estimate of the risk of the weighted average estimate:

E K k=1 α k fk (x) -f (x) 2 = E   K k=1 α k fk (x) -Y(x) 2 + 2σ 2 K k=1 α k ∂ fk (x) ∂ε(x)   -σ 2 .
Therefore, neglecting the term that does not depend on α in the last equation, the optimal weights α * can be obtained by solving the following:

α * (x) = arg min α∈R K k,k ′ α k α k ′ E( fk (x) fk ′ (x))+ K k=1 α k 2σ 2 E ∂ fk (x) ∂ε(x) -2E Y(x) fk (x) .
Using the first order optimality conditions, the last problem amounts to solving a linear system in α. This type of estimator is known to perform quite well for wavelet thresholding estimation where it is referred to as the SURE-LET (for SURE-Linear Expansion of Thresholds, see [START_REF] Blu | The SURE-LET approach to image denoising[END_REF]). However, in our framework, this method is not applicable. Indeed, since we work pixelwise, we have only one sample of the SURE for each pixel. This leads to an ill-conditioned system to solve, and the calculated α * behaves poorly. Algebraically, this is due to the fact that the matrix

[ fk (x) fk ′ (x)] k,k ′ =1,.
..,K is of rank one, and thus the system to solve is ill-conditioned. Our attempts to regularize the problem (e.g. with Tikhonov regularization known as "Ridge Regression" in statistics) or to solve the problem blockwise instead of pixelwise did not yield satisfying results. More generally, the problem of using SURE to take a local decision for each pixel x is difficult since this estimator has large oscillations (see for instance Fig. 5). In the next paragraph, we present how to regularize the risk maps, i.e., the "images" r 1 , • • • , r K .

Regularizing the risk maps

To make the risk estimates more robust, it is necessary to locally regularize the risk maps, in order to approximate at each pixel the expectations used in Eq. ( 13).

In [16], the convolution of the risk map is an efficient way to estimate the local risk in view of setting h since on both sides of an edge a large value of h should be used. Here, the anisotropy of the shapes implies that on one side of an edge the risk may be low whereas it may be very high on the other side. Since the convolution diffuses the risks on both sides of the edges, any comparison of the risks associated with each shape becomes unstable.

In order to diffuse the risks on each side of edges, we have adopted a heat equation with spatially and timely dependent coefficients (inspired by the Perona-Malik equation [START_REF] Perona | Scale space and edge detection using anisotropic diffusion[END_REF]). So, we let the risk maps r 1 , • • • , r K evolve according to:

   ∂r k ∂t (x, t) = div (g(|∇u(x, t)|)∇r k (x, t)) , r k (x, 0) = ( fk (x) -Y(x)) 2 + 2σ 2 ∂ fk (x) ∂ε(x) -σ 2 , ( 17
)
where g(x) = exp(-x 2 /κ 2 ), the parameter κ controls the anisotropy of the diffusion (the larger κ, the more isotropic the diffusion), and u is the smoothed noisy image which evolves at the same time using the Perona-Malik equation:

∂u ∂t (x, t) = div (g(|∇u(x, t)|)∇u(x, t)) , u(x, 0) = Y(x) . (18)
Curiously, we have noticed that we obtain better risk maps by diffusing √ r k instead of r k itself. Fig. 5 shows that this regularization procedure provides smooth risk maps, following edges of the underlying noise-free image, and finer than without regularization or with convolution.

Minimizer of the risk estimates (MRE)

With our measure of performance, the most natural way to address the selection problem (S) is to select the shape that minimizes the local risk estimate we have at hand:

fMRE (x) := fk * (x) where k * = arg min k r k (x) . ( 19
)
Such a selection rule is all the more relevant as the shapes are really different. In most cases, selecting the locally optimal shape yields satisfying results, but it is sometimes more appropriate to combine some of the best performing estimators as in the next paragraph.

Exponentially weighted aggregation (EWA)

In many cases, it might be better to combine several estimators rather than just selecting one. Especially, it happens to be relevant if the best estimators (in term of evaluated risk) are diversified enough or if the risk of the MRE was wrongly under-estimated. Thus, we have also used the statistical method of Exponentially Weighted Aggregation as introduced by Leung and Barron [START_REF] Leung | Information theory and mixing least-squares regressions[END_REF]. This method has been theoretically studied in [START_REF] Dalalyan | Aggregation by exponential weighting, sharp oracle inequalities and sparsity[END_REF][START_REF] Dalalyan | Aggregation by exponential weighting, sharp pac-bayesian bounds and sparsity[END_REF] and adapted for patch-based denoising in [START_REF] Salmon | NL-Means and aggregation procedures[END_REF]. It consists in aggregating the estimators by performing a weighted average with weights based on the confidence attributed to each estimator, measured in term of the risk. More precisely:

fEWA (x) := K k=1 α k fk (x) , (20)
with

α k = exp(-r k (x)/T ) K k ′ =1 exp(-r k ′ (x)/T )
.

The temperature parameter T > 0 is a smoothing parameter, that controls the confidence attributed to the risk estimates. If T → ∞, then the EWA is simply the uniform aggregate fUWA defined before. Conversely, when T → 0, then fEWA → fMRE . Most theoretical works about EWA (see [START_REF] Leung | Information theory and mixing least-squares regressions[END_REF][START_REF] Dalalyan | Aggregation by exponential weighting, sharp oracle inequalities and sparsity[END_REF][START_REF] Dalalyan | Aggregation by exponential weighting, sharp pac-bayesian bounds and sparsity[END_REF]) recommend a large temperature parameter T = 4σ 2 under few assumptions (like independence) on the estimators f1 (x), • • • , fK (x).

In practice, since assumptions on the estimators family may not be satisfied, we have used a smaller value, T = 0.02σ 2 .

Numerical and Visual Results

This section presents quantitative and qualitative results obtained on four images synthetically corrupted by AWGN. The corrupted images are obtained from four 256 × 256 noise-free images presented on Fig. 6: the famous cameraman image and city, windmill and lake1 . These images are particularly interesting in the study of our proposed Non-Local Means with Adaptive Patch Shapes (NLM-SAP) since they present high contrasted edges for which the classical NL-Means suffer from the rare patch effect (see Fig. 1).

Table 1 Gain by using multi-scale isotropic shapes in terms of PSNR and SSIM values (PSNR/SSIM). Circular patch shapes respectively with area of 12.5, 25 and 50 px 2 are used and their results are compared to the one obtained when using a combination of these three scales of patches. In all the experiments, unless otherwise specified, the NLM-SAP is used with the following default parameters:

search window: width ℓ = 11 px, -shape family: 15 shapes from families F1 and F2 on Fig. 2.a with shape areas of 12.5, 25 and 50 px 2 , -aggregation: EWA with T = 0.02σ 2 , and risk regularization: 50 iterations with timestep dt = 1/8 and κ = 30.

We set the bandwidth parameter to h = 0.7σ for patches of area 50 px 2 as suggested in [5]. However, as soon as we consider shapes of different areas, the parameter h has to adapt to the sizes of the shapes. To adapt the bandwidth to other sizes, we thus have extended a rule proposed by [START_REF] Kervrann | Optimal spatial adaptation for patch-based image denoising[END_REF][START_REF] Polzehl | Propagation-separation approach for local likelihood estimation[END_REF], based on quantiles of a χ 2 distribution, to the norm defined in Eq. ( 5). The quantile is selected such that h = 0.7σ for patches of area 50 px 2 . For the central pixel, we set its central weight as recommended in [START_REF] Salmon | On two parameters for denoising with Non-Local Means[END_REF]. The main limitation in computing time is due to the number K of shapes required by our NLM-SAP algorithm. We need to perform K times an NL-Means like algorithm, i.e., one for each shape. Thanks to our FFT acceleration, the computing time required for one shape, whatever the shape, is of about 2s for a 256 × 256 image with a Matlab implementation on an Intel Pentium 64-bit, 3.00 GHz. By comparison, the naive Matlab implementation of NL-Means takes about 100s, for square patches of area 7 × 7 px 2 . The computation of one local SURE map, using Perona-Malik diffusion, takes about 1s per shape. Finally, NLM-SAP 2 using 15 shapes leads to a computing time of about 45s which is less than the naive Matlab implementation of NL-Means.

Visual results are given to assess the denoising qualities relative to the different settings of NLM-SAP and to compare NLM-SAP with other denoising approaches (see Fig. 12). Numerical criteria support our claims: the 2 our Matlab implementation will be available on-line when this paper is published Peak Signal to Noise Ratio (PSNR) defined by

PSNR( f , f ) = 10 log 10 255 2 1 |Ω| x∈Ω (f (x) -f (x)) 2 , (21)
and the Structural SIMilarity (SSIM) defined in [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF].

The SSIM is between 0 and 1 and a value closed to 1 means that the estimated image has a similar structure to the noise-free image.

Behavior of NLM-SAP

In this section, we will study the behavior of NLM-SAP according to some parameters such as the type of family, the type of aggregation and the type of risk regularization. Each noisy image is corrupted by AWGN with standard deviation σ = 20. Table 1 and Figure 7 illustrate the gain of performance to use multi-scale patch shapes instead of using only one fixed size. In this experiment, we consider three circular shapes of areas: 12.5, 25 and 50 px 2 . Comparatively, for the original version of NL-means, the authors suggest to use square patches of fixed size 7 × 7 = 49 px 2 . Surprisingly, using the smallest shapes provides always the best PSNR. It means that the bias and the noise halos introduced by using large patches are actually more penalizing than the remaining noise left by the use of small patches. The aggregation of these three scales of shapes with our NLM-SAP methods improves Table 2 Gain by using anisotropic or mixture of isotropic and anisotropic shapes in terms of PSNR and SSIM values (PSNR/SSIM). The studied patch shapes are the isotropic disks, the half-pies, the quarter-pies, the bands and some combination of them (see Fig. 2 slightly the PSNR. Visually speaking, using only small isotropic patches already decreases the rare patch effect while using too large isotropic shapes produces a strong noise halo effect. However, the level of noise is much more decreased by using large shapes than small ones in homogeneous areas. Finally, combining different scales of isotropic patch shapes leads to a diminution of both the level of noise and the halo effect. This is well reflected in Tab. 1 by the gain in term of the SSIM criterion which provides quality measurements closer to our perception system. We will see in the following that the results can still be improved by considering both multi-scale and anisotropic patch shapes. ones presented on Fig. 2.a, i.e., the disks, the half-pies, the quarter-pies and the bands. Combination of these families are also studied. Our experiments show that most suitable shape families, both in terms of PSNR and SSIM, have to contain isotropic shapes, directional shapes and various scales of shapes. Increasing the number of shapes does not necessarily improve the denoising quality. Using 15 shapes from families F1 and F2 with the three different scales, seems to be a good trade-off between computing time and denoising quality. Figure 8 illustrates why using directional shapes is important to reduce the rare patch effect. Indeed, each oriented patch shape enables the restoration of edges in the target direction but is inappropriate in the other directions. Then a fine aggregation of them leads to high quality restoration of edges in all directions. Figure 9 displays weight maps induced by using patches with only one fixed square shape (i.e., NL-Means) compared to patches with adaptive scales and orientations (i.e., NLM-SAP). For NL-Means, all the weights are concentrated around the target pixel: this is the rare patch effect. For NLM-SAP, the weights are more spread, and other similar pixels are detected thanks to multi-scale and anisotropic patch shapes. It is clear that the limitation of the rare patch effect around high contrasted edges leads to a good reduction of the noise halo. Table 3 presents the numerical performances associated with four aggregation procedures: UWA, WAV, MRE and EWA. As expected, EWA provides best results, in terms of PSNR and SSIM, since compared to the other three it combines estimates with the best biasvariance trade-off. The local behaviors of NLM-SAP for WAV, MRE and EWA are presented on Fig. 10. The average areas and the average orientations of the selected shapes are given for the cameraman image. The selected patch shapes clearly adapt to the local scale and orientation of the image geometry. The chosen sizes of the shapes are smaller around edges and textured areas than in homogeneous areas. The chosen orientation follows the orthogonal orientation of the shapes which is consistent with the remarks given in Section 3. Compared to EWA, MRE suffers from brutal transitions, since it selects only one shape per pixel, while EWA evolves in a smoother way due to the weighted combination of shapes for each pixel.

Finally, we have studied the influence of the regularization of the risk maps on the aggregation results. Three methodologies are compared: aggregation using the noisy risk maps (i.e., SURE maps), the convoluted risk maps (using a disk kernel of radius 4) and the risk maps obtained by anisotropic diffusions. Table 4 gives the corresponding numerical performances and Fig. 11 illustrates the behavior of each type of risk map regularization. Surprisingly, the risk maps based on isotropic diffusions of SURE maps provide the best results in terms of PSNR and SSIM. However, the choice of the local sizes and orientations of the patch shapes is more relevant in the maps obtained by anisotropic diffusions, in terms of scale adaptivity, feature directions and spatial coherency. Using anisotropic diffusion, the NLM-SAP acts as expected by selecting big sizes of shapes, even around edges, since the shape orientations have been chosen properly to reduce the rare patch effect. By comparison, isotropic diffusion forces the size of shapes to be small around edges since it cannot select properly the suitable orientations. This difference of behaviors can be noticed around the cameraman's neck. In the case of the convolution, a small halo is present due to wrong decisions while for anisotropic diffusion this halo is strongly reduced. Then, even if isotropic diffusion provides the best numerical results in practice, anisotropic regularizations seems better theoretically founded to resolve this issue (see also Section 4.3). Other regularization strategies could then be investigated. To support this idea, Tab. 4 shows that there is still a gap of numerical performances between regularizations of the risk maps and the true risk maps (i.e., the local mean square errors provided by an oracle).

Comparisons with state-of-the-art methods

In this section, the proposed NLM-SAP approach is compared to state-of-the-art denoising methods. Comparisons have been performed with the classical (pixelwise) NL-Means [5], the blockwise NL-Means using UWA reprojection [5], the blockwise NL-Means using WAV reprojection [43], the pixelwise NL-means using SURE-based adaptive bandwidth selection [16], the Block-Matching and 3D filtering (BM3D) denoiser [7], and our proposed NLM-SAP approach.

Table 5 shows that NLM-SAP out-performs all other NL-Means improvements. NLM-SAP brings a gain of PSNR of about 1 dB compared to the classical NL-Means. The SSIM is also usually increased. The BM3D approach leads to better numerical results than all Non Local based approaches. Figure 12 gives the visual results. While the presence of the rare patch effect is very well illustrated by the noise halos for NLmeans and blockwise NL-Means using UWA aggregation, BM3D and NLM-SAP have reduced a lot this phenomenon. Our NLM-SAP provides smooth results with accurate details, such as the cameraman's head, the house windows, the windmill blades, the tree-trunk and the car. Visually, the quality of images obtained with NLM-SAP challenges those obtained with BM3D.

Conclusion

In this paper, we have addressed the problem of the rare patch effect arising in the NL-Means procedure and responsible of the noisy halos created around high contrasted edges. The proposed solution consists in substituting the square patches of fixed size by spatially adaptive patch shapes. A fast implementation of NL-Means, based on FFT calculations, has been proposed in this context to handle any kind of patch shape with arbitraty scale. Thanks to this acceleration, different estimates are obtained by using different patch shapes, typically one isotropic patch shape and four edge oriented patch shapes, all of them with three different scales. We have extended SURE-based approaches to aggregate properly these different shape-based estimates in a spatially adaptive way. To get an efficient locally adaptive filter, we have shown that the SURE-based risk maps require to be regularized and that anisotropic diffusion can be used to this purpose. Simulations have shown that exponentially weighted aggregation based on the regularized risk maps of the different shape-based estimates could lead to both numerical and visual improvements (the noise halo is suppressed around edges). Our method out-performs all the NL-Means improvements we have considered in our comparisons but is still outperformed by BM3D in terms of PSNR and SSIM. Future work could be to reduce computation time of the algorithm by choosing more convenient shapes or parallel implementations. Other method to regularize the risk maps could also be considered. Another extension might be to apply our framework to more general kinds of noise, by adapting results given by Hudson [START_REF] Hudson | A natural identity for exponential families with applications in multiparameter estimation[END_REF] on unbiased risk estimates for general exponential families. 12 Comparisons of the visual denoising performances of the proposed NLM-SAP approach and other state-of-the-art methodologies. From top to bottom, the input noisy images, the results obtained by the classical (pixelwise) NL-Means [5], the blockwise NL-Means using UWA reprojection [5], the BM3D denoiser [7], and our proposed NLM-SAP approach.

Fig. 1

 1 Fig.1Illustration of the noise halo appearing around high contrasted edges on images denoised by NL-means. The input noisy images was corrupted version of the noise-free images presented on Fig.6damaged by AWGN with standard deviation σ = 20. Noise halo arises from an abrupt lack of redundancy around edges sometimes referred as the rare patch effect.

Algorithm

  2D-FFT NL-Means for an arbitrary shape Inputs: noisy image Y, 2D-FFT of the shape F (S) Parameters: search window W , bandwidth h Output: estimated image f Initialize accumulator images A and B to zero for all shift vector δ in the search window W do Compute the square difference image ∆ δ ∆ δ (x) := (Y(x) -Y(x + δ)) 2 for all pixels x

Fig. 3

 3 Fig. 3 NL-Means pseudo-code for an arbitrary patch shape S. Pre-computations (based on 2D-FFT) of distances between shapes from the noisy image and shapes from its shifted version leads to a smaller complexity of O(|W | • |Ω| • log |Ω|), independent of the shape S.

Fig. 2

 2 Fig. 2 (a) Examples of shapes considered. The "central" pixel is shown in red, dark pixels illustrate high weights. Shapes are grouped in four categories: F1. the isotropic disk family, F2. the half-pies family, F3. the quarter-pies family and F4. the bands family. (b) Parametrization of the pie slices and bands.

Fig. 4

 4 Fig.4Examples of decentered patches near edges. If the patch is centered (a) fewer similar patches candidates are found than if the patch is decentered (c). The pixel of interest is in red, and black pixels are some similar pixels. (b) and (d) gives an illustration of the weights obtained (in black) by the two methods to denoise the red pixel.

Fig. 5

 5 Fig.5Maps of the estimated risk associated to patches with a small circular shape. From left to right, the underlying noise-free image, the map of the risk without regularization, with convolution (isotropic regularization) and with regularization based on anisotropic diffusion. Low risks are black, high ones are white.

Fig. 6

 6 Fig.6Chosen 256 × 256 noise-free images for our experiments, from left to right: cameraman, city, windmill and lake. These images present high contrasted edges for which the classical NL-Means suffer from the rare patch effect.

Fig. 7

 7 Fig. 7 Circular patch shapes: visual results obtained with NLM-SAP when using only circular shapes of different scales on a noisy realization of the cameraman image. The patch shape areas are respectively, from left to right, 12.5 px 2 , 50 px 2 and a combination of 12.5, 25 and 50 px 2 .

Fig. 8

 8 Fig. 8 Eight denoised images obtained for different oriented patch shapes. The proposed final aggregate is in the center.Each denoiser provides good performances in a specific target direction but suffers from noise halos in the other directions. The final (central) aggregate takes advantage of every oriented-denoiser to provide high quality restored edges. The patch shape used is indicated in white.

Fig. 9

 9 Fig. 9 Illustration of the diminution of the rare patch effect. (a) The noisy image with an highlighted target pixel and its neighborhood. (b) The associated maps of weights obtained by using only square patches of fixed size (i.e., NL-Means). (c) The associated maps of weights aggregating multi-scale and anisotropic patch shapes (i.e., NLM-SAP). Low weights are black, high ones are white.

Fig. 10 (

 10 Fig. 10 (top) Average areas and (bottom) average orientations of selected shapes for different aggregation procedures on a noisy realization of the cameraman image. From left to right: WAV, MRE and EWA aggregations. The average areas and the average orientations are represented using gray level colors whose legends are given on the top right corners.

Fig. 11 (

 11 Fig.11(top) Average areas and (middle line) average orientations of selected shapes for different risk maps on a noisy realization of the cameraman image. (bottom) Corresponding results focused on the cameraman's neck. From left to right, results using the noisy risk maps (i.e., SURE maps), the convoluted risk maps and the risk maps obtained by anisotropic diffusions. The average areas and the average orientations are represented using gray level colors whose legends are given on the top right corners.

  Fig.12Comparisons of the visual denoising performances of the proposed NLM-SAP approach and other state-of-the-art methodologies. From top to bottom, the input noisy images, the results obtained by the classical (pixelwise) NL-Means[5], the blockwise NL-Means using UWA reprojection[5], the BM3D denoiser[7], and our proposed NLM-SAP approach.

  .a).

			Cameraman	City	Windmill	Lake
	Noisy input image		22.13/0.400	22.13/0.567	22.13/0.385	22.13/0.456
	Patches with disk shapes	(family F1)	29.45/0.832	28.16/0.885	30.97/0.904	28.68/0.863
	Patches with half-pie shapes	(family F2)	29.43/0.832	28.08/0.886	30.97/0.906	28.60/0.863
	Patches with quarter-pie shapes	(family F3)	29.31/0.831	27.87/0.883	30.95/0.909	28.49/0.862
	Patches with band shapes	(family F4)	29.46/0.832	28.05/0.885	31.05/0.906	28.61/0.862
	Combination of shape families F1 and F2	29.50/0.833 28.21/0.887	31.11/0.907	28.73/0.865
	Combination of shape families F1, F2, F3 and F4	29.50/0.833	28.20/0.887	31.19/0.909	28.72/0.865

Table 2

 2 gives numerical results obtained by using different families. The compared families are the

Table 3

 3 Comparisons of different aggregation procedures in terms of PSNR and SSIM values (PSNR/SSIM). The compared aggregation types are UWA, WAV, MRE and EWA.

	Cameraman	City	Windmill	Lake

Table 4

 4 Comparisons of regularization procedures of the risk maps in terms of PSNR and SSIM values (PSNR/SSIM). The compared regularization procedures are the ones using the noisy risk maps directly (i.e., SURE maps), the convoluted risk maps and the risk maps obtained by anisotropic diffusions.

	Cameraman	City	Windmill	Lake

three sub-images extracted from the Laurent Condat's database: http://www.greyc.ensicaen.fr/ ~lcondat/

Acknowledgements The authors would like to thanks Jean-François Aujol, Loïc Denis, Yann Gousseau, Erwan Le Pennec and Florence Tupin for their interesting comments and criticisms.


Jean-Baptiste Charlier
email: jb.charlier@gmail.com

Catherine Bertrand

Stéphane Binet

Jacques Mudry

Nicolas Bouillier

Use of continuous measurements of dissolved organic matter fluorescence in groundwater to characterize fast infiltration through an unstable fractured hillslope (Valabres rockfall, French Alps)

Keywords:

Continuous measurements of natural fluorescence of dissolved organic matter (DOM) in groundwater have been used to characterize infiltration through an unstable fractured hillslope at the event time scale. Within the gneissic Valabres rockfall area in the southern French Alps, two sites, in the unstable massif and in the collapse area, were continuously monitored for electrical conductivity and for DOM natural fluorescence using a field fluorometer, and analyzed weekly for hydrochemical data. Two main results were found. From a methodological point of view, DOM natural fluorescence was a relevant tracer of fast infiltration in fissured media, knowing that continuous measurements improved the study of infiltration processes at the event time scale. From a hydrogeological point of view, the unstable fractured massif showed delayed dilution phenomena and the collapse areas showed fast and slow infiltration by piston-type flows via more or less open fractures. Consequently, from this conceptual model one may suppose that, unlike the collapse zones, the fractured hillslope threatening the valley is not submitted to strong pore-water pressure variations.

Introduction

Water infiltration is a key process in landslide triggering, and hydrogeochemical studies can give tools to characterize physical properties of rock slopes [START_REF] Bogaard | Hydrogeochemistry in landslide research: a review[END_REF]. In a fracture network within low porosity rock, a small infiltrated water volume can create large pressure. In this setting, fast infiltrations need to be well characterized in order to determine slope stability accurately.

Because infiltrated water contains large amounts of organic soil substances, the natural fluorescence (NF) of dissolved organic matter (DOM) and total organic carbon (TOC) are good tracers to discriminate water origin and flow paths in hydrosystems (e.g. [START_REF] Batiot | Use of total organic carbon (TOC) as tracer of diffuse infiltration in a dolomitic karstic system: the Nerja Cave (Andalusia, southern Spain)[END_REF][START_REF] Carstea | Continuous fluorescence assessment of organic matter variability on the Bournbrook River, Birmingham, UK[END_REF][START_REF] Meteo France | Percolation and particle transport in the unsaturated zone of a karst aquifer[END_REF]. For TOC, a non-conservative tracer, mineralization occurs at the scale of a few weeks; laboratory experiments showed TOC has an exponential decay with time [START_REF] Batiot | Étude expérimentale du cycle du carbone en régions karstiques[END_REF], meaning that the higher the TOC content, the faster the transit time in the system. TOC is thus adapted to study transfer time in hydrosystems at the event time scale.

However, the measurement of TOC is technically difficult due to continuous mineralization in sampling bottles and to the difficult automation of such devices in the field. To overcome these drawbacks, one may use a field fluorometer, developed for artificial dye tracers by [START_REF] Schnegg | A new field fluorometer for multi-tracer tests and turbidity measurement applied to hydrogeological problems[END_REF], which gives continuous measurements of water NF. Moreover, knowing that [START_REF] Cumberland | The freshwater dissolved organic matter fluorescence-total organic carbon relationship[END_REF] showed DOM fluorescence and TOC are linearly correlated in freshwater, it is hypothesized that continuous NF measurements could be used as a fine infiltration tracer. This article aims to characterize infiltration flows through an unstable basement massif. The original contribution of this work was to use continuous measurements of DOM NF in groundwater to monitor infiltration processes at the event time scale. The experimental site is the Valabres unstable hillslope in the southern French Alps. Two sites were monitored weekly for hydrochemical status and continuously for NF using a field fluorometer. The first site was a gallery draining the unsaturated zone of a fractured aquifer in the unstable massif, and the second was a spring draining a perched aquifer located in a highly fractured collapse zone. On the basis of contrasting results for these sites, a conceptual model is presented of infiltration flows at the event time scale.

Study area and sampling site

Study site

The Valabres hillslope is located in the Tinée valley in the southern Alps (Fig. 1a andb). The relief is marked by deep valleys with steep slopes, and the land cover is mainly composed of scattered shrubs and bare rock outcrops. In this area, the geographical relation between spring location and geology defines three main aquifer types (Fig. 1b): (1) a shallow aquifer in the colluvium (glacial deposits and screes), (2) an aquifer in the sedimentary rocks of Jurassic and Triassic ages, and (3) an aquifer in the gneissic rocks, whose springs are located near fractures or faults. The hillslope deformation is characterized by a continuous collapse area of about 2 km² (Fig. 1b) showing a double ridge oriented north-south that extends 2 km in the massif and bounded by counter-slopes. A decompressed area with open fractures lies between a vertical tectonic fracture towards the north-west and a fracture dipping about 45° towards the south-east. This collapsed area is bounded by undisturbed fractured rock forming stiff, currently unstable cliffs above the valley.

If one focuses on gneissic formations (with quartz, mica, plagioclase, K-feldspar, and amphibole), where main slope instabilities occur, most of springs are located at the middle of the slope, downstream of collapse zones (Fig. 1b). In the Valabres sector, the Tinée River crosses the regional geological structures. The bottom of the major trench zone is drained in the south-east by one of the most important springs of the valley, the Luicetta spring (Fig 1c).

The presence of a gallery, which enables performing investigations inside the rock mass, has served to monitor flows through the unstable massif. A sampling site was equipped in the gallery, and was coupled with the monitoring of the Luicetta spring, which drains the collapse area. These two sites in different geomorphological contexts were monitored physically and chemically. They are presented herein:

• Gallery. Flows from an undisturbed fracture in the gallery were collected to record the hydrological response of the unsaturated zone in the unstable fractured massif. The gallery has a geometry of 2.5 m width and of 2.3 m height, and is located under 100 m of overburden gneiss at the sampling site. According to [START_REF] Binet | Unstable rock slope hydrogeology: insights from the large-scale study of the western Argentera-Mercantour hillslopes (south-eastern France)[END_REF], the estimated infiltration transit velocity is between 1.10 -5 and 1.10 -4 m/s (response time of 10-20 days).

• Spring. Flows from the Luicetta stream at 750 m above sea level (ASL), downstream from the Luicetta spring (900 m ASL), are representative of flows in a shallow collapsed zone. The spring drains the decompressed area at the convergence of the main fractures (Fig. 2c), which is the boundary between a collapsed zone and the fractured gneiss. The catchment has an estimated area of 10 ha (+/-1 ha) and the spring has a minimum yield of 5 l/s. According to [START_REF] Binet | Unstable rock slope hydrogeology: insights from the large-scale study of the western Argentera-Mercantour hillslopes (south-eastern France)[END_REF], the estimated infiltration transit velocity ranges between 2.10 -2 and 8.10 - 3 m/s (response time of 1-4 days).

Method

Measurements were conducted from 28 October 2007 to 04 April 2008, that is, 160 days of monitoring.

Physical measurements

Daily rainfall intensity is an arithmetic mean of the time series of the weather stations of Rimplas and Saint-Etienne de Tinée (METEO France 2009), located 17 km to the south and 25 km to the north, respectively. Electrical conductivity (EC) and temperature measurements were carried out every 15 min using automatic field equipment (OTT, Orpheus, Aix-en-Provence, France).

DOM natural fluorescence (NF) was measured using a flow-through field fluorometer, GGUN-FL30 (Albillia, Neuchâtel, Switzerland), initially designed for tracer tests to measure dye tracer concentrations in the field [START_REF] Schnegg | An inexpensive field fluorometer for hydrogeological tracer tests with three tracers and turbidity measurement[END_REF][START_REF] Schnegg | A new field fluorometer for multi-tracer tests and turbidity measurement applied to hydrogeological problems[END_REF]. The field fluorometer is equipped with four optical components selected according to the absorption-emission spectra of the chosen artificial dye. Each component includes a light source (light emitting diode LED) exciting the fluorescence of the dye, and a photo-detector which yields an electric voltage of amplitude proportional to the fluorescence. The measured fluorescence by the photo-detector is filtered by an optical filter. Dye concentrations (in ppb) are determined by calibration with solutions of the dye at various concentrations. In this study, for DOM NF measurement, a 370 nm UV LED (suited for Tinopal artificial dye) was used for excitation, and the fluorescence was measured between 440 and 540 nm. The calibration of the field fluorometer was carried out using Tinopal at 10, 100, and 1,000 ppb. The detection limit for this tracer ranged between 0.14 and 0.2 ppb, and the error on the measurement range was estimated to be 8% [START_REF] Schnegg | An inexpensive field fluorometer for hydrogeological tracer tests with three tracers and turbidity measurement[END_REF]. Measurements were scheduled every 15 min. The light source of Tinopal was used because its excitation/emission wavelengths are in the range of organic matter fluorescence at high intensities: excitation between 200 and 400, and emission between 300 and 500 (i.e. Baker 2001; [START_REF] Katsuyama | Determining the sources of stormflow from the fluorescence properties of dissolved organic carbon in a forested headwater catchment[END_REF][START_REF] Blondel | Traçage spatial et temporel des eaux souterraines dans les hydrosystèmes karstiques par les matières organiques dissoutes[END_REF][START_REF] Carstea | Continuous fluorescence assessment of organic matter variability on the Bournbrook River, Birmingham, UK[END_REF]. In order to calibrate this probe, [START_REF] Savoy | Use of natural and artificial reactive tracers to investigate the transfer of solutes in karst systems[END_REF] has shown that a strong linear relationship between TOC and DOM NF (using a GGUN field fluorometer), exists in karstic systems. Moreover, [START_REF] Cumberland | The freshwater dissolved organic matter fluorescence-total organic carbon relationship[END_REF] showed that DOM fluorescence (using a laboratory spectrofluorometer) and TOC are linearly correlated in freshwater. From these last results, in various environments, it is hypothesised that NF was a consistent signal of DOM in surface and groundwaters in the presented gneiss formations.

Chemical measurements

Chemical measurements were carried out using an automatic sampler (Teledyne ISCO, Lincoln, NE) comprising 24 bottles previously rinsed with distilled water. Automatic samples were collected every 3 months (15 October 2007, 22 January 2008, and 11 March 2008). At these dates, these samples were completed by manual samplings to analyze HCO 3 -and to calculate the ionic balances, which were better than 8%, validating water analysis carried out by automatic samplers. Major ions were analyzed in the hydrogeological laboratory of Chrono-Environment (University of Franche-Comté, Besançon, France). For automatic samples, HCO 3 -was calculated in order to obtain a null ionic balance.

Results

Analysis of hydrogeological time series

The cumulative rainfall depth, P, for the whole monitoring period, was 374.4 mm. As shown in Fig. 2, this period was characterized by two main events: (1) event 1 with a P of 105 mm over the 5 days associated with 22 November 2007, for which there was a maximum rainfall intensity, Px, of 55 mm/day; (2) event 2 with a P of 131 mm over the 7 days associated with 11 January 2008, for which there was a Px of 62 mm/day.

The time series of EC in Fig. 2 showed opposite responses to rainfall events in the gallery and the spring sites. In the gallery, EC values were globally high, ranging from 391 to 429 μS/cm. The signal was stable, and after rainfall event 2, there was a slight decrease of 38 μS/cm attributed to a response time of 14 days and a transit time of 8.10 -5 m/s in the unsaturated zone of the massif. After that decrease, EC reached a stabilized plateau of about 404 μS/cm but never rose to the initial signal.

At the spring, the baseline of the EC signal in Fig. 2 ranged from 160 to 180 μS/cm, meaning that the residence time was lower than that of the gallery. One day after each rainfall event a drop was observed, followed by a recovery. For events 1 and 2, the decreases in EC were 18 and 46 μS/cm, respectively. After both events, the recovery of the EC signal formed a slight dome (10 days after rainfall peak) of more mineralized water; the maximum EC dome was higher than in the pre-event conditions.

As for EC, the NF time series plotted in Fig. 2 showed opposite variations for the two sites. For the gallery, there was a quasi-null NF signal (0-0.2 ppb). At the spring, there was a baseline between 0.25 and 0.5 ppb and peaks of NF reaching 8.5 and 12.0 ppb after rainfall events 1 and 2, respectively. The response time was 1 day.

Analysis of hydrochemical time series

The results of hydrochemical analysis are synthesized in Table 1. The two main observations were that (1) the highest concentrations were for HCO 3 -, Ca 2+ , and SO 4 2-, which were thus plotted in Fig. 2 and (2) as expected by the EC analysis in section Analysis of hydrogeological time series, concentrations in the gallery were 2-to 4-fold higher than for the spring. First, for the gallery, variations of Ca 2+ , SO 4 2-, and HCO 3 -were similar. Before event 2, Ca 2+ progressively rose from 55.0 to 72.5 mg/l, SO 4 2-from 70.0 to 74.5 mg/l, and HCO 3 -from 94.0 to 147.2 mg/l followed after this event by a decrease of 15.5, 8.6, and 40.5 mg/l, respectively. This pattern showed there was durable dilution of groundwater flows, as observed for the EC signal. Second, for the spring, Ca 2+ and HCO 3 -variations were different from the variation of SO 4 2-. Ca 2+ and HCO 3 -showed a rise of 6.5 mg/l from 21.5 to 28.0 mg/l, and of 22.3 mg/l from 51.1 to 73.4 mg/l, respectively, which was related to the EC dome 10 days after rainfall peak, whereas SO 4 2-was very stable with a mean value of 18.5 +/-0.5 mg/l.

Discussion on fast and diffuse infiltrations

The analysis of EC and NF fluctuations, shown in Fig. 2, revealed the dynamics differed for the two sites: a buffered response to rainfall (flat chemograph) in the gallery and a reaction with high intensity for the spring. First, for the gallery, the durable dilution of EC by event rainwater means that the system is highly inertial, due for example to a low permeability. The quasi-null NF signal showed that a pluri-week transit enabled mineralization of all organic substances in water. These observations showed that infiltration flows were diffuse.

Second, for the spring, EC and NF baselines indicated a shorter residence time than in the gallery, with higher NF and lower EC (i.e. lower mineralization) at the spring. To better characterize infiltration processes, NF vs. EC was plotted in Fig. 3. A clockwise hysteresis loop for each rainfall event was observed, shown by arrows. These hysteresis loops reflected differences in water-transit velocity. Three phases can be distinguished: (phase I) a decrease of EC (arrival of low mineralized water) without any change in DOM concentration; (phase II) an enrichment in organic substances (DOM peak) of less mineralized water, with a response time inferior to 1 day; (phase III) a higher increase of mineralization than the preevent water, and a decrease of DOM concentration, with a response time of 10 days, ending the hysteresis. These phases are coherent with the three following infiltration processes: (1) the occurrence of non-infiltrated event water (neither highly mineralized nor enriched with soil organic matter) by surface runoff, (2) the occurrence of infiltrated water by fast piston flows (not mineralized yet but enriched with soil organic matter), and (3) the occurrence of infiltrated water with a longer transit time allowing its mineralization (more mineralized but less enriched in soil organic matter).

The results are summarized in an infiltration model of the Valabres unstable rocks. Presented in Fig. 4, this model differentiates infiltration flows through the unsaturated zone of the deep aquifer localized in the unstable fractured massif, as well as in the perched aquifer localized in the weathered collapsed area. First, in the unsaturated zone of the unstable fractured medium (gallery), the durable decrease of chemographs and EC signal show mixing and dilution processes, as well as low storage. The small fluctuations of these signals also show a buffering effect of the massif characterized by low velocities of about 8.10 -5 m/s as stated by [START_REF] Binet | Unstable rock slope hydrogeology: insights from the large-scale study of the western Argentera-Mercantour hillslopes (south-eastern France)[END_REF] for the same site. Second, in the perched aquifer of the collapsed area (spring), three types of organized flows were distinguished: (1) a fast transit of surface runoff upstream from the Luicetta outlet (response time <24 h), (2) a fast transit via fast piston effect of water infiltrated in permeable fractures (transit velocity of about 10 -2 m/s) with low water storativity, and (3) a 10-fold slower transit of water enriched in Ca 2+ in less permeable media with higher storativity.

To study flows in fractured media, the main advantage of using continuous NF measurements is to expand the accuracy of EC when very recent seepage water is observed (about 1 week). In one case, EC signal at the outlet of an unstable rock is the result of two features for infiltrated water: its initial mineralization (rainwater), and its mineralization occurring during the transit time, with possibilities of water-rock interactions. In the other case, NF signal is the result of infiltration through the soil (enrichment of DOM) and is limited in time, because the mineralization of organic substances generates a rapid decrease of NF concentration [START_REF] Batiot | Étude expérimentale du cycle du carbone en régions karstiques[END_REF]. Thus, NF may be an irreplaceable supplement of EC because NF is a relevant tracer of fast infiltration, whereas EC is an integrating tracer of several concomitant processes as fast infiltration, or not, and residence time. This method serves well to discriminate flowpaths and appears to be widely adapted for fractured and karstic media.

The infiltration model developed using a hydrochemical approach may serve to better understand repartitions and variations of pore-water pressures in unstable hillslopes. In fact, due to the contrasted infiltration processes in the various fracture networks highlighted in the model, one may suppose that pore-water pressures, which increase during the infiltration events, are dissipated by the collapsed zone and that the fractured hillslope threatening the valley is not submitted to strong pore-water pressure variations. This study calls for further investigations in order to improve hydromechanical models developed for unstable hillslopes, including the heterogeneity of hydrogeological responses highlighted in such a context.

Conclusion

The aim of this study was to use continuous measurements of natural fluorescence (NF) of DOM in groundwater to characterize infiltration at two sites of an unstable hillslope: a fractured medium and a collapsed zone. The results showed first that, from a methodological point of view, a continuous NF signal is a relevant tracer of fast infiltration that improves the study of flow processes at the event time scale. Second, from a hydrogeological point of view, the unstable fractured massif showed delayed dilution phenomena and the collapse areas showed fast and slow infiltration with piston flow types via more or less open fractures. Finally, the infiltration model designed through this hydrochemical approach gives tools to better characterize physical properties in an unstable hillslope.

Figures and Tables

Fig. 1 a

 1 Fig. 1 a Location of the Valabres hillslope in the southern French Alps. b Hydrogeological map of the study site indicating the location of the gallery and the Luicetta spring. c Deep cross-section A-B showing the geomorphology of the Valabres hillslope

Table 1

 1 Temporal variability of major elements measured weekly from 22 October 2007 to 04 April 2008

	Ca 2+ Mg 2+ Na + K + F -Cl -NO 3	-SO 4	2-HCO 3 a	-
	[mg/l]
	Gallery
	min 53.50 3.20 4.40 0.93 0.01 4.95 2.47 65.33 92.23
	mean 59.70 3.45 4.69 1.03 0.08 5.87 3.80 70.18 110.41
	max 72.50 3.80 5.20 1.24 0.10 6.90 5.39 75.15 147.21
	SD 5.59 0.17 0.24 0.09 0.02 0.63 1.33 3.51 15.43
	Spring
	min 21.00 2.40 2.20 0.77 0.02 4.01 0.00 17.28 47.95
	mean 23.67 2.75 2.60 0.81 0.03 5.06 2.06 18.46 60.03
	max 28.00 3.20 3.00 0.92 0.05 6.53 3.89 18.96 73.41
	SD 2.29 0.20 0.21 0.04 0.01 0.60 1.12 0.51 7.92

a HCO 3 -calculated to obtain a null ionic balance

Acknowledgements This research was financed by INERIS (Institut National de l'Environnement Industriel et des Risques) and the French Ministry of Ecology and Sustainable Development, through STABROCK "Impact of climate changes on the stability of rock slopes: observation, experimentation and modeling". Many thanks to Bruno Regent for the field work.