Camille Sutour
email: camille.sutour@math.u-bordeaux.fr.

Charles-Alban Deledalle
email: cdeledal@math.u-bordeaux.fr.

Jean-François Aujol
email: jaujol@math.u-bordeaux.fr.

Estimation of the noise level function based on a non-parametric detection of homogeneous image regions *

Keywords: Noise level estimation, non-parametric detection, signal-dependent noise, least absolute deviation AMS subject classifications

published or not. The documents may come

Estimation of the noise level function based on a non-parametric detection of homogeneous image regions

Camille Sutour, Charles-Alban Deledalle, Jean-

1. Introduction. Most image processing applications, including denoising [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF][START_REF] Sutour | Adaptive regularization the NL-means: Application to image and video denoising[END_REF], super-resolution [START_REF] Freeman | Example-based super-resolution[END_REF], or feature extraction [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF], require the knowledge of the noise level. Noise estimation from multiple images is an over-constrained problem that was solved in [START_REF] Healey | Radiometric ccd camera calibration and noise estimation[END_REF], but its estimation from a single image remains a challenging problem. We review in the following some of the main approaches that have been proposed for the estimation of the noise level when the noise is spatially uncorrelated. Methods dedicated to correlated noise can be found, e.g, in [START_REF] Colom | A non-parametric approach for the estimation of intensity-frequency dependent noise[END_REF][START_REF] Portilla | Blind non-white noise removal in images using gaussian scale mixtures in the wavelet domain[END_REF][START_REF] Lebrun | The noise-clinic: a universal blind denoising algorithm[END_REF].

A first approach used for the estimation of the noise level consists in applying a linear transform that maps the image into a suitable domain where the signal and the noise are relatively well separated. The noise level can then be estimated by performing some statistics in its associated subspace. For instance, when the noise is additive white Gaussian, a popular estimator of the noise standard deviation is given by the mean absolute deviation of the finest wavelets coefficients, see e.g. [START_REF] Donoho | De-noising by soft-thresholding[END_REF]. Nevertheless, this approach tends to overestimate the noise level as it assumes there is enough separation between the signal and the noise coefficients. Indeed, singularities such as edges usually produce high responses in all sub-bands of the transformed domain. To avoid this issue, some authors [START_REF] Olsen | Estimation of noise in images: an evaluation[END_REF][START_REF] Tai | A fast method for image noise estimation using laplacian operator and adaptive edge detection[END_REF] have proposed to pre-filter the image before estimating the noise in order to remove such high frequency components, e.g., using a Laplacian-based filter together with a Sobel edge detector [START_REF] Rank | Estimation of image noise variance[END_REF]. Unfortunately, such approaches still lead to a subsequent over-estimation of the noise level. Rather than attempting to remove the high frequency components, some authors have proposed to take them into account by assuming some prior knowledge on the distribution of the signal coefficients. In [START_REF] De Stefano | Training methods for image noise level estimation on wavelet components[END_REF], the authors use a method of "trained" moments from an assumption of Laplacian distribution of the signal coefficients, while the authors of [START_REF] Zoran | Scale invariance and noise in natural images[END_REF] exploit the invariance of the kurtosis when assuming a generalized Gaussian distribution. Patch-based approaches have also been proposed for such a task, e.g., in [START_REF] Shin | Block-based noise estimation using adaptive gaussian filtering[END_REF][START_REF] Liu | Single-image noise level estimation for blind denoising[END_REF][START_REF] Pyatykh | Image noise level estimation by principal component analysis[END_REF]. While such techniques can be very powerful, they are nevertheless limited to constant noise levels. In real-world applications, the noise is signal-dependent and its level varies spatially according to the underlying intensity. Variance stabilizing transforms [START_REF] Anscombe | The transformation of poisson, binomial and negative-binomial data[END_REF] can reduce the dependency between the signal intensity and the noise. The authors of [START_REF] Pyatykh | Image sensor noise parameter estimation by variance stabilization and normality assessment[END_REF] iteratively perform the estimation of Poisson-Gaussian noise by applying a variance stabilizing transform, then they estimate the noise variance of the stabilized noise using a principal component analysis. In [START_REF] Boulanger | Patch-based nonlocal functional for denoising fluorescence microscopy image sequences[END_REF], Poisson-Gaussian noise is stabilized using the generalized Anscombe transform. However, this transform requires a prior knowledge of the noise parameters, obtained by linear regression on areas with homogeneous variance [START_REF] Gasser | Residual variance and residual pattern in nonlinear regression[END_REF]. The generalized Anscombe transform is also optimized for the stabilization of Poisson-Gaussian noise in [START_REF] Makitalo | Noise parameter mismatch in variance stabilization, with an application to poisson-gaussian noise estimation[END_REF].

Separation techniques have been extended to specific signal-dependent models, e.g., using a wavelet transform for a Poissonian-Gaussian model [START_REF] Foi | Practical Poissonian-Gaussian noise modeling and fitting for single-image raw-data[END_REF][START_REF] Azzari | Gaussian-cauchy mixture modeling for robust signal-dependent noise estimation[END_REF] or using a Gaussian mixture model of patches for additive noise with affine variance [START_REF]Indirect estimation of signal-dependent noise with nonadaptive heterogeneous samples[END_REF]. Unfortunately, in the general case, it appears pointless to design a linear transform able to separate the noise from the signal. The noise can also be distinguished from the signal components by principal component analysis [START_REF]Analysis and extension of the pca method, estimating a noise curve from a single image[END_REF] or by selecting blocks with lowest variance [START_REF] Colom | Analysis and extension of the percentile method, estimating a noise curve from a single image, Image Process[END_REF]. An alternative consists in pre-filtering the image to analyze the residue as if it was the noise component. For instance, [START_REF] Hensel | Modeling and real-time estimation of signal-dependent noise in quantum-limited imaging[END_REF][START_REF] Chen | An automated technique for image noise identification using a simple pattern classification approach[END_REF] have proposed different pre-filtering techniques for some specific types of signal-dependent noise. However, for our general problem, as such prefiltering techniques require a certain knowledge of the noise, this leads to a chicken and egg problem.

Another popular approach encountered in [START_REF] Chehdi | A new approach to identify the nature of the noise affecting an image[END_REF][START_REF] Beaurepaire | Identification of the nature of noise and estimation of its statistical parameters by analysis of local histograms[END_REF] and that we follow here relies on the fact that natural images contain homogeneous areas, where the signal to noise ratio is very weak, so only the statistics of the noise intervene. The statistics of the noise can then be deduced by first detecting homogeneous areas and next estimating the noise parameters in those regions. However, classic detectors also require assumptions on the noise statistics (e.g, Gaussian), preventing their use when the nature of the noise is completely unknown.

Our main contribution is to propose a non-parametric detection of homogeneous areas that only requires the noise to be spatially uncorrelated. This means that the detection performance of our detector, based on Kendall's τ coefficient [START_REF] Kendall | A new measure of rank correlation[END_REF], is independent of the noise statistical distribution. In Section 3, we introduce our homogeneous detector based on the rank correlation. We show that the performance of the detector relies on the underlying signal to noise ratio, with an error rate that decreases with the number of samples available. Once homogeneous areas are detected, we estimate the noise level function (NLF) [START_REF] Liu | Noise estimation from a single image[END_REF], i.e., the function of the noise variance with respect to the image intensities. The NLF estimation is detailed in Section 4. At this step we assume that the noise is signal-dependent with finite first and second order moments linked together by a second order polynomial. This can encompass a wide range of noise statistics encountered in image processing such as additive, multiplicative or Poissonian noise, as well as hybrid models. In the last step we perform a robust estimation of this polynomial NLF as the best ℓ 1 approximation of the mean/variance relationship of homogeneous regions. Section 5 details how the estimation is performed in practice, then Section 6 shows the efficiency of the detection of homogeneous areas and of the parameters estimation for different noise models, by comparing them to the state of the art. The interest of performing a reliable noise estimation is also illustrated with an application to image denoising: the non-local means algorithm [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF] is adapted to deal with hybrid noises, thanks to the estimated NLF parameters. Finally, Section 7 provides discussions and conclusions.

2. Context and notations. We denote by g : Ω → R an image defined on a 2D discrete domain Ω ⊂ Z × Z such that g can be identified to an element of R N 0 with N 0 = #Ω the number of pixels of the image. Hence, for an image g ∈ R N 0 , i ∈ [1, . . . , N 0 ] is the index of its i-th pixel and g i ∈ R the value of g at pixel i.

We assume that g is an observation of a clean unknown image g 0 ∈ R N 0 corrupted by a spatially uncorrelated signal dependent noise. More precisely, we assume that g can be modeled as the realization of a random vector G that verifies (2.1)

E[G] = g 0 and Cov[G] =      NLF(g 0 1 ) 0 NLF(g 0 2 ) . . . 0 NLF(g 0 N 0 )      .
where NLF : R → R + is coined the noise level function. As Cov[G] is assumed to be diagonal, this model restricts to spatially uncorrelated noise. But since the i-th diagonal element depends on g 0 i , the noise is said to be signal dependent. Remark, that the same noise level function applies for each diagonal elements, so the noise is assumed to be stationary. Noise models with spatially varying noise level functions are out of the scope of this study.

Example 2.1 (Additive white Gaussian noise).

A classic example is to consider an additive white Gaussian noise such that G = g 0 + ε, with ε ∼ N (0, σ 2 Id). In such a case, G satisfies eq. (2.1) with NLF(•) = σ 2 which reveals that the noise is signal independent.

Example 2.2 (Poisson noise).

Another example is to consider a Poisson noise such that g 0 ∈ N, G 0 i ∼ P(g 0 i ) and, for all i = j, G i and G j are independent. In such a case, G satisfies eq. (2.1) with NLF(g 0 i ) = g 0 i which reveals that the noise is signal dependent.

Example 2.3 (Multiplicative noise). Another example is to consider a multiplicative noise such that

G i = g 0 i × S i , with E[S] = 1 and E[SS ⊤ ] = γ Id, γ > 0.
In such a case, G satisfies eq. (2.1) with NLF(g 0 i ) = γ(g 0 i ) 2 which reveals again that the noise is signal dependent. A typical example is the gamma noise of parameter L ∈ N * , for which γ = 1 L , whose probability density function f is defined for

g i ∈ R + as f (g i | g 0 i , L) = L L g L-1 i Γ(L)(g 0 i ) L exp - Lg i g 0 i . (2.2)
The goal of our work is to estimate the unknown noise level function NLF. To do so, we consider K small disjoint blocks (ω k ) k∈[1,...,K] ⊂ Ω of size N such that KN ≤ N 0 . We denote by ω any of such block ω k and by b = g ω ∈ R N (resp. B = G ω and b 0 = g 0 ω ) the restriction of g in ω (resp. the restriction of G and g 0 in ω). Our motivation relies on the fact that most natural images, provided that their resolution is high enough, contain small blocks ω for which b 0 is constant. In this case we say that B is homogeneous, and otherwise when b 0 varies, B is said to be inhomogeneous. Hence, according to eq. (2.1), the empirical expectation and variance inside such a block b should provide a punctual estimation of the noise level function.

In practice b 0 is unknown as well as the homogeneous blocks. We hence propose in Section 3, a statistical test that can decide weather a block is likely to be homogeneous or not without more assumptions regarding the noise model introduced in (2.1). Once homogeneous regions are detected, several punctual estimations of the noise level function are available. Section 4 describes how these punctual estimates can be used to robustly estimate a noise level function with few additional assumptions, by fitting a second order polynomial function to the scatterplot.

3. Detection of homogeneous areas. The noise statistics can be estimated from several small homogeneous areas, or in fact areas that are close to be homogeneous, i.e. where the signal to noise ratio is weak. Indeed, in those regions the signal fluctuations can be assumed to be negligible compared to the noise fluctuations, so only the statistics of the noise intervene. Hence, provided that a large number of such regions has been detected, the statistics of each regions can be used to infer the underlying noise level function. The goal of this section is then to develop a method that automatically selects homogeneous areas in the image, based on the fact that natural images contain such regions.

Homogeneity and uncorrelation.

In the light of detection theory, we aim at designing a detector of inhomogeneous blocks that answers the following statistical hypothesis test from the observation b only

H 0 : B is homogeneous, i.e., b 0 constant (null hypothesis), H 1 : B is inhomogeneous, i.e., b 0 varies (alternative hypothesis). (3.1)
As b 0 is unknown, we cannot directly provide an answer to this test. To alleviate this difficulty, we suggest making use of this first proposition.

Proposition 3.1. Let B ∈ R N be homogeneous. We denote by B I the sub random vector of B whose entries are indexed by I ⊂ {1, . . . , N }. Then, for any disjoint subsets I and J of n elements of {1, . . . , N }, 2n ≤ N , X = B I and Y = B J are uncorrelated, i.e.:

Corr(X, Y ) = 0 .

This proposition tells us that if two random disjoint sequences X and Y of a block are correlated, then the block is inhomogeneous. In practice, the random variables are unknown and only their single realizations, the noisy sequences x = b I and y = b J , are available. As a consequence, the correlation between X and Y is also unknown and the arising question is "Are the realizations x and y significantly correlated?".

Hence, the goal is to design a score s : R n × R n → R that can be used to answer correctly to our hypothesis problem. More precisely, the score s should ensure that, for any constant P FA ∈ [0, 1[, there exists a threshold α > 0 such that, whatever the distribution of X and Y , the followings are satisfied

lim n→∞ P(|s(X, Y )| ≥ α false alarms | H 0 ) = P FA , (C1) lim n→∞ P(|s(X, Y )| < α miss detections | H 1 ) = 0. (C2)
The first condition is very important and states that as the size n gets large enough, the same threshold α can maintain the same probability of false alarms P FA as soon as X and Y are uncorrelated (regardless of the distribution of X and Y ). The test is said to be an asymptotically constant false-alarm rate (CFAR) detector. The second condition is natural and states that with the same threshold α the probability of miss detections vanishes as the size n of the sequences increases.

In the following, we propose to build such a statistical test between pairs of disjoint sequences of a block. As in practice, we cannot test every disjoint sequence, we focus on the correlation between neighboring sequences that we denote by x = (b 2k ) and y = (b 2k+1 ), where 2k and 2k + 1 represent neighbor pixel indexes for a given scan path of the underlying block ω. An example of the construction of the two neighboring sequences x = (b 2k ) and y = (b 2k+1 ) is displayed on Figure 1. If these two variables are found to be significantly correlated, this means that there is some dependencies between the pixels of the blocks and their neighbors, so we can assume that there is some structure and all fluctuations are not only due to noise.

A naive approach for answering the above hypothesis test is to use the empirical correlation between the two sequences x and y, i.e., s(x, y) = Corr(x, y). Unfortunately, the significance of such a test depends on the parameters of the random process that generates the observation, i.e., it cannot lead to a CFAR detection without extra assumptions on the distribution of X and Y . As our goal is to estimate the nature of the noise, it is important for our test of homogeneity not to rely on such parameters. We therefore consider a non-parametric approach whose statistical answer is independent of the noise model.

Non parametric test of independence.

As mentioned above, the empirical correlation test does not provide a CFAR detector as the threshold α would depend on the unknown parametrization of the random process that generates the observation. The key idea of non-parametric tests is to focus on the rank (i.e., on the relative order) of the values rather than on the values themselves. In this vein, a block is considered homogeneous if the ranking of the pixel values is uniformly distributed, regardless of the spatial organization of the pixels. Such correlation between the sequences x and y can be assessed based on the ranking of their values using Kendall's τ coefficient.

3.2.1. Kendall's τ coefficient. The Kendall's τ coefficient is a rank correlation measure [START_REF] Kendall | A new measure of rank correlation[END_REF] that provides a non-parametric hypothesis test for statistical dependence. Before turning to the definition of the Kendall's τ coefficient, let us first introduce the definitions of concordant, discordant and tied pairs. Definition 3.2. Two pairs of observations (x i , y i ) and (x j , y j ) are said to be concordant if the ranks of both elements agree, i.e., if (x i < x j and y i < y j ) or (x i > x j and y i > y j ). They are said to be discordant if (x i < x j and y i > y j ) or if (x i > x j and y i < y j ). If x i = x j or y i = y j , the pair is neither concordant nor discordant; it is said to be tied.

We can now turn to the definition of the Kendall's τ coefficient as introduced in [START_REF] Kendall | A new measure of rank correlation[END_REF]. Definition 3.3. Let x ∈ R n and y ∈ R n be two sequences without tied pairs. Kendall's τ : R n × R n → [-1, 1] coefficient is the quantity given by

(3.2) τ (x, y) = n c -n d n(n-1) 2 = 1 n(n -1) 1≤i,j≤n sign(x i -x j ) sign(y i -y j ),
where n c is the number of concordant pairs and n d the number of discordant pairs. Remark that Kendall's τ (x, y) coefficient depends only on the relative order of the values of x and y. It lies in the range -1 ≤ τ (x, y) ≤ 1. If the agreement between the two rankings is perfect, all pairs are concordant, then τ (x, y) = 1. If the disagreement between the two rankings is perfect, all pairs are discordant, then τ (x, y) = -1. A value τ (x, y) = 0 indicates the absence of significant correlation between x and y.

Definition 3.3 is however limited to sequences without tied pairs. This restrict the measure of correlation to the case where the unknown random process is continuous, e.g., Gaussian, and the apparition of tied pairs is hence of zero measure. In our imaging context, for instance with optical systems, the image formation processes typically include a photon counting step that is by nature discrete. Then ties in the sequence might appear with non-zero probability so it is important to take them into account properly.

When the data suffers from tied pairs, the definition of Kendall's τ coefficient might be extended according to [START_REF]The treatment of ties in ranking problems[END_REF] as follows.

Definition 3.4. Let x ∈ R n and y ∈ R n be two arbitrary sequences. Define from x and y the quantities

n 0 = n(n -1)/2
the total number of pairs, n 1 = i t i (t i -1)/2 the total number of tied values for the first quantity, n 2 = j u j (u j -1)/2 the total number of tied values for the second quantity, t i the number of tied values for the i th group for the first quantity, u j the number of tied values for the j th group for the second quantity.

Kendall's τ : R n × R n ∈ [-1, 1] coefficient is the quantity given by

(3.3) τ (x, y) = n c -n d (n 0 -n 1 )(n 0 -n 2 ) .
Remark that in the absence of ties, Definition 3.4 matches with Definition 3.3. Equipped with such a measure of correlation, we are now ready to study its ability to answer our hypothesis test under the requirements (C1) and (C2).

3.2.2. Distribution of τ in homogeneous blocks (under the null hypothesis H 0 ). In this section, we study the distribution of Kendall's τ coefficient under the null hypothesis H 0 and show accordingly that it can lead to a statistical test satisfying (C1). In particular as it is based only on the relative order of the values of x and y, it provides a non-parametric measure of correlation leading to a CFAR detector.

We first give the distribution of τ in the absence of tied pairs. Proposition 3.5. Let X and Y be two arbitrary sequences under H 0 without tied pairs. The random variable τ (X, Y ) has an expected value of 0 and a variance of 2(2n+5) 9n(n-1) . Moreover, in case of large samples n, its distribution is approximated by the normal distribution. More precisely, for any sequences x and y, let z : R n × R n ∈ R be the z-score defined as

(3.4) z(x, y) = 3(n c -n d ) n(n -1)(2n + 5)/2 .
The z-score is asymptotically distributed according to the standard normal distribution z(X, Y ) ∼ n→∞ N (0, 1) .

The proof can be found in [START_REF]Rank correlation methods[END_REF]. We now turn to the extension in the case of tied pairs. Proposition 3.6. For any sequences x ∈ R n and y ∈ R n , let z : R n × R n ∈ R be the z-score defined as

(3.5) z(x, y) = n c -n d √ v with v = (v 0 -v t -v u )/18 + v 1 + v 2 v 0 = n(n -1)(2n + 5) v t = i t i (t i -1)(2t i + 5) v u = j u j (u j -1)(2u j + 5) v 1 = i t i (t i -1) j u j (u j -1)/(2n(n -1)) v 2 = i t i (t i -1)(t i -2) j u j (u j -1)(u j -2) 9n(n-1)(n-2) .
. Then, for any sequences X and Y under H 0 , the z-score is asymptotically distributed according to the standard normal distribution

z(X, Y ) ∼ n→∞ N (0, 1) .
The proof can be found in [START_REF]The treatment of ties in ranking problems[END_REF]. A direct consequence of Proposition 3.6 is the following. Corollary 3.7. Let X and Y be two arbitrary sequences under H 0 , then

(3.6) lim n→∞ P(z(X, Y ) > α | H 0 ) = +∞ α 1 √ 2π exp - t 2 2 dt = 1 -φ (α) ,
where φ is the cumulative distributive function (cdf ) of the Gaussian distribution

(3.7) φ(x) = 1 √ 2π x -∞ e -t 2 2 dt.
The proof is immediate and comes from the definition of the probability of false-alarm given in (C1), combined with the normal distribution of z(X, Y ) under H 0 . We now turn to the first theorem that shows that the z-score satisfies the condition (C1) as its distribution does not rely on any assumption regarding the distribution of X and Y .

Theorem 3.8. Let X and Y be two arbitrary sequences under H 0 . Let P FA ∈ [0, 1[ and choose

(3.8) α = φ -1 (1 -P FA /2) = √ 2 erf -1 (1 -P FA ) ,
where erf is the Gauss error function, erf(x) = 1 √ π

x -x e -t 2 dt, then

(3.9) lim n→∞ P(|z(X, Y )| > α | H 0 ) = P FA .
The proof is given in appendix A. Besides, apart from the standard distribution approximation, the z-score still satisfies the CFAR condition (C1) for low numbers of samples n, as the distribution remains the same whatever the noise distribution.

In practice, the detection is performed by computing the associated p-value p(x, y) = 2 -2φ(|z(x, y)|), (3.10) i.e., the probability under the null hypothesis of having |z(X, Y )| greater than |z(x, y)|.

We then reject the null hypothesis if the p-value is smaller than the predetermined significance level of false alarm P FA . In general, the smaller the p-value, the larger the significance of the rejection of the null hypothesis H 0 .

3.2.3. Distribution of τ subject to a step-signal (an alternative hypothesis H 1 ). The computation of the distribution of τ under the null hypothesis reflects the non-parametric properties of Kendall's test, since no assumption on the distribution of X and Y (except the absence of spatial correlation) has been made.

When X and Y are somehow correlated, a model of dependence has to be assumed in order to derive some statistical properties for τ . Accordingly, we propose to study the ability of Kendall's τ coefficient to accept or to reject the alternative hypothesis H 1 when X and Y correspond to sequences extracted from a 1D step signal corrupted by additive white Gaussian noise. See Fig. 1 for the illustration of such a signal. Definition 3.9. Let n be an even size of sequence and H 1 be the hypothesis of dependence between X and Y parametrized by a > 0 and σ 2 as follows

∀ 1 ≤ k ≤ n, X k = ε k if k ≤ n/2 a + ε k if k > n/2 and Y k = η k if k ≤ n/2 a + η k if k > n/2 ,
where ε and η are independent random variables distributed along N (0, σ 2 ). Proposition 3.10. Let X and Y be two sequences satisfying H 1 . Then

(3.11) E[τ (X, Y ) | H 1 ] = n 2(n -1) 1 -2φ a √ 2σ 2 ,
and in particular

(3.12) lim a σ →+∞,n→+∞ E[τ (X, Y ) | H 1 ] = 1 2 .
The proof is detailed in appendix B. This proposition shows that, asymptotically to the dimension n, the expected value of τ in the case of a step tends towards 1/2 when the signal to noise ratio (a versus σ) is large enough. In other words, since the expected value of τ under the null hypothesis is zero, the two distributions are expected to be well separated as soon as the sequences are large enough and the signal to noise ratio is sufficient. It remains to show that the variance of the sampling distribution of τ decreases with the length of the sequence. Proposition 3.11. Let X and Y be two sequences satisfying H 1 . Then

Var[τ (X, Y ) | H 1 ] = O 1 n .
The proof is detailed in appendix C. This proposition guarantees that as the number of samples n increases, the variance of the sample distribution of τ under H 1 decreases towards zero. Since the expected value of τ under H 1 is strictly positive for a > 0, this statement provides that for a sufficient number of samples, the distributions of τ under both cases can be distinguished. The next corollary extends this result to the case of the z-score defined in (3.4),(3.5).

. Empirical distributions of Kendall's z-score under the null hypothesis H 0 (blue) and under the H 1 hypothesis of statistical dependence (red) for a growing number of samples (n = 32, 128, 512). The threshold α associated to the probability of false alarm P FA = 0.15 is displayed as a green line. On the top line, the noise power is small (σ = 40) compared to the step size (a = 128), resulting in a signal to noise ratio a/σ = 3.2, whereas on the bottom line, the noise power is stronger (σ = 60) compared to the step size (a = 64), resulting in a lower signal to noise ratio a/σ = 1.07, which makes the detection trickier.

Corollary 3.12. Let X and Y be two sequences satisfying H 1 . Then

(3.13) E[z(X, Y ) | H 1 ] = O( √ n) and Var[z(X, Y ) | H 1 ] = O (1) .
The proof is detailed in appendix D.

As the distribution of z under H 0 is zero mean with a standard deviation of 1, this corollary guarantees that as the number of samples n increases, the z-score can be used to discriminates between H 0 and H 1 . Indeed, as a consequence of Chebyshev's inequality, the next theorem shows that the z-score satisfies condition (C2).

Theorem 3.13. Let X and Y be two sequences satisfying H 1 . For any α > 0, we have

(3.14) P(z(X, Y ) < α | H 1 ) ≤ Var(z(X, Y )|H 1 ) (E[z(X, Y )|H 1 ] -α) 2 = O 1 n ,
and in particular, as

P(|z(X, Y )| < α | H 1 ) ≤ P(z(X, Y ) < α | H 1 ), we have (3.15) P(|z(X, Y )| < α | H 1 ) = O 1 n .
The proof is detailed in appendix E.

Illustrations. Figure 2 displays the empirical distributions of Kendall's z-score under the null hypothesis H 0 (blue) and the alternative hypothesis H 1 of statistical dependence in case of a step (red) for a growing number of samples (n = 32, 128, 512). The threshold α associated to the probability of non detection P ND = 0.15 is displayed as a green line. On the top line, the noise power is small (σ = 40) compared to the step size (a = 128) so the signal to noise ratio a/σ = 3.2 is important, whereas on the bottom line, the noise power is stronger (σ = 60) compared to the step size (a = 64), resulting in a lower signal to noise ratio a/σ = 1.07, which makes the detection trickier. These empirical distributions show that the error rate is related to the signal to noise ratio, i.e. the relationship between the noise power σ and the step size a, due to the proximity of both distributions. Indeed, on the bottom line where both values are close, according to Proposition 3.10 the expected value of τ , hence z, under H 1 is close to zero, hence to the expected value under H 0 , which makes the two distributions harder to distinguish. However, when the number of samples n increases, Corollary 3.12 guarantees that the expected value of z under H 1 grows towards the infinity, enabling to tell them apart more precisely.

Figure 3 illustrates condition (C2). The empirical probability of miss-detection is displayed as a function of the number of samples, in a log-log scale. The theoretical bound we derived is shown to decrease asymptotically as 1/n, while the empirical probability of miss-detection seems to offer an even better decrease rate.

Extension to homogeneous areas in images.

The above study shows that Kendall's coefficient is able to distinguish uniform regions from areas that contain an edge. Based on the fact that natural images can be approximated by a piece-wise constant model, this guarantees that homogeneous regions can be detected, and the error is controlled by the number of samples and the signal to noise ratio.

In practice, we run K = 4 Kendall's τ test for four pair of sequences (x (1) , y (1) ), (x (2) , y (2) ), (x (3) , y (3) ), (x (4) , y (4) ) corresponding respectively to horizontal, vertical and the two diagonal neighbors, as shown on Figure 4. We could also investigate other relationships, for example not directly related neighbors but k-neighbors (located at a distance of k pixels), for every k, or look at these correlations in another domain, for example in the Fourier domain, to expose frequency dependencies. The area is considered to be homogeneous if each of the K obtained p-values p(x (k) , y (k) ) reaches a given level of significance α. By doing so, the overall level of detection α eq after aggregation is no longer α but smaller and given by (3.16)

α eq = P K k=1 p(X (k) , Y (k) ) > α .
In order to control the overall level of detection α eq , we empirically estimated offline the relation between α eq and α.

Figure 5 illustrates the selection of homogeneous areas. We synthesize a hybrid noise as the sum of Gaussian, Poissonian and multiplicative gamma noise, resulting in an initial PSNR of 17.93dB. We compute Kendall's τ coefficient within non-overlapping blocks of size N = 16 × 16, by comparing horizontal, vertical, and diagonal neighbors. Figure 5b displays the p-value associated to the rank correlation computed between horizontal neighbors. Then we perform the selection by thresholding the p-value, in order to obtain an overall level of detection given by the probability of detection P D = 1 -P FA = 0.7. The blocks that have been considered homogeneous are displayed in red on Figure 5-c.

Model estimation.

Once the homogeneous blocks B = G ω are detected, the noise statistics can be estimated on those blocks. More precisely, for an homogeneous block B of size N , Eq. (2.1) implies that there exists µ ∈ R and

σ 2 > 0 such that E[B i ] = E[B j ] = µ and Var[B i ] = Var[B j ] = σ 2 for all i, j ∈ [1, . . . , N ].
Again as B is unknown, we can only relies on empirical statistics based on the observation b that unbiasedly estimate µ and σ 2 as follows

µ(b) = 1 N N i=1 b i and σ 2 (b) = 1 N -1 N i=1 (b i -µ) 2 , (4.1)
meaning that E[ µ(B)] = µ and E[ σ 2 (B)] = σ 2 . Accordingly, each homogeneous block gives a punctual estimation of the the noise level function defined in (2.1) as σ 2 (b) ≈ NLF( µ(b)). In this section, we explain how to robustly estimate the noise level function given P pairs of statistics ( µ(b p ), σ 2 (b p )), p ∈ [1, . . . , P ] from P detected homogeneous blocks b p = g ω p .

Before turning to the estimation of the noise level function, it is important to understand the behavior of the above variance estimator inside a single homogeneous region.

4.1. Punctual estimation from a single homogeneous area. The unbiased estimator of the variance σ 2 computed inside a block B of N pixels might suffer from statistical errors. This section aims at studying the performance of the variance estimator in the best possible case, i.e., when all fluctuations in b are only due to the noise (no false alarms).

More precisely, we are interested in quantifying the relative error

| σ 2 (B)-σ 2 | σ 2
under H 0 . To do so, we first introduce the following proposition.

Proposition 4.1. Let B be a homogeneous block satisfying H 0 . We have

(4.2) Var( σ 2 (B) | H 0 ) ∼ N →∞ σ 4 N β 2 - N -3 N -1 ,
where β 2 is the kurtosis given, for any i ∈ [1, . . . , N ], by the following formula

(4.3) β 2 ∼ N →∞ E B i -µ σ 4 H 0 .
The proof can be found in [START_REF] Cho | Variance of sample variance[END_REF]. Not surprisingly, this proposition shows that as N increases, the variance vanishes, and as σ 2 is an unbiased estimator, σ 2 (B) tends towards σ 2 . However, it shows also that the variance and hence the efficiency of the unbiased estimator of σ 2 depend on the variance σ 2 itself and the kurtosis β 2 of the unknown noise process that generates b.

We are now ready to introduce our third theorem. Theorem 4.2. Let B be a homogeneous block satisfying H 0 . The relative error in estimating the variance σ 2 has an expected value and a variance given by

E σ 2 (B) -σ 2 σ 2 H 0 ∼ N →∞ 1 N β 2 - N -3 N -1 × 2 π , (4.4) Var σ 2 (B) -σ 2 σ 2 H 0 ∼ N →∞ 1 N β 2 - N -3 N -1 × 1 - 2 π . (4.5)
The proof is detailed in appendix F. This formula gives the expected value and the variance of the estimation error when computing the variance on a single homogeneous block.

Interestingly, like for the variance of σ 2 (B), this theorem shows that the relative error vanishes with high probability as the block size N increases. But, unlike the variance of σ 2 (B), the relative error under H 0 solely depends on the kurtosis and the size of the block irrespectively of the variance σ 2 to be estimated. This shows that, for a given size of block N , the estimation error can then be predicted as a function of the kurtosis, as displayed on Figure 6-a where N = 16 × 16. Hence, without any assumptions on the underlying variance, one can bound the expected relative error made by our approach under H 0 provided some mild assumptions on the kurtosis. For instance, assuming the kurtosis is lower than 40 (that corresponds to a very high peakedness of the noise random process and so such high kurtosis is very unlikely to appear in practical imaging problem), the expected relative error will be lower than 30%.

The following examples support our claim in the case of noise distributions classically encountered in imaging problems.

Example 4.3 (Additive white Gaussian noise).Consider the example of an additive white Gaussian noise of parameter σ 2 as given in Example 2.1. In this case, the kurtosis β 2 is constant with respect to σ 2 and equal to 3. Provided the block is of size N = 16 × 16, Theorem 4.2 shows that the expected relative error in estimating σ 2 is of about 7%. the expected value that tends towards zero at the same rate, and the green line shows the variance estimation error at each bloc size for different noise realizations.

Example 4.4 (Multiplicative gamma case).Consider the example of the multiplicative gamma noise of parameter L as given in Example 2.3. In this case, the kurtosis is given by the formula β 2 = 6 L + 3. Here, contrary to the Gaussian case, the estimation error depends on the parameter of the distribution L. In the worst case when L = 1, Theorem 4.2 shows that for a block of size N = 16 × 16 the expected relative error in estimating the variance is of about 14%.

Example 4.5 (Poisson case).Consider the example of the Poisson noise as given in Example 2.3. In this case, the kurtosis varies spatially with the image intensity, and is given by β 2 = 1 g 0 i + 3, where g 0 i is the underlying intensity. for example, in a pessimistic case, where min i∈ω g 0 i = 0.04 for a block of size N = 16 × 16, which corresponds in an imaging problem of a pixel with a very low expected number of photons, Theorem 4.2 shows that the expected relative error in estimating the variance is of about 28%.

Besides, the estimation error computed here gives a bound on the error that affects the estimation of the variance on a single block, i.e., for one point of the NLF. The next step consists in fitting a NLF to the set of estimated points, so the overall error based on all the selected blocks should be reduced by the number of selected blocks, which depends on the image size and the desired level of false alarms.

However, when some false detections arise, i.e., if inhomogeneous blocks are selected, the variance will be overestimated, resulting in a higher estimation error.

The bounds computed above hence give a safety check for the error one might expect: if the overall error ends up being higher than the predicted bound, this means either that too many outliers were taken into account or that the NLF estimator did not fit properly.

4.2.

Noise level function estimation from all homogeneous areas. Provided that uniform regions have been detected, the noise model can be estimated from the mean/variance couples extracted from these areas. We now explain how to robustly estimate a noise level function (NLF) as defined in Eq. 2.1 given P pairs of statistics (μ p , σ2 p ) = ( µ(b p ), σ 2 (b p )), p ∈ [1, . . . , P ] from P detected homogeneous blocks b p = g ω p . Within the p-th homogeneous block, the intensity variations should be only ascribed to noise fluctuations, hence σ 2 p ≈ NLF( µ p ). As each block provides only a punctual infor-mation about the NLF, we are confronted to an under constraint approximation problem from P points of R × R + to a function of R → R + . It is then necessary to inject some extra regularity assumptions. We suggest that the NLF is a positively increasing second order polynomial of the image intensity, parametrized by θ = (a, b, c) t ∈ (R + ) 3 and given by

σ 2 p = NLF θ (µ p ) = aµ 2 p + bµ p + c. (4.6)
According to examples 2.1, 2.2 and 2.3, this model encompasses Gaussian noise, with (a, b) = (0, 0), in which case the NLF is constant, Poisson noise, with (a, c) = (0, 0), in which case the NLF is linear, and multiplicative noise, with (b, c) = (0, 0), in which case the NLF is parabolic, as well as hybrid models, e.g., Poissonian-Gaussian. Note that the positivity constraints can be dropped to embrace more general models to allow non increasing noise level function.

The goal is to find the polynomial coefficients (a, b, c) such that the estimated variances σ2

p can be well approached as NLF θ (μ p ) = aμ 2 p + bμ p + c. Denoting by μ ∈ R P , σ2 ∈ (R + ) P and NLF θ (μ), the vectors obtained by stacking each of these estimations, our approximation problem can be expressed in a vectorial form as

σ2 ≈ NLF θ (μ) = Aθ where A =      μ2 1 μ1 1 μ2 2 μ2 1 . . . . . . . . . μ2 P μP 1      (4.7)
Of course, if the noise is known to be Gaussian, Poissonian or multiplicative respectively, the operator A should be adapted, keeping only a single column of ones, of the means or of the means to the square respectively. In case of the Poisson-Gaussian noise, only the two columns of ones and of the means should be kept.

We focus in the following on two different estimators in order to robustly estimate θ.

4.2.1. Least squares estimation. The least squares (LS) method consists in minimizing a ℓ 2 -norm of the residue NLF θ (μ) -σ2 . The estimation is formalized as the solution of the follow constraint optimization problem (4.8) θLS = argmin

θ∈(R + ) 3 NLF θ (μ) -σ2 2 = argmin (a,b,c)∈(R + ) 3 aμ 2 + bμ + c -σ2 2
whose solutions can be obtained using quadratic programming tools. Remark that without the positivity constraint on (a, b, c) the least square solution is explicit and the following holds (4.9)

θ LS = (A t A) -1 A t σ2 and NLF LS θ (μ) = Π A σ2
where Π A = A(A t A) -1 A t is the so-called hat function that is the projector on the space of second order polynomial. Assuming σ2 would have approximately a normal distribution with covariance proportional to the identity, the above least square estimator would correspond to the maximum likelihood estimator that is guaranteed to converge in probability towards the underlying second order polynomial NLF function. However, even though σ2 has a diagonal covariance matrix (as the blocks are independent), it is not proportional to the identity. Worst, the distribution of σ2 deviates from normality as the data might suffer from outliers. Indeed, large errors can arise when some non-homogeneous samples happen to be selected. This leads us to studying other estimation methods more robust to outliers.

Algorithm 1 A preconditioned primal-dual algorithm of [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] for least absolute deviation Set γ ∈ R P,1 and τ ∈ R 1,P such that γ i = 1/ 3 j=1 |A ij |, and

τ j = 1/ P i=1 |A ij |. Π |z|≤1 the projector of R P such that |z i | ≤ 1 ∀1 ≤ i ≤ p, Π θ≥0 the projector of R 3 such that θ j ≥ 0 ∀1 ≤ j ≤ 3, Initialize z 0 = 0 ∈ R P , θ0 = 0 ∈ R 3 and θ 0 = 0 ∈ R 3 , Iterate for k ≥ 0    z k+1 = Π |z|≤1 z k + γ(Aθ k -σ) , θ k+1 = Π θ≥0 (θ k -τ A * z k+1 ) , θ k+1 = 2θ k+1 -θ k .

Least absolute deviation.

If miss-detections happen as inhomogeneous samples are selected, the presence of structures in the signal might result in a highly over estimated variance resulting from the summation of the noise variations with the signal variations. In order to reduce the impact of such outliers, we suggest using the least absolute deviation (LAD) method. Unlike the LS estimation, the LAD estimation is shown to be more robust to outliers [START_REF] Pollard | Asymptotics for least absolute deviation regression estimators[END_REF].

The LAD estimation is based on the minimization of the ℓ 1 -norm of the residue (instead of the ℓ 2 -norm) as the solution of the following constraint problem (4.10) θLAD = argmin

θ∈(R + ) 3 NLF θ (μ) -σ2 1 = argmin (a,b,c)∈(R + ) 3 aμ 2 + bμ + c -σ2 1 .
There is no close form solution for such a problem even in the un-constrained case. Fortunately, we can derive a solution based on a fast iterative algorithm, using the preconditioned primal-dual algorithm of Chambolle and Pock [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF], as displayed in Algorithm 1 and where Π is the projection operator on a convex set.

While the LS estimator corresponds to the maximum likelihood estimator in the case of estimated variances following a specific normal distribution, the LAD estimator is the maximum likelihood estimator assuming a specific Laplacian distribution. Unlike the Gaussian distribution, the Laplacian distribution has heavy tails meaning that it is suitable to describe data subject to extreme values or outliers. According to the maximum likelihood interpretation of the LAD estimator, errors are assumed to follow the same Laplacian distribution, which is not verified in practice, see Theorem 4.2. And even if it were Laplacian, the scale of the Laplacian should be rescaled with respect to the unknown underlying variance as suggested by Eq. 4.4. However, due to the performance of the homogeneous detector the scatterplot is reliable (the amount of false detections is reasonable and the false detections are not too aberrant). The LAD estimator hence offers a fast and reliable estimation that is robust to outliers.

A more realistic model would be to use distributions with heavier tails, such as the Cauchy distribution, as proposed in [START_REF] Azzari | Gaussian-cauchy mixture modeling for robust signal-dependent noise estimation[END_REF]. This leads to optimization problems that are typically non convex, and hence harder to minimize. The Laplacian distribution is in this sense a good trade-off between robustness and computational efficiency.

Examples of other robust estimators include Theil-Sen's estimator [START_REF] Theil | A rank-invariant method of linear and polynomial regression analysis[END_REF] (also referred to as Kendall robust line-fit method) used to fit affine functions based on median slopes, its generalization as the quantile regression [START_REF] Koenker | Quantile regression[END_REF] for finding curves corresponding to the conditional median or other quantiles of the data and Hubert's M-estimators [START_REF] Huber | Robust estimation of a location parameter[END_REF] designed to be robust under partial knowledge of the data distribution. Figure 8. Mean relative error as a function of the desired probability of detection P D = 1 -P FA . On the left, the noise is Poissonian, with different levels. It can be estimated with the prior knowledge of the noise model, using the least squares or the least absolute deviation adapted to the noise model (Poisson LS or Poisson LAD), or it can be described by the general second order model with the hybrid estimators, LS or LAD. On the right, the noise is hybrid with different NLF parameters. Hybrid noise are estimated with the general second order model using the least square and the least absolute deviation estimators. Whereas the noise model is known or not, the LAD estimator is shown to be more accurate that the LS estimator.

Comparison of estimators.

To compare the least absolute deviations estimator with the least squares estimator, and for all the following experiments, we have generated a set of noisy images with different noise levels, using a database of 150 clean grayscale high quality images of resolution 540×720, taken from a 7 Mpix Canon Powershot S70 or a 8 Mpix Canon IXY Digital 910IS1 . The noise on these images can be considered to negligible, as shown on the sample of this database displayed on Figure 7. Based on the knowledge of the real noise parameters θ = (a, b, c), we can compute the mean relative error

(4.11) MRE( θ) = 1 |I| f i ∈I NLF θ (f i ) -NLF θ (f i ) NLF θ (f i ) ,
where I is a discretization of the interval of image intensities. Figure 8 displays the evolution of the mean relative error as a function of the desired level of detection P D = 1 -P FA . On Figure 8-a while the other two parameters a and c are set to zero. We then compare such estimators to the least square (LS) and least absolute deviation (LAD) estimators that estimate the three parameters (a, b, c) of a general second order model without the knowledge of a Poissonian model. On Figure 8-b, the images have been corrupted by hybrid noise with varying NLF parameters (a, b, c). In this case, a general second order model is directly estimated.

4.4. Parameter estimation results. Figure 9 displays the estimation of the NLF with the least absolute deviations estimator for different noise models (Gaussian, Poisson, multiplicative gamma and hybrid). The mean/variance couples extracted from the homogeneous patches are shown in gray crosses, the real NLF is displayed in green and the estimated one using the least absolute deviation estimator appears in red. There is no prior on the noise model, so the estimation is performed by seeking a second order polynomial NLF. The results fit closely, and the shape of each estimated curve shows that even with no prior knowledge of the noise level function, the model is accurately guessed. The knowledge of the noise model can make the parameters estimation even more accurate.

Practical settings.

5.1. Automatic adaptation of the block size N . The optimal block size has to satisfy a compromise between the reliability of the detection and the performance of the estimation. Indeed, both Kendall's test (see Section 3.2.3) and the computation of the statistics on each block (see Section 4.1) need a sufficient number of samples, dictated by the block size. But if the block size is too big then the number of homogeneous zones might be too small or even null, hence limiting the performance of the NLF estimation. To make sure that the estimation remains efficient or even possible, we impose a minimum number of detected areas. We divide the intensities of the input image into three bins of equal size, and we force the detection step to find at least three regions in each bin. If this condition is not satisfied, then the block width is reduced by two pixels in each direction, until the number of detected regions is sufficient. This guarantees that most images can be processed, and that the initial block size can be set to be large enough to limit the error in Kendall's test and in the computation of the statistics.

Selection of the maximum block size.

The maximal block size can then be set as the size that balances the compromise between the reliability of the detection and the performance of the estimation. We have studied on a dataset of images corrupted by Gaussian noise the evolution of the mean relative error as a function of the block width W , where the block size N = W × W , displayed on Figure 10. Since we cannot in this case automatically reduce the block size in order to guarantee a minimum number of regions, we impose as before a minimum of three detected regions in each bin by selecting the ones that have the best Kendall score, even if the associated p-value is below the threshold. Again, this was only done to determine offline on our dataset the maximum size window. This guarantees a sufficient number of data for the estimation. However, if the block size is too big the data might be irrelevant, meaning that the statistics would arise from non-homogeneous areas. The optimum appears to be reached around W = 16 pixels.

Selection of the probability of detection.

To set the optimum probability of detection, rules by the threshold α, we have used the results obtained from Figure 8. Indeed, the convex shape of the error curves shows two opposite influences. At low levels of detection, the number of selected areas is small. False detections do not arise very likely, but the NLF estimation is not finely estimated as the number of data is insufficient. When the level of detection increases, the number of samples for the estimation increases, which reduces the estimation error, but the number of false detections increases as well, leading to the presence of outliers. The outliers are responsible for the increasing error at higher levels of detection, and they also explain the superiority of the least absolute deviation estimator over the least square estimator, at high levels of detection in particular. This leads to a trade-off between detection and outliers, indicating that for a level of detection ≈ 60%, a relative error lower than 10% can be reached without prior knowledge of the noise model. Accordingly, in all the following experiments, we chose the parameter P D = 0.6.

Experiments.

In this section, we discuss and compare the efficiency of our proposed approach in terms of NLF estimation errors and image denoising applications. For the sake of replicability, a Matlab implementation for the automatic Table 1 Estimation error under low Gaussian noise, with increasing variance. The estimators are MAD [START_REF] Donoho | De-noising by soft-thresholding[END_REF], the Gaussian-Cauchy mixture model [START_REF] Azzari | Gaussian-cauchy mixture modeling for robust signal-dependent noise estimation[END_REF], the PCA method [START_REF]Analysis and extension of the pca method, estimating a noise curve from a single image[END_REF], the percentile method [START_REF] Colom | Analysis and extension of the percentile method, estimating a noise curve from a single image, Image Process[END_REF], Noise Clinic [START_REF] Colom | A non-parametric approach for the estimation of intensity-frequency dependent noise[END_REF], and our algorithm. In our method, we automatically adapt the block size to provide enough homogeneous zones, so the average block size for each level is displayed. σ = 1 σ = 2 σ = 5 σ = 10 MAD [START_REF] Donoho | De-noising by soft-thresholding[END_REF] 3.263 1.546 0.621 0.314 Gaussian-Cauchy [START_REF] Azzari | Gaussian-cauchy mixture modeling for robust signal-dependent noise estimation[END_REF] 55.767 11.281 1.619 0.4669 PCA [START_REF]Analysis and extension of the pca method, estimating a noise curve from a single image[END_REF] 0.376 0.163 0.136 0.235 Percentile [START_REF] Colom | Analysis and extension of the percentile method, estimating a noise curve from a single image, Image Process[END_REF] 0.421 0.181 0.077 0.0378 Noise Clinic [START_REF] Colom | A non-parametric approach for the estimation of intensity-frequency dependent noise[END_REF] 0

noise estimation and its application to image denoising is available for download at https://github.com/csutour/RNLF. We validate the proposed approach with respect to the state of the art algorithms that perform noise estimation, in terms of efficiency, accuracy, and validation for other processing tasks. Based on the database of 150 natural images, we generate a set of noisy images, for different noise statistics (Gaussian, Poisson, Poisson-Gaussian, hybrid...) and levels. To our knowledge, few state of the art methods can deal with such a general noise model as the one we propose here, so we estimate the noise parameters with the different estimators, according to the noise model involved and the properties of each estimator: the mean absolute deviation (MAD) estimator [START_REF] Donoho | De-noising by soft-thresholding[END_REF], which is only suitable for the estimation of Gaussian noise, the Gaussian-Cauchy mixture model [START_REF] Azzari | Gaussian-cauchy mixture modeling for robust signal-dependent noise estimation[END_REF] which is the most general model, the PCA method [START_REF]Analysis and extension of the pca method, estimating a noise curve from a single image[END_REF], the percentile method [START_REF] Colom | Analysis and extension of the percentile method, estimating a noise curve from a single image, Image Process[END_REF], and Noise Clinic [START_REF] Colom | A non-parametric approach for the estimation of intensity-frequency dependent noise[END_REF][START_REF] Lebrun | The noise-clinic: a universal blind denoising algorithm[END_REF], that estimate frequency-dependent noise but that we use here for the estimation of affine noise, the estimation based on the variance stabilization transform (VST) [START_REF] Pyatykh | Image sensor noise parameter estimation by variance stabilization and normality assessment[END_REF] that applies for Poisson-Gaussian noise, and our algorithm that can estimate either a given model or a general second order one.

Estimation of low Gaussian noise.

A fine estimation of very low Gaussian noise is hard to achieve. In our case, this is due to the difficulty of finding homogeneous areas. Indeed, according to the characterization of homogeneous areas that we define in the introduction, a region is homogeneous if its signal to noise ratio is weak. However, with such low noise levels, the signal fluctuations cannot be assumed to be negligible compared to the noise variations, so the detection is harder to achieve.

Table 1 gives the mean relative error when estimating very low Gaussian noise (σ = 1, 2, 5, 10) with the different estimators adapted for Gaussian noise: MAD [START_REF] Donoho | De-noising by soft-thresholding[END_REF], the Gaussian-Cauchy mixture model [START_REF] Azzari | Gaussian-cauchy mixture modeling for robust signal-dependent noise estimation[END_REF], the PCA method [START_REF]Analysis and extension of the pca method, estimating a noise curve from a single image[END_REF], the percentile method [START_REF] Colom | Analysis and extension of the percentile method, estimating a noise curve from a single image, Image Process[END_REF], Noise Clinic [START_REF] Colom | A non-parametric approach for the estimation of intensity-frequency dependent noise[END_REF], and our algorithm. It does indeed increase when the noise level gets too small. In the proposed method, we automatically adapt the block size to provide enough homogeneous zones, which is more difficult at smaller noise levels. Table 1 also displays in the bottom line the average optimal block size, which has to be reduced under lower noise levels.

Estimation of strong Poisson noise.

Estimating pure Poisson noise can also be challenging, due to the behavior of the noise, especially in the darkest areas when it can appear as shot noise. Figure 11

Table 2

Estimation error under strong Poisson noise, with different powers. The estimators are MAD [START_REF] Donoho | De-noising by soft-thresholding[END_REF], the Gaussian-Cauchy mixture model [START_REF] Azzari | Gaussian-cauchy mixture modeling for robust signal-dependent noise estimation[END_REF], the PCA method [START_REF]Analysis and extension of the pca method, estimating a noise curve from a single image[END_REF], the percentile method [START_REF] Colom | Analysis and extension of the percentile method, estimating a noise curve from a single image, Image Process[END_REF], Noise Clinic [START_REF] Colom | A non-parametric approach for the estimation of intensity-frequency dependent noise[END_REF], the VST based method [START_REF] Pyatykh | Image sensor noise parameter estimation by variance stabilization and normality assessment[END_REF], and our algorithm. [START_REF]Analysis and extension of the pca method, estimating a noise curve from a single image[END_REF] 0.681 0.788 0.933 1.000 Percentile [START_REF] Colom | Analysis and extension of the percentile method, estimating a noise curve from a single image, Image Process[END_REF] 0.279 0.371 0.477 0.834 Noise Clinic [START_REF] Colom | A non-parametric approach for the estimation of intensity-frequency dependent noise[END_REF] 0.430 0.457 0.411 0.998 VST [START_REF] Pyatykh | Image sensor noise parameter estimation by variance stabilization and normality assessment[END_REF] 0.281 \ \ \ Proposed (Poisson) 0.053 0.053 0.072 0.238 Proposed (hybrid) 0.121 0.131 0.106 0.433 on the top line. On the bottom line, the noise parameter are estimated first assuming a Poisson model (on the left), then with the general second order model (on the right). In both cases, the relative error is below 5%, showing the reliability of the estimation. Table 2 gives the mean relative error when estimating strong Poisson noise (b = 60, 100, 200, 1000) with the different estimators adapted for Poisson noise: the Gaussian-Cauchy mixture model [START_REF] Azzari | Gaussian-cauchy mixture modeling for robust signal-dependent noise estimation[END_REF], the PCA method [START_REF]Analysis and extension of the pca method, estimating a noise curve from a single image[END_REF], the percentile method [START_REF] Colom | Analysis and extension of the percentile method, estimating a noise curve from a single image, Image Process[END_REF], Noise Clinic [START_REF] Colom | A non-parametric approach for the estimation of intensity-frequency dependent noise[END_REF], the VST based method [START_REF] Pyatykh | Image sensor noise parameter estimation by variance stabilization and normality assessment[END_REF], which does not work for very high levels, and our algorithm, which offers the most reliable estimation results.

6.3. Comparison on general noise models. Tables 3 to 6 illustrate the estimation performance of the suitable estimators for commonly encountered noise models: Gaussian 3, Poisson-Gaussian 4, multiplicative gamma noise 5 and hybrid noise 6. The different estimators that are used are the mean absolute deviation (MAD) estimator [START_REF] Donoho | De-noising by soft-thresholding[END_REF], which is only suitable for the estimation of Gaussian noise, the Gaussian-Cauchy mixture model [START_REF] Azzari | Gaussian-cauchy mixture modeling for robust signal-dependent noise estimation[END_REF] which is the most general model, the PCA method [START_REF]Analysis and extension of the pca method, estimating a noise curve from a single image[END_REF], the percentile method [START_REF] Colom | Analysis and extension of the percentile method, estimating a noise curve from a single image, Image Process[END_REF], and Noise Clinic [START_REF] Colom | A non-parametric approach for the estimation of intensity-frequency dependent noise[END_REF][START_REF] Lebrun | The noise-clinic: a universal blind denoising algorithm[END_REF], that can estimate frequency-dependent noise but are used here for the estimation of affine noise, and the estimation based on the variance stabilizing transform (VST) [START_REF] Pyatykh | Image sensor noise parameter estimation by variance stabilization and normality assessment[END_REF] that applies for Poisson-Gaussian noise. We compare these methods to our propose algorithm that can estimate either a given model or a general second order one.

In each case, we can either estimate the suitable model with the involved parameters, e.g. constant for Gaussian noise, affine for Poisson-Gaussian, etc, or estimate a second order polynomial NLF, which corresponds to the 'hybrid' case. Results show that the proposed method offers an accurate noise level estimation. It is quite competitive for affine noise, and it offers state of the art estimation results for second order noise models, in particular multiplicative gamma noise.

6.4. Application to image denoising. Automatic noise estimation can be useful for many applications that require the knowledge of the noise level, in particular image denoising. In order to illustrate the use of the proposed method, the estimated NLF is re-injected for blind image denoising.

The non-local means algorithm (NL-means) [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF] has been adapted in order to use only the knowledge of the noise parameters (a, b, c) provided by the estimated NLF. The NLmeans algorithm is based on the natural redundancy of the image structures, not just locally but in the whole image. While local filters average pixels that are spatially close to each other, the NL-means algorithm compares patches, i.e., small windows extracted around each pixel, in order to average pixels whose surroundings are similar. For each pixel i ∈ Ω, the solution of the NL-means is (6.1)

u NL i = j∈Ω w i,j g j j∈Ω w i,j
, where the weights w i,j ∈ [0, 1] are computed in order to select the pixels j whose surrounding ρ j ⊂ Ω is similar to the surrounding ρ i ⊂ Ω of the central pixel i, for instance, as proposed in [START_REF] Sutour | Adaptive regularization the NL-means: Application to image and video denoising[END_REF] by (6.2)

w i,j = exp - |d(g ρ i , g ρ j ) -m ρ d | s ρ d where m ρ d = E[d(G ρ i , G ρ j )] and s ρ d = Var[d(G ρ i , G ρ j )]
are respectively the expectation and the standard deviation of the dissimilarity d computed between two noisy patches of size |ρ| following the same distribution.

In order to take into account the signal-dependence of the noise, we approximate the noise to additive Gaussian noise with polynomial variance, and we adapt the dissimilarity measure d as follows

(6.3) d(g ρ i , g ρ j ) = 1 |ρ| |ρ| k=1 g ρ i k -g ρ j k 2 NLF θ (g ρ i k ) + NLF θ (g ρ j k )
.

Table 3 Mean relative error (MRE) when estimating Gaussian noise with MAD [START_REF] Donoho | De-noising by soft-thresholding[END_REF], the Gaussian-Cauchy mixture model [START_REF] Azzari | Gaussian-cauchy mixture modeling for robust signal-dependent noise estimation[END_REF], the PCA method [START_REF]Analysis and extension of the pca method, estimating a noise curve from a single image[END_REF], the percentile method [START_REF] Colom | Analysis and extension of the percentile method, estimating a noise curve from a single image, Image Process[END_REF], Noise Clinic [START_REF] Colom | A non-parametric approach for the estimation of intensity-frequency dependent noise[END_REF][START_REF] Lebrun | The noise-clinic: a universal blind denoising algorithm[END_REF], and our algorithm when estimating either a Gaussian model or a second order one, and PSNR achieved after denoising with the NLF-means algorithm, using the estimated NLF parameters.

Estimator

| [START_REF] Azzari | Gaussian-cauchy mixture modeling for robust signal-dependent noise estimation[END_REF], the PCA method [START_REF]Analysis and extension of the pca method, estimating a noise curve from a single image[END_REF], the percentile method [START_REF] Colom | Analysis and extension of the percentile method, estimating a noise curve from a single image, Image Process[END_REF], Noise Clinic [START_REF] Colom | A non-parametric approach for the estimation of intensity-frequency dependent noise[END_REF][START_REF] Lebrun | The noise-clinic: a universal blind denoising algorithm[END_REF], the VST based method [START_REF] Pyatykh | Image sensor noise parameter estimation by variance stabilization and normality assessment[END_REF] and our algorithm when estimating either an affine model or a second order one, and PSNR achieved after denoising with the NLF-means algorithm, using the estimated NLF parameters. The application of the automatic noise estimation to image denoising shows both the performance of the estimation and the necessity of the knowledge of the noise statistics to perform image processing tasks. Table 7 displays the PSNRs achieved after blind denoising using the NLF-means algorithm proposed above. The images have been corrupted with hybrid noise with different levels. In the first column, the denoising results are obtained using the standard NL-means algorithm (that assumes Gaussian noise). Only the noise variance (i.e. the σ parameter) has been estimated, using the proposed estimation method for a specific Gaussian model. Both the low PSNRs and the artifacts on Figure 9.a) show that the denoising is quite inappropriate; hence the necessity to properly estimate the noise parameters and/or model. The other two columns use the NLF-means algorithm. In the second column, the NLF parameters have been estimated with our proposed method while on the third column we have used to true parameters used to simulate the noise. The superiority of these methods over the standard Gaussian NL-means shows the interest of assuming a relevant noise model, while the feeble difference between the last two columns Computational time (in seconds) for the estimation of the noise parameters on an image of size 540×720 corrupted with Gaussian noise with MAD [START_REF] Donoho | De-noising by soft-thresholding[END_REF], the Gaussian-Cauchy mixture model [START_REF] Azzari | Gaussian-cauchy mixture modeling for robust signal-dependent noise estimation[END_REF], the PCA method [START_REF]Analysis and extension of the pca method, estimating a noise curve from a single image[END_REF], the percentile method [START_REF] Colom | Analysis and extension of the percentile method, estimating a noise curve from a single image, Image Process[END_REF], Noise Clinic [START_REF] Colom | A non-parametric approach for the estimation of intensity-frequency dependent noise[END_REF][START_REF] Lebrun | The noise-clinic: a universal blind denoising algorithm[END_REF], the VST based method [START_REF] Pyatykh | Image sensor noise parameter estimation by variance stabilization and normality assessment[END_REF], and our algorithm.

MRE PSNR Gaussian-Cauchy

Estimator

Time (s) MAD [START_REF] Donoho | De-noising by soft-thresholding[END_REF] 0.58 Gaussian-Cauchy [START_REF] Azzari | Gaussian-cauchy mixture modeling for robust signal-dependent noise estimation[END_REF] 32 PCA [START_REF]Analysis and extension of the pca method, estimating a noise curve from a single image[END_REF] 1 Percentile [START_REF] Colom | Analysis and extension of the percentile method, estimating a noise curve from a single image, Image Process[END_REF] 0.5 Noise Clinic [START_REF] Colom | A non-parametric approach for the estimation of intensity-frequency dependent noise[END_REF] 20 VST [START_REF] Pyatykh | Image sensor noise parameter estimation by variance stabilization and normality assessment[END_REF] 530 Proposed 1.4

illustrates the reliability of our estimation procedure for denoising purposes. Figure 12 displays the result of such denoising on a natural image, corrupted with a hybrid noise of parameters (a, b, c) = (0.0312, 0.625, 100), with an initial PSNR of about 20dB. On Figure 12-a, the standard NL-means have been used, assuming a Gaussian distribution. The single noise parameter σ 2 has been estimated and injected in the NLmeans algorithm. Since the noise model is inappropriate, the denoised image suffers from artifacts on uniform areas and loss of details, on the hat for example. On Figure 12-b, the NLF has been estimated using the proposed method then the estimated parameters (a, b, c) are used for the NLF-means algorithm, while on Figure 12-c the true NLF parameters have been used. Both denoising results offer a satisfying result on which details are preserved, showing that the NLF parameters offer a good model and that the estimation is accurate enough to perform efficient denoising.

Besides, the estimations performed by each estimator in Tables 3 to 6 are re-injected in the NLF-means algorithm to evaluate the noise estimation through the denoising performance, using the PSNR computation. Results show once again that the proposed method offers efficient denoising thanks to the accurate noise level estimation. 6.5. Computational time. We use a fast implementation [START_REF] Knight | A computer method for calculating kendall's tau with ungrouped data[END_REF] that computes Kendall's coefficient in O(N log N ) for a block of size N , resulting in a global complexity in O(N 0 log N ), for an image of size N 0 divided into non-overlapping blocks of size N . For an image of size N 0 = 512 × 512 and blocks of size N = 16 × 16 = 256 pixels, the fast implementation allows to compute Kendall's coefficient for one direction in 0.12 seconds, and the detection on the whole image as well as the computation of the noise statistics on the selected homogeneous areas are performed in about 0.6 seconds. Table 8 displays the computational time for the estimation of the noise parameters on an image of size 540×720 corrupted with Gaussian noise with MAD [START_REF] Donoho | De-noising by soft-thresholding[END_REF], the Gaussian-Cauchy mixture model [START_REF] Azzari | Gaussian-cauchy mixture modeling for robust signal-dependent noise estimation[END_REF], the PCA method [START_REF]Analysis and extension of the pca method, estimating a noise curve from a single image[END_REF], the percentile method [START_REF] Colom | Analysis and extension of the percentile method, estimating a noise curve from a single image, Image Process[END_REF], Noise Clinic [START_REF] Colom | A non-parametric approach for the estimation of intensity-frequency dependent noise[END_REF][START_REF] Lebrun | The noise-clinic: a universal blind denoising algorithm[END_REF], the VST based method [START_REF] Pyatykh | Image sensor noise parameter estimation by variance stabilization and normality assessment[END_REF], and our algorithm, which is quite competitive.

7. Discussion and conclusion. We have developed a fully automatic noise estimation method that relies on the non-parametric detection of homogeneous blocks. These blocks are selected by measuring Kendall's τ coefficient between neighboring pixels, a non-parametric test that requires few hypothesis on the distribution of the samples. In other words, it means that our detector of homogeneous blocks has the exact same performance without extra tuning for any kind of uncorrelated signal-dependent noise including weak/strong Poisson noise, weak/strong Gaussian noise, multiplicative, Rician, hybrid, impulsive, or salt and pepper noise. More precisely, we have shown that the detection errors are controlled with respect to the signal to noise ratio and the number of samples. Next, the noise level function is assumed to follow a second order polynomial model, and is estimated thanks to an estimation method of least absolute deviation. Results on synthetic images show that the proposed method gives an accurate noise level estimation with few assumptions. This NLF estimation can then be used for further image processing applications, for example image denoising, as demonstrated with the use of the NLF estimation in the NL-means algorithm to perform adapted denoising.

Future work might lead to the study of a more general noise model, beyond second order polynomial. First, it would be interesting to encompass spatially varying noise level functions. Indeed, in many image devices the noise level varies with the distance to the sensor for example. Another issue lies in spatially correlated noise. Spatial correlation is not handled by our Kendall-based detection, yet it occurs in many natural image corruptions. Dealing with spatial correlation would extend the noise estimation to artifacts detection, for example compression artifacts. Beyond noise estimation, it might in fact be interesting to extend the proposed approach to other corruptions such as blur detection.

Proof. Using the definition of the expectation and the independence between the two variables, we have:

E[sign(ε i -ε j )] = R 2 sign(ε i -ε j )ϕ(ε i )ϕ(ε j ) dε i dε j = ε i ∈R ε i ε j =-∞ ϕ(ε j ) dε j - +∞ ε j =ε i ϕ(ε j ) dε j ϕ(ε i ) dε i = ε i ∈R ε i ε j =-∞ ϕ(ε j ) dε j - +∞ ε j =ε i ϕ(ε j ) dε j ϕ(ε i ) dε i = ε i ∈R (φ(ε i ) -(1 -φ(ε i ))) ϕ(ε i ) dε i = ε i ∈R (2φ(ε i ) -1) ϕ(ε i ) dε i = 1 4 (2φ(ε i ) -1) 2 +∞ -∞ = 1 4 - 1 4 = 0.
Proposition B.2. Let ε i , ε j and ε k be independent realizations of ε.

(B.3) E[sign(ε i -ε j ) sign(ε i -ε k )] = 1 3 .
Proof. Using once again the independence between the three variables and following the same idea as for Prop. B.1, we have:

E[sign(ε i -ε j )sign(ε i -ε k )] = R 3 sign(ε i -ε j ) sign(ε i -ε k )ϕ(ε i )ϕ(ε j )ϕ(ε k ) dε i dε j dε k = ε i ∈R ε j ∈R sign(ε i -ε j )ϕ(ε j ) dε j ε k ∈R sign(ε i -ε k )ϕ(ε k ) dε k ϕ(ε i ) dε i = ε i ∈R (2φ(ε i ) -1) 2 ϕ(ε i ) dε i = 1 6 (2φ(ε i ) -1) 3 +∞ -∞ = 1 6 + 1 6 = 1 3 .
Proposition B.3. Let ε i and ε j be two independent realizations of ε and a ∈ R.

(B.4) E[sign(ε i -ε j -a)] = 2 R φ(ε i -a)ϕ(ε i ) dε i -1.
Besides, if ε follows a standard distribution N (0, σ 2 ), then:

(B.5) E[sign(ε i -ε j -a)] = 2φ -a √ 2σ -1 = 1 -2φ a √ 2σ .
Proof. In the general case where no particular distribution is assumed for ε, we use the same idea as for the previous propositions.

E[sign(ε

i -ε j -a)] = R 2 sign(ε i -ε j -a)ϕ(ε i )ϕ(ε j ) dε i dε j = ε i ∈R ε i -a ε j =-∞ ϕ(ε j ) dε j - +∞ ε j =ε i -a ϕ(ε j ) dε j ϕ(ε i ) dε i = R (2φ(ε i -a) -1) ϕ(ε i ) dε i = 2 R φ(ε i -a)ϕ(ε i ) dε i -1.
If ε ∼ N (0, σ 2 ), then ε iε j = µ ij follows a Gaussian distribution N (0, 2σ 2 ), with density ϕ 2σ 2 , and

E[sign(ε i -ε j -a)] = E[sign(µ ij -a)] = R sign(µ -a)ϕ 2σ 2 (µ)dµ = R sign( √ 2σ 2 t -a)ϕ(t)dt, with t = µ √ 2σ 2 = - a/ √ 2σ 2 t=-∞ ϕ(t)dt + +∞ t=a √ 2σ 2 ϕ(t)dt = 1 -2φ a √ 2σ = 2φ -a √ 2σ -1.
This last result is confirmed by the following proposition: Proposition B.4. Let ε ∼ N (0, σ 2 ) be a random Gaussian variable. For all a ∈ R, we have:

(B.6) ε φ(ε -a)ϕ(ε) dε = φ -a √ 2σ .
Proof. This proposition can be shown using the property of the error function (defined in Theorem 3.8): for all a, b, c, d ∈ R,

(B.7) erf b -ac √ 1 + 2a 2 d 2 = R erf(ax + b) √ 2πd 2 exp - (x + c) 2 2d 2 dx.
Proof. [Back to the proof of Prop. 3.10] Based on Def. 3.3, and using the expression of X and Y defined in (3.9), we separate the cases according to the location of the indexes in the step signal, and we have:

τ = 1 n(n -1) n i=1 n j=1 sign(X i -X j ) sign(Y i -Y j ) = 1 n(n -1) i, j ≤ n/2 or i, j > n/2, i = j sign(ε i -ε j ) sign(η i -η j ) ( * ) + i ≤ n/2, j > n/2 sign(ε i -ε j -a) sign(η i -η j -a) ( * * ) + i > n/2, j ≤ n/2 sign(ε i -ε j +a) sign(η i -η j + a) ( * * * )
.

Each of the three sums can be evaluated separately using the propositions shown in the first part.

B.9) E[( * * )] = E[sign(ε i -ε j -a)] × E[sign(η i -η j -a)] = 1 -2φ a √ 2σ (
E[τ |H 1 ] = 1 n(n -1) 2 × n 2 2 1 -2φ a √ 2σ 2 = n 2(n -1) 1 -2φ a √ 2σ 2 .
Appendix C. Proof of Proposition 3.11.

Proof. In order to assess the variance of τ under the hypothesis H 1 , we need to evaluate E[τ 2 ]. Using the definition of τ in Def. 3.3, we compute τ 2 then distinguish the cases where X and Y are independent or not. This consists in considering the cases when the indexes can coincide or not.

τ 2 = 1 n 2 (n -1) 2 n i, j = 1 i = j n k, l = 1 k = l = i = j sign(X i -X j ) sign(Y i -Y j ) sign(X k -X l ) sign(Y k -Y l ) ( * ) +4 n i, j = 1 i = j n k = 1 k = i, j sign(X i -X j ) sign(Y i -Y j ) sign(X i -X k ) sign(Y i -Y k ) ( * * ) +2 n i, j = 1 i = j (sign(X i -X j ) sign(Y i -Y j )) 2 =n(n-1)
.

Each of the three sums is then evaluated separately, taking into account the location of the indexes in the step. ( * * ) : i = j = k. In a similar fashion, we distinguish different cases based on the location of the indexes in the step signal, resulting in several sums that we evaluate separately.

( * ) : k = l = i = j. ( * ) = n i, j = 1 i = j n k, l = 1 k = l = i = j sign(X i -X j ) sign(Y i -Y j ) sign(X k -X l ) sign(Y k -Y l ) = n i, j ≤ n/2 or i, j > n/2 n k, l ≤ n/2 or k, l > n/2 sign(ε i -ε j ) sign(η i -η j ) sign(ε k -ε l ) sign(η k -η l ) E=0 + n k ≤ n/2 l > n/2 sign(ε i -ε j ) sign(η i -η j ) sign(ε k -ε l -a) sign(η k -η l -a) E=0 + n k > n/2 l ≤ n/2 sign(ε i -ε j ) sign(η i -η j ) sign(ε k -ε l +a) sign(η k -η l + a)
In order to evaluate the left integral, we use the fact that a > 0 and φ is a nondecreasing function (begin a CDF), so φ(ε i -a) ≤ φ(ε i ). This leads to

ε i φ(ε i -a)φ(ε i )ϕ(ε i ) dε i ≤ ε i φ(ε i ) 2 ϕ(ε i ) dε i = 1 3
, so

E[sign(ε i -ε j ) sign(ε i -ε k -a)] ≤ 4 3 -2φ -a √ 2σ ≤ 1 3 + (1 -2φ -a √ 2σ ) = 1 3 + (2φ a √ 2σ -1).
With the product of the two terms, ( * 3) is also a product of two independent terms, so in a similar way we first evaluate

E[sign(ε i -ε j ) • sign(ε i -ε k +a)] = R 3 sign(ε i -ε j ) sign(ε i -ε k +a)ϕ(ε i )ϕ(ε j )ϕ(ε k ) dε i dε j dε k = ε i ∈R ε j ∈R
sign(ε iε j )ϕ(ε j dε j To evaluate this integral, we use that φ(ε) ≤ 1 for all ε ∈ R:

E[sign(ε i -ε j ) sign(ε i -ε k +a)] ≤ 4 R φ(ε i +a)ϕ(ε i ) dε i -2 R φ(ε i +a)ϕ(ε i ) dε i = 2φ a √ 2σ .
With the products of the two terms, this leads to Using φ(ε i -a) ≤ φ(ε i ), R φ(ε i -a) 2 ϕ(ε i ) dε i ≤ R φ(ε i ) 2 ϕ(ε i ) dε i = 1 3 , so

E[sign(ε i -ε j -a) • sign(ε i -ε k -a)] ≤ 4 3 -4φ -a √ 2σ + 1 ≤ 7 3 -4φ -a √ 2σ = 1 3 + 2(1 -2φ -a √ 2σ ) = 1 3 + 2(2φ a √ 2σ -1).
Taking into account the two multiplied terms, we have

(C.7) E[( * 4)] ≤ 1 3 + 2(2φ a √ 2σ -1) 2 .
Expression ( * 4) refers to a sum where j and k are in the same subset, so the number of terms is n .

Finally, ( * 5) is evaluated as follows:

E[sign(ε i -ε j +a) • sign(ε i -ε k +a)] = R (2φ(ε i +a) -1) 2 ϕ(ε i ) dε i = 4 R φ(ε i +a) 2 ϕ(ε i ) dε i -4 R φ(ε i +a)ϕ(ε i ) dε i + R ϕ(ε i ) dε i .
Since φ(ε) ≤ 1 for all ε ∈ R, φ(ε) 2 ≤ φ(ε), and

E[ sign (ε i -ε j +a) • sign(ε i -ε k +a)] ≤ 4 R φ(ε i +a)ϕ(ε i ) dε i -4 R φ(ε i +a)ϕ(ε i ) dε i + R ϕ(ε i ) dε i =1 ≤ 1.
This leads to (C.8) E[( * 5)] ≤ 1.

As for ( * 4), the number of terms here is n .

Recalling the first result on ( * ) and adding it to ( * * ), we have: In the end, we can derive a bound for the variance of τ under the hypothesis H 1 : This bound effectively shows the decrease rate of τ as the number of samples n increases, hence concluding the proof.

E[τ 2 |H 1 ] ≤ 1 n 2 (n -1) 2 n 2 (n -2) 2 4
Var[τ |H 1 ] = E[τ 2 |H 1 ] -E[τ |H 1 ] 2 ≤ (n -2) 2 4(n -1) 2 1 -2φ a √ 2σ 4 + 1 (n -1)
= µ 4 σ 4 is the kurtosis of the noise distribution and N is the number of pixels in an image block.

If the block size N is large enough, we can use the Law of large numbers to assume that σ 2 follows a standard distribution, i.e.: (F.2) σ 2 ∼ N (σ 2 , Var( σ 2 )) ⇔ σ 2 -σ 2 ∼ N (0, Var( σ 2 ))

Based on the properties that if X follows a normal distribution N (0, σ 2 ), then the central absolute moments are given by the following formula:

(F.3) E [|X| p ] = σ p 2 p 2 Γ p+1 2 √ π .
We can then deduce the expected value and the variance of the estimation error.

E σ 2 -σ 2 σ 2 = Var( σ 2 ) σ 2 × √ 2 √ π
, and (F.4)

Var σ 2 -σ 2 σ 2 = E σ 2 -σ 2 2 σ 4 -E σ 2 -σ 2 σ 2 2 = Var( σ 2 ) σ 4 - Var( σ 2 ) σ 4 × 2 π = Var( σ 2 ) σ 4 × 1 - 2 π . (F.5)
Injecting the variance of the estimator σ 2 of Proposition 4.1, this concludes the proof.

Figure 1 .

 1 Figure 1.Step signal associated to the alternative hypothesis H 1 (see Def. 3.9) and selection of the neighboring sequences.

Figure 3 .

 3 Figure 3. Evolution of the empirical probability of miss-detections as a function of the number of samples n, with probability of false alarm P FA = 0.15. Displayed in a log-log scale, the bound in 1/n is verified.

Figure 4 .

 4 Figure 4. Neighborhood selection for the independence tests.

Figure 5 .

 5 Figure 5. Detection of homogeneous areas in an image corrupted with a hybrid noise as the sum of Gaussian, Poissonian and multiplicative gamma noise whose NLF parameters are (a, b, c) = (0.0312, 0.75, 400), resulting in an initial PSNR of 17.93dB. a) Noisy image (range [0, 255]), b) p-value (range [black = 0, white = 1]) of the associated Kendall's τ coefficient computed within blocks of size N = 16 × 16, and c) selected homogeneous blocks (red) by thresholding the p-value to reach a probability of detection of P D = 1 -P FA = 0.7.

Figure 6 -Figure 6 .

 66 Figure 6. a) Evolution of the expected value and standard deviation of the relative error of the variance estimator as a function of the kurtosis β 2 for a size of block N = 16 × 16, and b) Evolution of the expected value and standard deviation of the relative variance estimation error as a function of the block size for an additive white Gaussian noise.

Figure 7 .

 7 Figure 7. Images extracted from the database of 150 high quality grayscale images of resolution 540×720, taken from a 7 Mpix Canon Powershot S70 or a 8 Mpix Canon IXY Digital 910IS, where the noise can be considered to be negligible. a) Poisson noise b) Hybrid noise

Figure 9 .

 9 Figure 9. Parameter estimation for different types of noise with the LAD estimator. The noise model is supposed to be unknown so the NLF is estimated as a general second order polynomial function.

Figure 10 .

 10 Figure 10. Evolution of the mean relative error as a function of the block width W .

Figure 11 .

 11 Figure 11. Estimation of strong Poisson noise, with parameter b = 1000 and initial PSNR≈2dB.The original and noisy images are displayed on the top line. On the bottom line, the noise parameter are estimated first assuming a Poisson model (on the left), then with the general second order model (on the right). In both cases, the relative error is below 5%, showing the reliability of the estimation.

b

  = 60 b = 100 b = 200 b = 1000 Gaussian-Cauchy [2

Figure 12 .

 12 Figure 12. Denoising of a hybrid noise with true parameters (a, b, c) = (0.0312, 0.625, 100), initial PSNR = 20.34dB. The noisy image is displayed on Fig. 5-a. a) Standard NL-means assuming and estimating Gaussian noise, b) NLF-means with the estimated NLF and c) NLF-means with the true NLF.

(B. 8 )

 8 E[( * )] = 0 using prop. B.1 and the independence of ε i and ε j .

(

a) 2 , 2 ,

 22 using the independence between the two variables ε,η and Proposition B.3. The sum involving ( * * ) occurs when i and j are in the two different halves of the set of samples, so it contains n 2 × n 2 terms.(B.10) E[( * * * )] = E[sign(ε iε j +a)] × E[sign(η iη j + a)] = 2φ a √ 2σ(a) -1 using Prop. B.3 and the independence. As for ( * * ), ( * * * ) occurs when i and j are in different subsets so the sum contains n 2 × n 2 terms. Hence, combining the three terms, (B.11)

( * 3 )√ 2σ 4 , 2 n 2 - 1 2 2 n 2 - 1 2 1 -2φ a √ 2σ 4 .

 3422122114 ≤ n/2 or k, l > n/2 sign(ε iε j -a) sign(η iη ja) sign(ε kε l ) sign(η kη l ) E=0 + n k ≤ n/2 l > n/2 sign(ε iε j -a) sign(η iη ja) sign(ε kε l -a) sign(η kη la) ( * 1) + n k > n/2 l ≤ n/2 sign(ε iε j -a) sign(η iη ja) sign(ε kε l +a) sign(η kη l + a) ≤ n/2 or k, l > n/2 sign(ε iε j +a) sign(η iη j + a) sign(ε kε l ) sign(η kη l ) iε j +a) sign(η iη j + a) sign(ε kε l -a) sign(η kη la) iε j +a) sign(η iη j + a) sign(ε kε l +a) sign(η kη l + a)Since k = l = i = j and using the independence, some of these expectations are null. The other terms are evaluated using Prop. B.3 and the independence, and we have:E[( * 1)] = 1 -2φa and ( * 3) come from sums where i, j and k, l are in the two distinct subsets, so each term occurs n times. ( * 2) appears twice and it comes from sums where i, j and k, l are in the two distinct subsets, so this term occurs 2 × n times. Adding up the expectations, (C.4) E[( * )] = n 2 (n -2) 2 4

( 2 , 2 n 2 n 2 - 1 = n 2 (

 222212 and the term ( * 2) appears 2 × n

4 R

 4 ε k ∈R sign(ε iε k +a)ϕ(ε k ) dε k ϕ(ε i ) dε i = R (2φ(ε i ) -1) (2φ(ε i +a) -1) ϕ(ε i ) dε i using Prop. B.3, = R 2φ(ε i +a) (2φ(ε i ) -1) ϕ(ε i ) dε i since R (2φ(ε i ) -1) ϕ(ε i ) dε i = 0 with Prop. B.1, = φ(ε i +a)φ(ε i )ϕ(ε i ) dε i -2 R φ(ε i +a)ϕ(ε i ) dε i .

( 2 , 2 n 2 n 2 - 1 = 4 R

 222214 and the expression ( * 3) arises 2 × n n 2 (n-2) 4 times.( * 4) is handled in the same manner:E[sign(ε iε j -a) • sign(ε iε k -a)] = R 3 sign(ε iε j -a) sign(ε iε k -a)ϕ(ε i )ϕ(ε j )ϕ(ε k ) dε i dε j dε k = ε i ∈R ε j ∈R sign(ε iε j -a)ϕ(ε j dε j ε k ∈R sign(ε iε k -a)ϕ(ε k ) dε k ϕ(ε i ) dε i = R (2φ(ε i -a) -1) 2 ϕ(ε i ) dε i using Prop. B.3, = φ(ε i -a) 2 ϕ(ε i ) dε i -4 R φ(ε i -a)ϕ(ε i ) dε i

2 n 2 n 2 - 1 =

 2221 n 2 (n-2) 8

2 n 2 n 2 - 1 = n 2 (n- 2 ) 8 .

 2221228 Using the results computed for ( * 1) → ( * 5), the expectation of ( * * ) is bounded by: E[( *

Table 4

 4 Mean relative error (MRE) when estimating Poisson-Gaussian noise with the Gaussian-Cauchy mixture model

	MAD [15] Gaussian-Cauchy [2] PCA [12] Percentile [11] Noise Clinic [13] [27]	a| \ \ \ \ \	| b| \ 0.55 0.05 0.040 |c-c| MRE PSNR c 0.08 0.081 26.33 26.41 0.51 0.513 24.12 \ 0.11 0.105 25.69 \ 0.39 0.389 23.69 \ (24.29)
	Proposed, Gaussian Proposed, hybrid	\ 0.00 1.56 0.12 0.056 0.03 0.030 26.72 \ 26.64

Table 5

 5 Mean relative error (MRE) when estimating multiplicative gamma noise with the Gaussian-Cauchy mixture model[START_REF] Azzari | Gaussian-cauchy mixture modeling for robust signal-dependent noise estimation[END_REF] or and our algorithm, and PSNR achieved after denoising with the NLF-means algorithm, using the estimated NLF parameters.

	[2] Noise Clinic [13] PCA [12] Percentile [11] [28]	\ \ \ \	0.06 0.11 0.040 0.54 2.71 0.88 0.28 0.44 0.12 0.45 0.34 0.390 26.44 26.10 23.67 26.36 (25.18)
	VST [36] Proposed, affine Proposed, hybrid	0.07 0.15 0.034 26.23 0.11 0.25 0.055 26.07 0.01 0.13 0.26 0.054 \ \ 26.02
	Estimator Gaussian-Cauchy [2] 0.80 2.01 120.62 10.521 17.70 |a-a| a | b| MRE PSNR | c|
	Proposed, gamma Proposed, hybrid	0.05 0.06 1.14 \	\ 8.02	0.046 25.69 0.086 25.14

Table 6

 6 Mean relative error (MRE) when estimating hybrid second order noise with the Gaussian-Cauchy mixture model[START_REF] Azzari | Gaussian-cauchy mixture modeling for robust signal-dependent noise estimation[END_REF] or and our algorithm, and PSNR achieved after denoising with the NLF-means algorithm, using the estimated NLF parameters.

	Estimator	|a-a| a	|b-b| b	|c-c| c	MRE PSNR
	Gaussian-Cauchy [2] 0.27 2.53 0.60 0.070 27.81
	Proposed	0.20 1.75 0.36 0.078 27.74

Table 7

 7 Mean PSNRs achieved after denoising a set of noisy images corrupted with hybrid noise, with the standard NL-means assuming a Gaussian model, then with the NLF-means with the estimated NLF and the true NLF.

	Initial PSNR NL-means	NLF-means	NLF-means
		(Gaussian) (estimated NLF) (true NLF)
	19.52	26.26	28.85	28.96
	18.87	25.88	28.48	28.58
	18.04	26.24	27.93	28.09
	17.51	25.86	27.73	27.87

Table 8

 8

http://www.gipsa-lab.grenoble-inp.fr/ ~laurent.condat/imagebase.html

|b-b| b |c-c| c

This study has been carried out with financial support from the French State, managed by the French National Research Agency (ANR) in the frame of the "Investments for the future" Programme IdEx Bordeaux -CPU (ANR-10-IDEX-03-02). C. Sutour would like to thank the DGA and the Aquitaine region for funding her PhD. J.-F. Aujol acknowledges the support of the Institut Universitaire de France.

Appendix A. Proof of Theorem 3.8. Proof. Theorem 3.8 gives the relationship between the level of detection α and the probability of false alarm. In fact, the detection threshold α is actually chosen in order to achieve the desired probability of false alarm associated to condition (C1). Indeed, using the definition of the probability of false alarm and Corollary 3.7, we have:

Now using the choice of α given by α = φ -1 (1 -P FA /2), we get:

Appendix B. Proof of Proposition 3.10. We compute in the following proofs the expectation and variance of Kendall's coefficient under the alternative hypothesis H 1 .

Under H 1 , due to the nature of Gaussian noise, ties are not an issue. We can use the definition given in Def. 3.3 that does not take ties into account. Besides, to simplify the notations, we drop the dependency on (X, Y ) in the notation of τ or z. We also drop the conditional hypothesis in the expressions E[τ |H 1 ] and Var[τ |H 1 ], recalling that the results we demonstrate hereafter apply for the H 1 hypothesis.

Preliminary results. This section shows some basic properties that are useful for the computation of the mean and variance of τ . Let ε be a random variable, with density ϕ. Unless mentioned otherwise, no distribution is supposed for ε. We note ε i , ε j realizations of the random variable ε, and φ the cumulative distributive function (cdf) of ε:

Let ε i and ε j be two independent realizations of ε. Then

.

E[( * 1)] = 1 9 using Prop. B.3, and the expression appears in two sums where i, j, k are in the same subset, resulting in 2

( * 2) is the product of two independent terms, so we first evaluate

with Prop. B.4

Appendix D. Proof of Corollary 3.12. Corollary 3.12 uses the results on the expectation and variance of τ under H 1 to describe the behavior of the associated z-score.

Proof. Under the hypothesis H 1 , due to the nature of Gaussian noise, ties are not an issue. Hence, z is defined as z = . According to Prop. 3.10 and 3.11 that rule the expectation and variance of τ under H 1 , the corollary is then immediate.

Appendix E. Proof of Theorem 3.13. Theorem 3.13 controls the probability of miss-detection, as defined in condition (C2). Proof. Chebyshev's inequality states that if X is a random variable with finite expected value µ and non-zero variance σ 2 , then

Using this inequality for Kendall's z-score, we have

The bound in O 1

n is then ensured by Corr. 3.12. Appendix F. Proof of Proposition 4.2. Proposition 4.2 rules the estimation error between the estimator σ 2 and the true value σ 2 .

Proof. Following Proposition 4.1, the estimator variance is

where β 2 = E X-µ σ