Antoine Houdard
email: antoine.houdard@telecom-paristech.fr

Andrès Almansa
email: andres.almansa@parisdescartes.fr

Julie Delon
email: delon@parisdescartes.fr

Andrés Almansa

A Almansa

Demystifying the asymptotic behavior of global denoising

Keywords: Diagonal Estimation,, Global Denoising,, Wiener Filtering,, Asymptotic study

published or not. The documents may come

Introduction

Image denoising usually amounts to estimate an image u from the observation

u = u + , (1)
where is an additive noise. Classically, is assumed to be the realization of a random vector following a Gaussian distribution N (0, σ 2 I N ). This noise model is widely used since the Anscombe transform permits to transform the more realistic Poisson noise in a nearly Gaussian noise with fixed variance. The story of image denoising, and more generally of image restoration, is probably as old as the story of image processing. If classical image denoising methods have always gathered various mathematical tools, such as neighborhood filters [START_REF] Yaroslavsky | Digital picture processing: an introduction[END_REF], variational models [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], non linear Partial Differential Equations [START_REF] Guichard | A review of pde models in image processing and image analysis[END_REF] or transform domain estimation [16], the last true revolution came with the introduction of patch-based methods in 2004 [START_REF] Awate | Image denoising with unsupervised information-theoretic adaptive filtering[END_REF][START_REF] Buades | A non-local algorithm for image denoising[END_REF]. At the time, trendy restoration approaches (total variation or wavelet thresholding) were not able to provide a satisfying trade off between texture preservation, flat areas restoration, detail reconstruction, and apparition of oscillating artifacts. By exploiting the selfsimilarity in images, patch-based approaches permitted to significantly improve the performance of image denoising algorithms in this regard.

Leading methods nowadays are all patch-based [START_REF] Dabov | Image denoising by sparse 3-d transform-domain collaborative filtering[END_REF][START_REF] Lebrun | A Nonlocal Bayesian Image Denoising Algorithm[END_REF][START_REF] Yu | Solving inverse problems with piecewise linear estimators: from Gaussian mixture models to structured sparsity[END_REF][START_REF] Pierazzo | Da3d: Fast and data adaptive dual domain denoising[END_REF]. These extremely popular approaches have been adopted in a huge range of applications. Their underlying assumption being that similar patches can be seen as independent realizations of the same distribution, the performance of a denoising algorithm should increase when the number of realization increases. Theoretically, this should lead to a form of asymptotic opti-mality when the image size tends toward infinity. Consistency results, under stationnarity hypotheses, have been shown for instance for the DUDE algorithm [START_REF] Ordentlich | A discrete universal denoiser and its application to binary images[END_REF][START_REF] Weissman | Universal discrete denoising: Known channel[END_REF] and for the Non Local Means [START_REF] Buades | A non-local algorithm for image denoising[END_REF]. Now, despite their non-local nature, most of these algorithms limit the search area for similar patches to a medium-sized neighborhood around each pixel. Doing otherwise would confront them to a dilemma [START_REF] Duval | A bias-variance approach for the nonlocal means[END_REF]. A larger search size means potentially more similar patches, reducing the variance of the denoising estimator. However increasing the search area in natural images also tends to increase the risk to consider dissimilar patches as similar, thus increasing the bias of the denoising estimator. Most authors found the best compromise in relatively small search areas. As a consequence, increasing the image size does not necessarily improve denoising performance. This observation was supported by extensive experimentation in [START_REF] Levin | Patch complexity, finite pixel correlations and optimal denoising[END_REF], who showed that even if an infinite database of natural image examples was available, non-local denoising performance would attain an asymptotic performance that does not tend to infinite signal to noise ratio. Non-local methods seemed to be doomed to fundamental limits that could not be overcome.

In 2012, Talebi and Milanfar [START_REF] Talebi | Global image denoising. Image Processing[END_REF][START_REF] Talebi | Global denoising is asymptotically optimal[END_REF] proposed a truly global denoising approach where each pixel is used to denoise every other pixel. They claimed in a subsequent paper [START_REF] Talebi | Asymptotic performance of global denoising[END_REF] that this approach is asymptotically optimal, in the sense that "the mean-squared error monotonically decays with increasing image size", regardless of image content, at least in an oracle scenario. In this context, this paper raises again the question: can denoising methods be fixed in such a way that they attain infinite PSNR when given an infinite number of examples (or an image of infinite size) ? They opened the debate by showing that given an oracle, such an asymptotic performance seems to be possible. However two questions are still left open:

1. What conditions has to satisfy an infinite image for the asymptotic result to hold? 2. Do these conclusions extend to the non-oracle case?

This article tries to give a precise answer to the first question, and some elements of response to the second one. To do so, we revisit the theory of diagonal estimation (refered to as Wiener filtering in Talebi's paper) that was first developed for wavelet bases. Considering images as vectors of R N , a diagonal estimator u = W u is a non linear estimator of u that is diagonal in a given orthonormal basis V = {V i } i=1,...,N , which means that it can be written

u = W u = V ΛV T u = N k=1 λ k ( u) • u, V k • V k , ( 2
)
where Λ is a diagonal matrix whose k th coefficient λ k ( u) depends on the observation u (otherwise, the resulting estimator would be linear). The diagonal estimation framework is widely used in image processing: the basis V is often chosen as a Fourier or a wavelet basis [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF][START_REF] Donoho | Ideal denoising in an orthonormal basis chosen from a library of bases[END_REF]16], or can for instance be built up as an orthonormal dictionnary from the image itself [START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF]. The success of diagonal estimation stems partly from the fact that if the image u is sparse in the orthonormal basis V , these "diagonal estimators are nearly optimal among all non linear estimators", as stated in [16]. The global denoising formalism introduced by Talebi and Milanfar [START_REF] Talebi | Global image denoising. Image Processing[END_REF][START_REF] Talebi | Global denoising is asymptotically optimal[END_REF] can be reinterpreted in this context. Indeed, the idea of global denoising boils down to build V as an orthonormal basis that diagonalizes a given denoising filter (such as NLmeans [START_REF] Buades | A non-local algorithm for image denoising[END_REF]) computed on u. In the context of diagonal estimation, we will derive a novel asymptotic study of global denoising. Basically, we introduce precise conditions both on the image and the global filter to ensure that the mean-squared error

MSE( u|u) def = 1 N E( u -u 2 ), (3)
for global image denoising decays toward zero for increasing image size. We will see that classical results of the diagonal estimation theory also permit to envision possible answers to the question of the extension of global denoising to the non oracle case. The paper is organized as follows. In Section 2 we provide a short reminder on the theory of diagonal estimation in an orthonormal basis. The first contribution of this paper is to revisit this framework to present the global denoising formalism and to put it into perspective relatively to classical diagonal estimation results. The second and main contribution is the novel asymptotic study of global denoising presented in Section 3. Finally, in Section 4, we also discuss and show experiments on several open issues, including the extension of these results to the non-oracle case.

Global filtering revisited

Diagonal estimation : a short reminder

We recall here the basic properties of a diagonal estimator in terms of quadratic risk minimization, before revisiting the theory of global denoising in this context.

Quadratic risk

Assume that W is deterministic, i.e. that the coefficients λ k are independent of the random noise and only rely on the unknown image u. In this case, the mean quadratic risk or mean squared error (MSE) of the diagonal estimator given by Equation ( 2) can be easily derived. Let us denote by b the projection of the unknown image u in the orthogonal basis V , that is b = V T u. The MSE between u and u can be written as a function of the eigenvalues (λ k ) and the projection b, as a sum of a variance and bias terms.

Proposition 1 Let u = V ΛV T u, with V ΛV T a deter- ministic filter. Then, MSE( u|u) = 1 N N j=1 (1 -λ j ) 2 b 2 j + σ 2 λ 2 j . (4)
Proof (Proof of Proposition 1)

N • MSE( u|u) def = E( u -u 2 ) = E( V Λ T u -V ΛV T u 2 ) variance term + E( V ΛV T u -u 2 ) bias term = E( V ΛV T 2 ) + E( (Λ -I N )V T u 2 ) = ( N i=1 λ 2 i )σ 2 + N i=1 (λ i -1) 2 E[(V T u) 2 i ].
Observe that the last equality holds only because the filter W = V ΛV T does not depend on the noise .

Oracle quadratic risk minimization

Minimization of the MSE w.r.t. the {λ i }'s For a fixed orthonormal basis V , the previous MSE is a convex function of the eigenvalues λ i , and reaches its global minimum for

λ i = b 2 i σ 2 + b 2 i . (5)
The corresponding minimal value of the MSE is

MSE := MSE(λ ) = σ 2 N N j=1 λ j = σ 2 N N j=1 b 2 i σ 2 + b 2 i . (6)
This formula shares similarities with Wiener filters, with the difference that the coordinates {b i } are not expected values but actually depends on the oracle image u, which is assumed to be deterministic. This oracle MSE cannot be attained in practice but only represents a lower bound for the quadratic risk of diagonal estimators in the basis V . However, it can be shown that some well chosen thresholding estimators have a risk which is not too far from the oracle one [16].

Minimization w.r.t. the {b i }'s The previous oracle diagonal estimation is done in a given basis V , which could for instance be chosen as a Discrete Cosine Transform basis or a Wavelet basis. Obviously, the final estimation strongly depends on this choice, and one might wonder in practice how to optimize the selection of the basis V for a given image u. The quantity MSE from equation ( 6

Ψ (b, µ) = σ 2 N N j=1 b 2 i σ 2 + b 2 i -µ   N i=j b 2 j -u 2 2   ,
where µ is a Lagrange multiplier. The derivation with respect to the b i yields

∂ bi ψ(b, µ) = 2σ 4 N b i (σ 2 + b 2 i ) 2 -2µb i , (7)
and the derivation w.r.t. µ

∂ µ ψ(b, µ) = u 2 2 - N j=1 b 2 j . (8)
Setting [START_REF] Donoho | Ideal denoising in an orthonormal basis chosen from a library of bases[END_REF] to zero implies the existence of i 0 such that b i0 = 0 (otherwise u 2 2 would be zero). Then, setting [START_REF] Dabov | Image denoising by sparse 3-d transform-domain collaborative filtering[END_REF] to zero yields for i 0

µ = σ 4 N 1 (σ 2 + b 2 i0 ) 2 ,
and for each i = i 0 b i σ 4 N (σ 2 + b 2 i0 ) 2 -(σ 2 + b 2 i ) 2 (σ 2 + b 2 i0 ) 2 (σ 2 + b 2 i ) 2 = 0, (9)
which implies b i = 0 or b 2 i = b 2 i0 . Using (8) again gives the following generic form for the critical points of

Ψ b i =      ± u 2 2 #I if i ∈ I 0 otherwise,
where I ∈ P ({1, . . . , N }) \ ∅ is the support of b.

Let b be a critical point and I its support, we have

MSE (b) = σ 2 N u 2 2 σ 2 + u 2 2 #I σ 2 N u 2 2 σ 2 + u 2 2
where the equality occurs if and only if #I = 1. Thus, among all critical points, the minimal ones are the b = ± u 2 e i . Finally, they are also global minima for Ψ as we have for all b

MSE (b) σ 2 N b 2 j σ 2 + u 2 2 = σ 2 N u 2 2 σ 2 + u 2 2 = MSE (b ).
The previous Proposition means that an optimal basis V should be such that V T u = u 2 e i , for a given i in {1, . . . , N }. It follows that V should be composed of the vector u u 2 and simply completed in an orthonormal basis. The resulting oracle filter would be

(W ) ij = u i u j σ 2 + u 2 2 . ( 10
)
Again, even if this filter is not reachable since it depends on the unknown oracle u, this results strongly support the intuitive idea that ideal bases should provide a sparse representation of u. In practice, diagonal estimation should be applied in a well-adapted basis for each image, typically a basis V that provides a fast decrease of the {b j }. The principle of global filtering [START_REF] Talebi | Global image denoising. Image Processing[END_REF], described in Section 2, is to rely on classical non linear filters from the denoising literature to choose V .

Non-oracle case

The oracle u and its projection b = V T u being unavailable, we need a way to approximate the previous estimation from the knowledge of b = V T u. A classical solution is to consider (hard or soft) thresholding estimators in a given orthonormal basis V , in order to discard irrelevant components. As illustrated by Figure 1, the result of hard thresholding is far from being as satisfying as the oracle estimation (5), at least on a Discrete Cosine Basis. However, for a specific value of the threshold T , the mean squared error obtained with a hard or soft thresholding can still be controlled by the one of the oracle attenuation (5). Proofs of this theorem can be found in [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF] or [16]. It helps to predict what kind of images can be well denoised by hard thresholding in a given basis. For a DCT basis for instance, we can expect a lower oracle MSE(λ ) for smoother images, and the same property should hold for thresholding. This is illustrated by Figure 1, which shows two noisy images and their respective denoised versions by oracle attenuation and hard thresholding. At the same noise level, the second image has a better oracle result (d) than the first one, and this is also true for the hard thresholding result (c). We can also conclude that a basis V nearly optimal for the oracle should also be a good choice for the thresholding estimation. Observe that the threshold T = σ √ 2 ln N is not really optimal in practice. A good way to fix T for soft thresholding is to resort to the SURE estimator of the MSE [16].

Global filtering in this context

Global denoising [START_REF] Talebi | Global image denoising. Image Processing[END_REF][START_REF] Talebi | Global denoising is asymptotically optimal[END_REF] draws on the concept of diagonal estimation in order to improve current denoising filters. As described in the previous section, a well chosen basis should provide a sparse representation of u and a general basis obviously cannot fit well for all natural images. Global denoising builds V as an orthonormal basis that diagonalizes a classical non linear denoising filter (such as NLmeans [START_REF] Buades | A non-local algorithm for image denoising[END_REF]), computed on u. The underlying assumption is that if the chosen denoising filter is well adapted to the image, the coefficients b j will decrease relatively quickly and the diagonal estimate will be all the more efficient.

Principle of global denoising

Assuming the same image formation model (1), numerous classic denoising filters, such as Gaussian or bilateral filters as well as NL-means [START_REF] Buades | A non-local algorithm for image denoising[END_REF] type filters, can be written under the form

u = W u, (11)
where W = D -1 K, with K the positive definite kernel from the filter and D a diagonal matrix with entries

D ii = j K ij , i ∈ {1, . . . , N } 1 .
Starting from a given denoising filter W , the idea of global denoising, made popular by Milanfar in [START_REF] Milanfar | Symmetrizing smoothing filters[END_REF], [START_REF] Milanfar | A tour of modern image filtering: New insights and methods, both practical and theoretical[END_REF], is to modify this filter W , in order to decrease the mean square error between u and u. For instance, if we assume that W can be diagonalized in an orthonormal basis V (this can be ensured by symmetrizing it, as described in the next section), the oracle attenuation (5) of the eigenvalues can be applied to improve the filter.

Symmetrizing the filter W

To ensure the fact that W can be diagonalized in an orthonormal basis, the authors of [START_REF] Talebi | Global image denoising. Image Processing[END_REF] propose to replace 1 For instance, for NL-means we would have

K i,j = e - P i -P j 2 2h 2
, with P i and P j the patches centered at i and j and h a parameter this filter by a symmetric doubly stochastic version W s of W which minimizes the cross-entropy

i,j W s ij log W s ij W ij . (12)
In practice, this minimization problem can be solved numerically with the Sinkhorn algorithm, which consists in iteratively normalizing the rows and the columns of W until convergence. Starting from a positive definite kernel K, it can be shown that the resulting filter W s is positive definite, symmetric and doubly stochastic, and that its eigenvalues are very close to those of W [START_REF] Milanfar | Symmetrizing smoothing filters[END_REF]. In practice, the denoising results obtained with this symmetrized filter appear to be equivalent or slightly better than the ones obtained with W [START_REF] Milanfar | Symmetrizing smoothing filters[END_REF][START_REF] Chan | Demystifying symmetric smoothing filters[END_REF]. In the following, we always consider the filter W in its symmetric and doubly stochastic version.

Deterministic filter

The mean-squared error formulation ( 4) is valid only if the filter W = V ΛV T is deterministic. This assumption is sensible when V is fixed, as a DCT or wavelet basis for instance. However when the filter W comes from the noisy image u, this hypothesis does not hold.

To illustrate this fact, we compare, for different choices of the filter W , the theoretical MSE theo computed by formula (4) with the experimental mean-squared error

MSE eval = 1 N N j=1 |u j -W u j | 2 .
Figure 2 shows the relative error

| MSE eval -MSE theo | MSE theo (13)
for the following filters W , computed on three different images:

1 The first filter is independent from the noise present in u, so the relative error is very small. On the contrary, the NL-means filter computed directly on the noisy image strongly depends on the noise realization, and the relative error between the theoretical and experimental MSE remains above 10% for all three images. Finally, observe that if the NLM is computed on a version of u that has already been denoised in a first step (by a NLM kernel or another denoising procedure), the resulting W seems to be partly decoupled from the noise, at least enough for the theoretical MSE theo to be a good predictor of MSE eval .

(a) (b) (c) images (a) (b) (c)
O-NLM 6.6 % (± 3.0) 0.3 % (± 1.9) 0.4 % (± 0.9) NLM 34.2 % (± 3.0) 24.6 % (± 1.9) 11.2 % (± 0.9) P-NLM 6.8 % (± 1.7) 2.7 % (± 0.3) 1.0 % (± 0.3)

Two oracle levels

In the previous section we used the noiseless image u as an oracle to compute the weights of the O-NLM filter W . Note that this use of the oracle is different from the one introduced in equation ( 5) (section 2.1.2) to compute the optimal eigenvalues λ j . In general it will be clear from the context which level of oracle we refer to. In ambiguous cases we shall refer to the first one as W -oracle and to the second one as λ-oracle.

Discussion

By producing a basis V that is well-adapted to the image we want to denoise, global image denoising usually produces better results than a diagonal estimation on a DCT or wavelet basis. However, global denoising still suffers from two major issues:

first, the λ-oracle u is needed in order to optimize the eigenvalues; second, memory cost and computation time are untractable because of the eigendecomposition of the filter W of size N × N .

In order to bypass the first issue, we saw in section 2.1.3 that hard or soft thresholding could provide MSE re-sults controlled by the optimal MSE * . Another possibility would be to try multiple sets of eigenvalues and keep the ones minimizing a SURE estimator of the MSE. This is the solution proposed by the GLIDE algorithm [START_REF] Talebi | Global denoising is asymptotically optimal[END_REF]. The second issue can be solved by computing only a small percentage of eigenvectors. In GLIDE, Talebi and Milanfar make use of the Nyström extension in order to approximate the filter W and its first eigenvalues.

Asymptotic study

In this part, we study the asymptotic behavior of the MSE given by formula (4) when the image size increases. In [START_REF] Talebi | Asymptotic performance of global denoising[END_REF], the authors claim that global denoising is asymptotically optimal, in the sense that the MSE in (3) tends to zero. Before going further, let us mention that this decay of the global MSE may occur while some local areas of the image remain poorly denoised even when the image size tends to infinity, as it is shown on Figure 3. In order to explore the precise conditions of this convergence, we define in Section 3.2 a reasonable model for an image whose size grows to infinity. We also assume a parametric model for the decay of the coefficients b j , and we derive in Section 3.3 different conditions of convergence for the MSE and its corresponding decay rate. Finally, in Section 4.2, we discuss and illustrate these different results and the realism of these models for different choices of images and filters W .

In the following, we always consider that the filter W is independent from the noise .

Upper bound on the optimal MSE

We have seen in Section 2.1 that the oracle risk for diagonal estimation was given by

MSE = σ 2 N N j=1 b 2 j σ 2 + b 2 j
, with b = V T u the projection of the oracle image u in the eigenbasis V . Now, this MSE can be upper bounded by the l 1 -norm of b divided by N : This shows that even with a global MSE converging to zero, the restoration can remain locally bad and even get worse when the image size increases.

MSE = σ 2 N N j=1 b 2 j σ 2 + b 2 j σ 2 N N j=1 b 2 j 2σ|b j | = σ 2N b 1 . (14)
The authors of [START_REF] Talebi | Asymptotic performance of global denoising[END_REF] suggest that this upper bound might converge towards 0 when N grows to infinity. In order to prove this convergence, they assume that the sorted coefficients |b j | drop off at a given rate α > 0

|b j | C j α .
We shall see below and in section 3.2 that this models requires C to depend on N to make sense. Nevertheless for a fixed image size, this hypothesis seems quite reasonable for different existing filters, as illustrated by Figure 4. When working with Fourier or space-frequency decompositions, the value of α was shown to be related to the regularity of the image [16], and values of α be- tween 0.5 and 1 were shown to be in agreement to actual image data [START_REF] Facciolo | Irregular to Regular Sampling, Denoising, and Deconvolution[END_REF]. Such models have also been used in asymptotic studies where the image resolution tends to infinity, but here we are interested in the asymptotic behavior when image size grows to infinity at constant resolution. In this particular kind of asymptotic study, we cannot expect the rate α and the constant C to remain constant when the image size grows towards infinity. Put another way, there is no reason that we can bound the b j coefficients independently of the image size N . To demonstrate this claim, we propose a model for an image whose size grows to infinity. Under this model we will show that the coefficients b N j actually depend on N. Then we propose a more complete parametric model for the coefficients decay

c 2 #{u N i = 0} N u N 2 2 N -→ N →∞ 0. ( 16
)
This shows that the ratio of non-zero pixels collapses when the image size goes to infinity. This case will not be considered in the following. Indeed if we consider an infinite image U such that u N 2 = o( √ N ) then, the upper-bound ( 14) tends to zero when N goes to infinity:

b N 1 √ N b N 2 = √ N u N 2 = o(N ),
which implies the convergence. We provide in Section 4 an experiment with an image padded with zeros illustrating this case. This leads us to define the widespread infinite image model as follows.

Hypothesis 1 (Widespread infinite image model) Let U be an infinite image and denote u N its truncation of size N . Then U is said to be non sparse if there exists m > 0 and M > 0 such that m

√ N u N 2 M √ N . ( 17
)
Domination decay model Now consider a sequence of orthogonal bases V N (the eigenbases of symmetric filtering operators W N ). Recall that we denote by b N = V N u N the projection of the image of size N on the corresponding eigenbasis. We need a realistic model on the asymptotic behaviour of b N j when N, j go to infinity. In this part we design an upper bound for |b N j | which is both simple and easy to manipulate to prove convergence results; adapted to the data, in the sense it has the same shape as the |b N j |. In order to design it we start with the model from [START_REF] Talebi | Global denoising is asymptotically optimal[END_REF] namely

|b j | C j α ,
with α > 0. If we consider such a model for all N with an image verifying Hypothesis 1 then we have

m 2 N u 2 2 = b 2 2 C 2 N j=1 1 j 2α .
This implies the divergence of the sum in the right term, which is thus equivalent to N 1-2α when N goes to infinity. This yields m 2 = O(N -2α ) and so α 0 which is a contradiction. As a consequence, the constant C should depend on N . In the following, we consider the model

|b N j | C N j α ,
and we discuss how to simplify it based on information given by numerical experiments. We need to define how the "constant" C N grows with the image size N .

> 0 if for all N, j ∈ N |b N j | C N γ j α . ( 18
)
where b N = V N u N is the projection of u N on the basis V N .

A case where hypotheses 1 and 2 are trivially satisfied is the case of constant images with a DCT filter. Indeed, the corresponding b is the optimal one from Proposition 2. In the next section, we study the convergence of the upper-bound of the MSE under these hypotheses.

Conditions of convergence

In Section 2.1 we showed that the optimal diagonal estimator on a given basis V N could be bounded in terms of the 1 -norm of the coefficients b N in that basis. In the following, we show that under Hypotheses 1 and 2, this 1 norm can in turn be upper-bounded by a decreasing function of N ,

MSE(λ ) σ 2N b 1 C 1 N r , (19)
thus ensuring convergence of the optimal MSE at a rate r that depends on the parameters α and γ of the decay model. When this rate is positive then we can use this second upper bound to prove the asymptotic optimality of diagonal estimation on that basis.

Theorem 2 (Asymptotic optimality) Consider

an infinite image U that satisfies Hypothesis 1 ( i.e. non-sparsity) and a sequence of orthogonal bases V such that the pair (U, V) satisfies Hypothesis 2 ( i.e. (C, α, γ) decay rate of the image projection on that basis).

If the decay rate is fast enough, i.e. if

1 2 γ < 1 and α > γ,
then the denoising provided by oracle optimization of diagonal estimation on that basis is asymptotically optimal meaning that the MSE tends to 0 when the image size N goes to infinity.

The proof of this result is the combination of the two following Lemmas. The first one shows that the hypothesis on image energy [START_REF] Milanfar | A tour of modern image filtering: New insights and methods, both practical and theoretical[END_REF] constrains the parameters α and γ of the domination criterion and the second one further restricts the values of these parameters to ensure convergence of the upper-bound of the MSE. This lemma emphasizes the fact that we actually cannot bound the |b N j | independently of N as long as we have images that are not loosing energy with increasing size. The only way to obtain γ = 0 (a bound independent of N ) is to impose α = 0 which leads to the pathological case b ∝ (1, . . . , 1). 17) on the image energy gives

Proof (Proof of Lemma 1) Because b N 2 = u N 2 , the model (
m 2 N b N 2 2 M 2 N.
Applying the decay criterion |b N j | C N γ j α in the previous equation yields

m 2 N C 2 N 2γ N j=1 1 j 2α .
The behavior when N goes to infinity of the sum in the right term differs depending on α:

if α < 1 2 the sum diverges and there exists a constant C such that

N j=1 1 j 2α ∼ N →∞ C N 1-2α -if α = 1
2 the sum diverges and there exists a constant C such that This lemma shows that the convergence can actually occur with all α > 0 as long as γ is not too large. We also notice that the model proposed in [START_REF] Talebi | Global denoising is asymptotically optimal[END_REF] in C j α satisfies the convergence hypothesis. However, we saw with that this model is not compatible with Hypothesis 1.

N j=1 1 j 2α ∼ N →∞ C ln N -if α >
Proof (Proof of Lemma 2) We have |b N j | C N γ j α so σ 2N b 1 Cσ 2N N j=1 N γ j α = Cσ 2 N γ-1 N j=1 1 j α .
The behavior when N goes to infinity of the sum in the right term differs depending on α:

if α < 1 the sum diverge and there exists a constant C such that

N j=1 1 j α ∼ N →∞ C N 1-α .
if α = 1 the sum diverges and there exists a constant C such that

N j=1 1 j α ∼ N →∞ C ln N.
if α > 1 the sum converges to a constant C ∈ R.

The first case leads to

σ 2 N b 1 = O N γ-α ,
and convergence occurs only if α > γ. The second and the third cases lead respectively to

σ 2 N b 1 = O N γ-1 ln N , and
σ 2 N b 1 = O N γ-1 ,
and convergence occurs only if γ < 1.

The proof of Lemma 2 also provides a decay rate that we summarize in the following corollary.

Corollary 1 (Decay rate) Under conditions of convergence in Theorem 2, that is 1 2 γ < 1 and α > γ the MSE of optimal diagonal oracle estimation satisfies

MSE(λ ) = N →∞ O 1 N r with r ∈]0, 1 2 ] defined by -r = α -γ when γ < α < 1 -r = 1 -γ when α > 1
The particular case α = 1 yields convergence in O log N N 1-γ .

This result shows that the decay is always slower than 1 √ N and it can be really slow when r is close to zero. Thus, even though we can have an asymptotic optimal filtering, the decay rate can be so small that we cannot actually see it even if we work with huge images. Moreover, this asymptotic study is performed on the oracle diagonal filter. This result is by itself essentially theoretical. However, in combination with Donoho-Johnstone Theorem 1, we might further use this result to prove, under specific conditions on the infinite image, the asymptotic optimality of non-oracle filtering.

Corollary 2 (Decay rate of thresholding) Assume that the convergence conditions of Theorem 2 are satisfied. From the Donoho-Jonhstone Theorem 1, the MSE obtained by thresholding the coefficients b N j satisfies

MSE(λ th ) = N →∞ O log N N r
with r ∈]0, 1 2 ] defined as in Corollary 1. The particular case α = 1 yields convergence in O (log N ) 2 Then by Corollary 1 there exists a constant C such that

MSE(λ ) C N r , with r ∈]0, 1 2 ]. It follows that MSE(λ th ) 2C log N N r + (2 log N + 1) σ 2 N + 2C N r .
The two last terms in the previous inequality are o log N N when N goes to infinity that yield the announced result. A similar proof can be done for the case α = 1.

Special cases

In the following two paragraphs we discuss some simple particular cases in which the asymptotic behaviour of global denoising can be directly deduced.

In more realistic cases we need to experimentally fit our image model to natural images for different bases in order to predict what would happen when the image size tends to infinity. This experimental study is deferred to section 4.

Optimal basis with optimal eigenvalues

When an oracle is used both to choose the optimal basis V and the optimal eigenvalues λ of the filter W = V ΛV T , we showed in Proposition 2 that the optimal M SE decays like σ 2 /N , so we have r = 1, a much faster convergence than in the more realistic cases based on an image model. In this case only b 1 is non zero for all values of N , so computing α and γ does not make any sense.

Gaussian textures on DCT basis

Another case of interest is the case when the image u is a Gaussian texture, meaning that u = h * m is generated by convolving a known kernel h with a white noise image m where m i ∼ N (0, τ 2 ) iid.

In this case when choosing V as a Fourier or DCT basis, a straightforward calculation (see appendix A)

shows that the MSE upper bound in equation ( 14) for global filtering with optimal λ j in that basis becomes

MSE * bound = σ 2 2N b 1 = 1 N 2 σ 2 k | h N (k)|| m N (k)| =: A N .
Thus asymptotically we have a strictly positive MSE bound

MSE * bound = A N ----→ N →∞ A ∞ = στ ĥ 1 √ 2π > 0
for Gaussian textures when using the Fourier or DCT basis. Our experiments (see Section 4.3) confirm this finding. Indeed, when choosing a Fourier or DCT basis V = F , then MSE * N remains constant when N → ∞, in this case r ≈ 0. Nevertheless, when choosing an adaptive basis V from the diagonalization of the non-local means filtering operator, then MSE * N does experimentally tend to 0 for Gaussian textures. This shows that the NLM basis may better exploit the self-similarity in Gaussian textures.

Oracle vs. non oracle filters

Hypothesis 2 on the domination decay model assumes that the family of orthogonal bases V is well adpated to the infinite image U. This happens in particular when the chosen filters are oracle filters, which means that they are computed in the image itself. Consider for instance the case of a very simple oracle filter which consists in denoising ũ by averaging at pixel i all values ũj such that |u i -u j | ≤ ε for a given threshold ε > 0. If the infinite image U is bounded, for instance with values in [0, 1[, we can show that the value MSE( u N |u N ) converges to a limit smaller than ε 2 when the image size increases. This result being satisfied for every ε > 0, the MSE of this oracle filter is naturally asymptotically optimal. A proof of this example is given in appendix B.

The case of non oracle filters is of course far more ambiguous. Consider for instance the case of a dead leaves infinite image model studied in [START_REF] Levin | Patch complexity, finite pixel correlations and optimal denoising[END_REF]. The previous argument shows that a well chosen oracle filter would denoise this image perfectly. However, because of the independence between the leaves, it is clearly not possible to achieve a null asymptotic MSE for a non oracle filter, since for a given leaf, the values observed outside of the leaf are useless to denoise the pixels inside the leaf.

Experiments

In the previous section we introduced a decay model for the b N j coefficients of natural image sequences decomposed on the orthonormal basis given by a symmetric filtering algorithm.

We gave precise conditions on the (γ, α) parameters of this decay model. These conditions may be used to determine whether optimal diagonal estimation on this basis can yield asymptotically optimal denoising performance when applied to a certain family of image sequences.

In practice, answering this question requires to estimate these coefficients from a particular filter/basis based on a truncated image sequence. The next Section 4.1 explains how these model parameters are estimated from real data. Then in Section 4.2 we analyze the asymptotic performance of several denoising algorithms based on the estimated parameters. Finally, in Section 4.3 we provide a discussion about some specific cases that are of particular interest to illustrate our analysis.

Estimating model parameters (C, γ, α)

Theorem 2 gave us a sufficient condition for asymptotic optimality of a filter on an image sequence. This condition is based on the assumption that the |b N j | coefficients follow a particular model, namely:

|b N j | ≈ C N γ j α . ( 20
)
Observing different curves j → |b N j | for various images, sizes N and orthonormal bases in loglog scale (see the first column of Figure 11 for an example), we notice that the model ( 20) holds except for the first few largest coefficients and for a significant proportion of the smallest coefficients. This behaviour can be expected, since we sorted the coefficients. It appears even when the |b N j | coefficients are only white noise (as illustrated in Figure 8). Thus we exclude the values of j < d = 5 and j > N p (for p = 0.6) from the bilinear regression that allows to fit the values of C, α and γ to the |b N j | coefficients.

Put another way we find α, γ and C that minimize

log(|b N j |) -(γ log(N ) -α log(j) + log(C)) 2 ,
with N from N min to N max and j from d to N p .

Experimental results

Table 1 shows the estimated model parameters for the test images from Figure 7 In all cases the oracle NLM basis seems to satisfy the conditions of Theorem 2 and to provide the fastest asymptotic convergence rate. On the other hand, on these experiments, the DCT and wavelet bases some-times seem to not satisfy the conditions of Theorem 2, and when they do, the asymptotic convergence rate is extremely slow (always smaller than r = 0.1) except for sparse images that do not verify Hypothesis 1 and trivially converge in many common bases. This means that if the oracle NLM basis was known for an arbitrarily large noisy image, then we could use hard thresholding as in Corollary 2 to obtain a denoised image with arbitrarily small MSE. Of course the same conclusion was known (since Donoho-Johnstone) for the non-adaptive wavelet and DCT bases, but convergence does not hold for all natural images, and when it does it may be too slow for the procedure to be practical. For oracle NLM asymptotic convergence seems to be faster with respect to image size but we are confronted to two difficulties:

1. the W -oracle is in principle unknown; and 2. diagonalizing an NLM filter is extremely expensive computationally (O(N 3 ) with respect to the number N of pixels).

In order to address the first difficulty we included in our tests the asymptotic performance of the prefiltered NLM. Directly computing the NLM filter on the noisy image is not acceptable as explained in Section 2.2.3. However, applying it to a pre-filtered version of the image helps both (a) to satisfy the requirement of independence of the filter and noise, and (b) to make the filter closer to the oracle one. We show in Table 2 the asymptotic convergence rate we estimated for the pre-filtered NLM basis and for the image lena.

For these experiments we used the denoising algorithm NL-Bayes [START_REF] Lebrun | A Nonlocal Bayesian Image Denoising Algorithm[END_REF] to obtain the pre-filtered image. Moreover, we tuned the parameters of NL-Bayes in order to slightly over-denoise the image. This trick allows to ensure that the filter is almost independent of the noise realization (at the expense of the potential loss of some subtle image structures). The experiment shows that for the P-NLM basis not only do we achieve asymptotic convergence, but the convergence rate is surprisingly close to the convergence rate for the oracle NLM basis. However, as it can be seen in Figure 11 the actual MSE and the bound are always larger for Prefiltered NLM than for Oracle NLM.

Unfortunately this asymptotic behaviour in the nonoracle case cannot be generalized to all natural images. Indeed under certain texture models a lower bound has been established for all possible image denoising algorithms as recalled in section 4.3.3.

However all these model estimates should be taken with a grain of salt, for several reasons: Clearly more experiments on larger images are required to confirm or infirm the conclusions of this initial experimental study. Doing so will require the use of more sophisticated and numerically efficient ways to compute the eigenbasis of the NLM filter on medium to large-size images. This could be achieved by means of randomized numerical linear algebra [START_REF] Halko | Finding Structure with Randomness: Probabilistic Algorithms for Constructing Approximate Matrix Decompositions[END_REF], but such techniques do assume a low rank structure of the filtering matrix, so they cannot be used to estimate the full spectrum of eigenvalues of W . Rather they should be used in conjunction with incremental schemes like in [START_REF] Brand | Fast low-rank modifications of the thin singular value decomposition[END_REF]. This shall be the subject of further research.

Sparse Image

When the image is not widespread (like in the image sparse in Figure 7), it shows trivial internal redundancy that can be exploited by both the DWT and the NLM bases. Hence the behavior of the MSE is dominated by the black part of the image, and for this reason we obtain a very fast decay of the MSE, hence almost matching the theoretical rate r = 1 that is achieved when we use not only optimal eigenvalues λ but also an optimal basis V .

Texture change

In the previous section, we saw experiments for images that do not change drastically with increasing size. But we can wonder what happens when the image suddenly changes with increasing size.

We show in Figure 9 an image composed of two textures and the corresponding curve presenting the behavior of the MSE for the O-NLM case. When we add the second texture, the MSE increases, but when the first texture reappears the MSE starts to decay again. This behavior can explain the fact that there is no need for a stationary hypothesis on the image to obtain convergence. For a sufficiently good basis, able to capture the self-similarity of images, such as NLM-O, we can hope for an asymptotically optimal denoising. This relies on the fact that when the scene size tends toward infinity, we can expect similar structures to reappear again and again in the image.

Gaussian & dead-leaves texture models

To emphasize the importance of the basis we provide a numerical experiment with a synthetic Gaussian texture in Figure 10. We proved in Section 3.4.2 that for such a texture the convergence does not hold for a DCT basis even if this texture presents a lot of self-similarity. The numerical experiment confirms that result and provides experimental evidence for the asymptotic decay of the MSE in the Oracle-NLM basis.

Unfortunately, this positive result for the O-NLM case cannot possibly be extended to the non-oracle case. Indeed Levin et al. [START_REF] Levin | Patch complexity, finite pixel correlations and optimal denoising[END_REF] established a strictly positive lower bound for any image denoising algorithm. This result hols for infinite images that do not present longdistance statistical dependencies. In that case the optimal denoiser for a given pixel x uses the values of the noisy image in pixels y within a neighborhood of x which does not exceed a certain maximal distance D. For y beyond that neighborhood, u(x) and u(y) are independent, so the values of u(y) provide no useful information to estimate u(x). This is the case for Gaussian textures generated by a compactly supported kernel h, and for the dead-leaves model [START_REF] Alvarez | The Size of Objects in Natural and Artificial Images[END_REF]. For images of this kind, Levin's positive lower bound implies that asymptotically zero MSE is impossible to achieve by any non-oracle denoising algorithm. Our experimental result on the Gaussian texture suggests asymptotic convergence of global denoising towards zero MSE, only in the W -oracle case, as seen in section 3.4.3 which explains this mechanism. But this result does not extend to the case where global denoising does not use an oracle to define the filter W .

Conclusion

In this paper, we analyzed the following question:

Can an image denoising algorithm attain asymptotically zero estimation error when the image size tends to infinity? This question was recently raised in [START_REF] Talebi | Global denoising is asymptotically optimal[END_REF][START_REF] Talebi | Global image denoising. Image Processing[END_REF] in the context of oracle-optimized non-local filtering schemes. That work suggests a positive answer but their reasoning is based on conditions on the infinite image that we show incompatible with reasonable assumptions. We refine these conditions to better account for natural images, and provide a more general theory of optimal asymptotic denoising performance. In particular our theory explores how to partially avoid the use of an oracle, it does not restrict itself to global image denoising, and establishes links to the older diagonal estimation theory, as well as with the optimality results of Donoho and Johnstone [START_REF] Donoho | Ideal denoising in an orthonormal basis chosen from a library of bases[END_REF].

More specifically, our work highlights the central role played by the oracle in the work of [START_REF] Talebi | Global denoising is asymptotically optimal[END_REF], and makes a clear distinction between two different ways in which the oracle is used, namely: First a W -oracle is used to construct the entries in the non-local filter W whose diagonalization provides a basis V . Then a λ-oracle is used to find optimal weights {λ * j } for a given basis V . The link we established with diagonal estimation theory means that the λ-oracle can be avoided using Donoho and Johnstone's theorem, meaning that we can study the convergence of a denoising algorithm that uses a λoracle, in order to predict the asymptotic convergence (at a slower rate) of an algorithm that does not use such a λ-oracle.

The W -oracle, however, is more difficult to avoid, since we do not have a tool equivalent to Donoho and Johnstone's theorem in this setting. Hence non-oracle convergence properties need to be directly tested on a version of the algorithm that does not use the Woracle. And this is quite problematic because, without an oracle, special care is required to ensure that the filter W and the noise n are independent. And this independence is required for our asymptotic analysis of the MSE to be valid. The quest for more general ways to define non-local and non-oracle filters W , in a way that their independence from image noise is ensured, is still an open subject for future research.

As a whole our generalized analysis of the asymptotic behaviour of global image denoising provides less optimistic conclusions than those in [START_REF] Talebi | Global denoising is asymptotically optimal[END_REF] but still leaves the door open for asymptotically zero denoising error. Our experimental study on small images seems to indicate that the oracle non-local means filter can be optimized to attain asymptotically zero error, and that a non-oracle version (i.e. without W -oracle) of that filter may have a similar behaviour, even though at a much slower convergence rate and on a more restricted number of examples. Clearly, more extensive experimentation on a wider variety of larger-sized images is required to determine whether these conclusions may have any practical interest. However, performing such an experimental evaluation requires huge amounts of computation, and can only be addressed if faster and more incremental matrix decomposition algorithms are developed.

These conclusions and the prospect of asymptotically zero MSE may appear to be in contradiction with the strictly positive lower bounds for image denoising established by Levin et al. [START_REF] Levin | Patch complexity, finite pixel correlations and optimal denoising[END_REF]. A careful inspection reveals that there is no such contradiction, rather different models and complementary viewpoints that we shall try to clarify below:

The positive lower bound of Levin et al. is valid for certain statistical image models such as Gaussian and dead-leaves textures For images of this kind, Levin's positive lower bound implies that asymptotically zero MSE is impossible to achieve by any non-oracle denoising algorithm. Our experimental result on the Gaussian texture suggests asymptotic convergence of global denoising towards zero MSE, only in the W -oracle case. But this result does not extend to the case where global denoising does not use an oracle to define the filter W .

For more decidedly self-similar images like lena or bricks our experiments indicate that even the P -N LM filter that does not use a W -oracle is compat-ible with asymptotically zero MSE. This shows that for this image Levin's assumption of absence of longdistance dependencies does not hold, otherwise there would be a contradiction.

The quest for a statistical model for natural images that takes self-similarity into account in a realistic way is still a very active area of research. Extending such a model for image sizes tending to infinity poses yet an additional challenge. Future research in that direction would hopefully allow to unify Levin's and Talebi's views on asymptotic behaviour of image denoising. Proof Consider for simplicity images of square size N = n 2 with n ∈ Z 2 . Now consider m N the restriction of m to

A Global denoising for Gaussian textures on Fourier basis

I N = [-N 2 , N 2
2 )∩ Z 2 . In the truncated case the convolution is understood in the periodic sense, so that we can write

u N (k) = h N (k) • m N (k) ∀k ∈ I N with h N (k) = ĥ( 2πk n ). Now if we decompose u N in the Fourier basis V = F * then b N k = (F u N ) k = h N (k) • m N (k)
So the upper bound given in equation ( 14) for the optimal MSE is Proof Indeed, let n = 1 ε , for 0

MSE * N ≤ 1 N 2 σ b N 1 2 = 1
≤ k < n let C k = {i ∈ Ω; u N i ∈ [ k n , k+1
n [}, and for every pixel i let

A ε i = {j; |u N i - u N j | ≤ ε}. It is clear that if i ∈ C k , then C k ⊂ A ε i . We can write MSE( u N |u N ) = 1 N n-1 k=0 i∈C k E[|u N i -u N i | 2 ] = 1 N n-1 k=0 i∈C k E[|u N i - 1 #A ε i j∈A ε i ũN j | 2 ] = 1 N   n-1 k=0 i∈C k |u N i - 1 #A ε i j∈A ε i u N j | 2 + n-1 k=0 i∈C k σ 2 #A ε i   ≤ ε 2 + σ 2 N n-1 k=0 i∈C k 1 #A ε i ≤ ε 2 + σ 2 N n-1 k=0 i∈C k 1 #C k ≤ ε 2 + nσ 2 N .
For a given value of ε, this term converges to ε 2 when N → ∞.

C Additionnal experiments

Figures 12 through 16 show the detailed asymptotic convergence results for the images in Figure 7 and Table 1.

2 = u 2 . 2

 222 ) depends only on b = V T u, the projection of the oracle image u on the basis V . The following Proposition describes the form of the b minimizing (6). The matrix V T being orthonormal, the minimization is constrained by b Proposition Minimizing b → MSE (b) under the constraint b 2 = u 2 provides the following 2N global minimums b = ± u 2 e i where e i is the i-th vector of R N basis. Proof (Proof of Proposition 2) Let define the function Ψ by

Fig. 1 (

 1 Fig. 1 (a) Original images u, (b) Noisy images u with σ = 15, (c) Images u denoised by hard-thresholding in a DCT basis with a threshold T = σ √ 2 ln N , (d) Images u denoised by diagonal estimation with the oracle λ * i in a DCT basis.

  . a Non-Local Means filter [5] computed on the original image u (called Oracle-NLM or O-NLM); 2. a Non-Local Means filter computed on the noisy image u (called NLM); 3. a Non-Local Means filter computed on a version of u already denoised by NL means (called pre-filtered NLM or P-NLM);

Fig. 2

 2 Fig. 2 Relative error (13) between the theoretical MSE provided by formula (4) and the experimental MSE, for three images and three different filters W (NL-Means computed on the u, NL-Means computed on the u, NL-Means computed on a prefiltered version of u. The mean and standard deviation have been computed on 5 different realisations of noise with σ = 15 for each image.

Fig. 3

 3 Fig. 3 Example of denoising with a global MSE decaying to zero when image size grows to infinity, while a local MSE increases. (b) is an image constructed by repeating N times a pattern of a constant image with a vertical line. (a)presents the behavior of the MSE between the denoised image obtained by optimal diagonal estimation in a DCT basis and the clean image. The MSE is shown in log 10 scale, when the pattern is repeated N times with N increasing. (c) presents a zoom on the structured part of the image and the result of the denoising for N = Nmax. The denoised image presents important ringing artifacts. Finally, (d) shows the behavior of the local MSE on the part presented in (c) when N grows. This shows that even with a global MSE converging to zero, the restoration can remain locally bad and even get worse when the image size increases.

Fig. 4

 4 Fig.4Decay rate of the coefficients b j for the image synthetic from Figure7, in a loglog scale graph. In blue: coefficients in a DCT basis. In red: coefficients in a wavelet basis. In yellow: coefficients in the eigenbasis of the oracle NLM filter.

2

 2

3. 2

 2 Proposed models Infinite image model Consider an image of infinite size U : Z 2 → {l, . . . , L}, taking values in a discrete set of gray levels {l, . . . , L} ⊆ N. For typical 8 bit images l = 0 and L = 255.From this image we construct an infinite sequence of images of growing size Nu N def = U ϕ(1) , . . . , U ϕ(N ) ,by truncating the infinite image to size N , for all N ∈ N. The function ϕ : N + → Z 2 sweeps the plane in spiral starting from the origin. Since the image gray level values are bounded, the L 2 -norm of u N satisfies the following inequality l means that the energy of the growing image increases at most like O( √ N ). This information on the L 2 -norm of u N is important because it constrains the behavior of b N as b N 2 = u N 2 . Because generally the lowest value l is zero, the energy of the image may not grow as fast as √ N . However we show that if u N 2 = o( √ N ), the image is becoming sparse with increasing size: because U is taking values in a discrete finite set, by setting c = min {U i = 0, i ∈ N} we have

Fig. 5

 5 Fig. 5 Behaviour of max j (|b N j |) with increasing size N for three different images from Figure 7 in loglog scale. Here b N is the image u projected in the DCT basis.

Fig- ure 5

 5 shows that max j |b N j | is increasing linearly with N in loglog scale. That leads us to consider N γ as a model for C N . Finally, we consider the following decay model for b N : Hypothesis 2 (Domination decay model) Let U be an infinite image and denote by u N its truncation of size N . Let V = (V N ) N be a family of orthogonal bases of increasing size. Then the pair (U, V) is said to fit the domination decay model with parameters C, α and γ

Fig. 6

 6 Fig. 6 Illustration of the domain of compatibility and the domain of convergence provided by the two lemmas 1 and 2. The intersection in red represent the set of parameters that provides the result in Theorem 2.

Figure 6

 6 Figure 6 illustrates the results provided by Lemma 1 and lemma 2 on the parameters α and γ. The resulting parameters for Theorem 2 are given by the intersection of the two domains. Lemma 1 (Compatibility with image model) Assume that U satifies Hypothesis 1. If the projection of U on V satisfies the decay model of Hypothesis 2 with parameters (C, α, γ), then either γ 1 2 and α ≥ 1 2 or γ α and α < 1 2 .

1 2 2 .

 12 the sum converges to a constant C ∈ R The first case leads to m 2 = O N 2γ-2α and so α γ. The second and the third cases lead to m 2 = O N 2γ-1 ln N and m 2 = O N 2γ-1 respectively and so γ 1 Lemma 2 (condition of convergence) Considering the model (18) we have convergence to zero of the bound (14) only if α 1 and γ < 1 or α < 1 and α > γ.

N 1

 1 -γ . Proof (Proof of Corollary 2) Let consider the case α = 1. By Donoho-Jonhstone Theorem we have MSE(λ th ) (2 log N + 1) σ 2 N + 2 MSE(λ ) .

Fig. 7

 7 Fig. 7 The images used for the experiments. The sub-images sizes are (a) 128 × 128, (b) 128 × 192, (c) 128 × 256, (d) 128 × 320, (e) 128 × 384, (f) 128 × 448, (g) 128 × 512. Images credits: lena and man are standard images used in image processing, simpson is from Julie Delon, brick, sparse and mixed are Brodatz textures [4] and synthetic is a random generated gaussian texture.

Fig. 8

 8 Fig. 8 Decay of the coefficients |b j | for white noise in the DCT basis (red) in loglog scale. The slope of the bound is αm ≈ 0.05 (blue).

  and for three orthonormal bases, namely: DCT: The DCT basis which diagonalizes convolution filters; Wavelet: The orthogonal Haar basis, implemented via the discrete wavelet transform; Oracle NLM: The orthogonal basis which diagonalizes the oracle (symmetrized) non-local means filter, i.e. with patch distances computed on the oracle clean image.

Table 2 -

 2 Fitted parameters α and γ for lena image in the three bases DCT, O-NLM and P-NLM. The parameter r is the decay rate of corollary 1. The cubic computational cost of exactly computing the eigenbasis of the NLM filters obliged us to limit our evaluation to relatively small image sizes. -Model (20) can not always be perfectly fit by all images and bases. The model seems to hold for "stationary" images or for images that contain a relatively small number of stationary components. Otherwise the task of fitting this model is particularly difficult. -Model (20) only gives a coarse upper bound for the actual MSE * . The second column in Figures 12 through 16 shows that even though this upper bound is relatively coarse, the actual MSE * does follow the same kind decay with N as the upper bound. Nevertheless, when comparing the actual MSE * of all four bases (Figure 11) we observe that the real performance of the prefiltered NLM is actually comparable to that of DCT or wavelet bases; even though the convergence rate r estimated on this model (0.234 for P-NLM vs 0.06 for DCT, 0.10 for DWT and 0.236 for O-NLM) seemed to indicate that the prefiltered NLM was much superior to DCT and rather close to the oracle NLM performance.

4. 3

 3 Fig. 9 Top: Image mixed with two different textures. Bottom: the corresponding MSE for diagonal estimation using O-NLM basis when the size grows from (a) to (g) in log 10 log 10 scale.

Fig. 10

 10 Fig. 10 Top: Image synthetic. Bottom: the corresponding MSE for diagonal estimation using DCT basis and O-NLM basis when the size grows from (a) to (g) in log 10 log 10 scale.

Fig. 11

 11 Fig. 11 Left column: the decay of the |b N j | for each size and the result of the model fitting (dotted lines) for the image lena for the different bases (from top to bottom) DCT, Wavelet, O-NLM and P-NLM. Right column: the decay of the MSE (blue), the upper bound from (14) (orange) and the fitted bound (yellow).

Proposition 3 2 2N b 1 = 1 N 2 σ 2 k

 32122 Let u = h * m be a Gaussian texture where m i ∼ N (0, τ 2 ) iid, and the kernel h∈ 1 (Z 2 ) has a smooth Fourier transform h ∈ L 1 ([-π 2 , π 2 ] 2 ) ∩ C ∞ .Then when choosing V as a Fourier or DCT basis, and λ as the optimal oracle eigenvalues for that basis, then the MSE for global denoising is upper bounded byMSE * bound = σ | h N (k)|| m N (k)| =: A NAsymptotically we have a strictly positive MSE boundMSE * bound = A N ----→

|Proposition 4

 4 h N (k)|| m N (k)| =: A NIn the asymptotic case when N → ∞, a simple Riemann sum argument based on the regularity of ĥ and on the knownvalue of dI | m| = |dI|E( m N (k)) = τ 2π |dI| over a small interval (mean absolute deviation of a gaussian), leads to the conclusion that the upper bound tends to a strictly positive constantMSE * ∞ ≤ lim N →∞ A N = στ h 1 √ 2π = A∞ > 0.B An Oracle filter that provides asymptotically zero MSE Consider an oracle filter which consists in denoising ũ by averaging at pixel i all values ũj such that |u i -u j | ≤ ε for a given threshold ε > 0 and an infinite image U bounded with values in [0, 1[. Then the value MSE( u N |u N ) converges to a limit smaller than ε 2 .

Fig. 12

 12 Fig. 12 Left column: the decay of the |b N j | for each size and the result of the model fitting (dotted lines) for the image simpson for the different bases (from top to bottom) DCT, Wavelet, and O-NLM. Right column: the decay of the MSE (blue), the upper bound from (14) (orange) and the fitted bound (yellow).

Fig. 13

 13 Fig. 13 Left column: the decay of the |b N j | for each size and the result of the model fitting (dotted lines) for the image bricks for the different bases (from top to bottom) DCT, Wavelet, and O-NLM. Right column: the decay of the MSE (blue), the upper bound from (14) (orange) and the fitted bound (yellow).

Fig. 14

 14 Fig. 14 Left column: the decay of the |b N j | for each size and the result of the model fitting (dotted lines) for the image sparse for the different bases (from top to bottom) DCT, Wavelet, and O-NLM. Right column: the decay of the MSE (blue), the upper bound from (14) (orange) and the fitted bound (yellow).

Fig. 15

 15 Fig. 15 Left column: the decay of the |b N j | for each size and the result of the model fitting (dotted lines) for the image synthetic for the different bases (from top to bottom) DCT, Wavelet, and O-NLM. Right column: the decay of the MSE (blue), the upper bound from (14) (orange) and the fitted bound (yellow).

Fig. 16

 16 Fig. 16 Left column: the decay of the |b N j | for each size and the result of the model fitting (dotted lines) for the image man for the different bases (from top to bottom) DCT, Wavelet, and O-NLM. Right column: the decay of the MSE (blue), the upper bound from (14) (orange) and the fitted bound (yellow).

Table 1

 1 Fitted parameters α and γ for the different images of figure 7 in the three bases DCT, DWT and O-NLM. The parameter r is the decay rate of corollary 1.

	image	basis	fitted α fitted γ	r
		DCT	0.827	0.767	0.060
	lena	DWT	0.858	0.752	0.10
		O-NLM	0.806	0.570	0.236
		DCT	0.941	0.741	0.201
	simpson	DWT	1.106	0.753	0.247
		O-NLM	1.025	0.572	0.428
		DCT	0.796	0.755	0.041
	bricks	DWT	0.913	0.877	0.035
		O-NLM	0.902	0.688	0.213
		DCT	0.842	0.375	0.467
	sparse	DWT	0.909	0	0.909
		O-NLM	1.099	0.021	0.979
		DCT	0.708	0.727	-
	synthetic	DWT	0.554	0.570	-
		O-NLM	0.646	0.540	0.106
		DCT	0.720	0.520	0.200
	man	DWT	0.802	0.578	0.224
		O-NLM	0.759	0.366	0.393

From now on, we will write b N j instead of b j to remember that the behavior of these coefficients strongly depend on the image size N .

The authors would like to thank H. Talebi and P. Milanfar for encouraging and insightful discussions, and the anonymous reviewers for critical observations that helped to clarify and improve the paper. This work has been partially funded by the French Research Agency (ANR) under grant no. ANR-14-CE27-001 (MIRIAM). Most of this work was done while A. Almansa worked with CNRS LTCI at Télécom ParisTech.