1.2.4. Step 4: Computing the Empirical Variance of the High-Frequency Coefficients. The set of transformed blocks {$\mathbf{D}_0, \dots, \mathbf{D}_m, \dots, \mathbf{D}_{M-1}$} is rewritten with respect to the corresponding value of $\mathbf{V}_m^L$ in ascending order. The block that gives the lowest low-frequency variance will be noted as $\mathbf{D}_{(0)}$, the next with the lowest low-frequency variance as $\mathbf{D}_{(1)}$ and so on. The block with the highest low-frequency variance is $\mathbf{D}_{(M-1)}$. Given the list of sorted blocks {$\mathbf{D}_{(m)}$}, the noise variance estimate associated with the high-frequency coefficient at $(i, j)$ is defined by

$$ \mathbf{V}^H(i, j) = \frac{1}{K} \sum_{k=0}^{K-1} [\mathbf{D}_{(k)}(i, j)]^2, $$

where $i+j \ge T$ and $K = \lfloor pM \rfloor, p < 1$ is the position of the p-quantile in the list {$\mathbf{D}_{(m)}$}$_{m \in [0, M-1]}$. Note that this empirical variance estimate is made with the list of the K transformed blocks whose empirical variance as measured in their low-frequencies is lowest. It is understood that these blocks are likely to contain only noise. Thus their high frequencies are good candidates to estimate the noise. Noise in high and low frequencies is uncorrelated and since most of the energy of the ideal image is concentrated in the low and medium frequency coefficients (because of the sparsity of most natural images), one can assume that $\mathbf{V}^H(i, j)$ gives an accurate estimation of the noise variance. However, if the image is highly textured, those high-frequency coefficients might give a variance that is explained by the textures of the image and not by the noise.

1.2.5. Step 5: Choosing the Best K and Obtaining the Final Noise Estimate. The final noise estimation is given by the median of the variance estimates $\mathbf{V}^H(i, j)$,

$$ \hat{\sigma} := \sqrt{\text{median}_{i,j} \left( \{\mathbf{V}^H(i, j) \mid i + j \ge T\} \right)}. $$

However, the values in the list {$\mathbf{V}^H(i, j)$} depend on the value of the quantile K. Ponomarenko et al. [16] propose to use the following adaptive strategy to find out the best value for K:

(1) Set $K = \sqrt{M}$. The original setting is $K = M/512$, because the algorithm is designed to work with AVIRIS images of size 512 × 677. In order to be able to use any size of image, we propose to set $K = \sqrt{M}$.

(2) Compute an upper bound of noise variance as $A = 1.3\mathbf{V}_{K/2}^L$.

(3) Determine a new $K = m_{\min}$, where $m_{\min}$ is the value of m that minimizes $|A - \mathbf{V}_m^L|$.

(4) Repeat seven times the steps 2 and 3.

(5) Set $A = A/5$.

Nevertheless, we found that fixing directly a small percentile equal to 0.5% of the set of variances gives more accurate and reliable results than the above procedure. This is the only place where our implementation differs from the original algorithm. The complete algorithmic description of the original method is summarized in algorithm 26. The modified version of the algorithm that uses a fixed percentile $p = 0.5\%$ instead of the iterations to find the value of K is given in algorithm 27.

For a review of several noise estimation methods we refer the reader to the *Secrets of image denoising cuisine* [1] and *Estimation of noise in images: an evaluation* [55] articles.
**Algorithm 26** Pseudo-code for the Ponomarenko et al. noise estimation algorithm.

**PONOMARENKO** - Returns the standard deviation of the white Gaussian noise of the input image.

**Input** $\tilde{U}$: noisy image.

**Output** $\hat{\sigma}$: estimated standard deviation of its noise.

1: $w = 8$.

2: $T = 9$.

3: $N_x = \text{width}(\tilde{U})$

4: $N_y = \text{height}(\tilde{U})$

5: $M = (N_x - w + 1)(N_y - w + 1)$

▷ number of (overlapping) blocks in $\tilde{U}$

6: $\mathbf{W} \leftarrow$ all $M$ possible $w \times w$ (overlapping) blocks in $\tilde{U}$.

7: $\mathbf{D} \leftarrow \text{DCT}(\mathbf{W})$.

▷ 2D orthonormal DCT-II of the $w \times w$ blocks in $\mathbf{W}$

8: $\delta[i, j] = 1$ if $(i + j < T) \land (i + j \neq 0)$ else 0, $\forall (i, j) \in [0, w - 1]^2$.

9: $\theta = \sum_{i=0}^{w-1} \sum_{j=0}^{w-1} \delta[i, j]$.

10: $\mathbf{V}_m^L = \frac{1}{\theta} \sum_{i=0}^{w-1} \sum_{j=0}^{w-1} [\mathbf{D}_m(i, j)]^2 \delta[i, j]$

11: $K = \sqrt{M}$.

12: **for** $n = 1...7$ **do**

13: $A = 1.3 \mathbf{V}_{K/2}^L$.

14: $K = \text{argmin}_m (|A - \mathbf{V}_m^L|)$.

15: **end for**

16: $K = K/5$.

17: $\mathbf{I} = \text{sort}_m(\mathbf{V}_m^L)$.

▷ $\mathbf{I}$ contains the sorting indices

18: $\mathbf{V}^H[i, j] = \frac{1}{K} \sum_{k=0}^{K-1} \mathbf{D}_{i[k]}^2(i, j)$

▷ $\mathbf{V}^H[i, j]$ is defined only for those $[i, j]$ such that $i + j \ge T$.

19: $\hat{\sigma} = \sqrt{\text{median}_{i,j} (\mathbf{V}^H(i, j))}$

## 2. Online Demo

**2.1. Example: traffic Image.** The results of this example can be reproduced by adding noise with parameters $A = 0$ and $B = 0.5$ to the traffic image. The rest of the parameters are the default parameters of the demo. Figure 1 shows the input noiseless image traffic before adding signal dependent noise with variance $\sigma^2 = 0.5\mathbf{U}$.

Figure 2 shows the noise estimated for the three first scales of the signal-dependent noise with variance $\sigma^2 = 0.5\mathbf{U}$ added to the traffic image. Because the noise was added to a noise-free image, we can compute the RMSE for the different scales $S_0, S_1$ and $S_2$, and the corresponding errors are 0.15, 0.18 and 0.16, respectively. Note that, as expected, the noise standard deviation is divided by approximately two when down-scaling the image by the same ratio.
**Algorithm 27** Pseudo-code for the Ponomarenko et al. noise estimation algorithm (using a fixed percentile).

**PONOMARENKO WITH FIXED PERCENTILE** - Returns the standard deviation of the white Gaussian noise of the input image.

**Input** $\tilde{U}$: noisy image.

**Output** $\tilde{\sigma}$: estimated standard deviation of its noise.

1: $w = 8$.

2: $T = 9$.

3: $p = 0.005$.

4: $N_x = \text{width}(\tilde{U})$

5: $N_y = \text{height}(\tilde{U})$

6: $M = (N_x - w + 1)(N_y - w + 1) \quad \triangleright \text{number of (overlapping) blocks in } \tilde{U}$

7: $\mathbf{W} \leftarrow \text{all } M \text{ possible } w \times w \text{ (overlapping) blocks in } \tilde{U}$.

8: $\mathbf{D} \leftarrow \text{DCT}(\mathbf{W})$. $\quad \triangleright \text{2D orthonormal DCT-II of the } w \times w \text{ blocks in } \mathbf{W}$

9: $\delta[i, j] = 1 \text{ if } (i + j < T) \land (i + j \neq 0) \text{ else } 0, \forall (i, j) \in [0, w - 1]^2$.

10: $\theta = \sum_{i=0}^{w-1} \sum_{j=0}^{w-1} \delta[i, j]$.

11: $K = pM$. $\quad \triangleright \text{Get the } p\text{-quantile position in the list of variances. Typically } p = 0.005 \Rightarrow \text{ the } 0.5\%-\text{percentile.}$

12: $\mathbf{V}_m^L = \frac{1}{\theta} \sum_{i=0}^{w-1} \sum_{j=0}^{w-1} [\mathbf{D}_m(i, j)]^2 \delta[i, j]$.

13: $\mathbf{I} = \text{sort}_m(\mathbf{V}_m^L)$. $\mathbf{I}$ contains the sorting indices.

14: $\mathbf{V}^H[i, j] = \frac{1}{K} \sum_{k=0}^{K-1} \mathbf{D}_{[k]}^2(i, j)$, $\quad \triangleright \mathbf{V}^H[i, j]$ is defined only for those $[i, j]$ such that $i + j \ge T$

15: $\hat{\sigma} = \sqrt{\text{median}_{i,j}(\mathbf{V}^H(i, j))}$

<figure><img src="image_102.png" /><figcaption>FIGURE 1. Noise free input image traffic before adding noise with variance $\sigma^2 = 0.5U$.</figcaption></figure>
<figure><img src="image_103.png" /><figcaption>FIGURE 2. The noise estimated with the IPOL Ponomarenko et al. demo, for the three first scales of the signal-dependent noise with variance $\sigma^2 = 0.5$ U added to the traffic image. From left to right: scales $S_0$ (original), $S_1$ and $S_2$. Note that the noise standard deviation is approximately divided by two when down-scaling.</figcaption></figure>
# The Percentile method

Given a white Gaussian noise signal $N_\sigma$ on a sampling grid, its variance $\sigma^2$ can be estimated from a small $w \times w$ pixels sample. However, in natural images we observe $\tilde{U} = U + N_\sigma$, the combination of the geometry of the scene that is photographed and the added noise. In this case, estimating directly the standard deviation of the noise from $w \times w$ samples of $\tilde{U}$ is not reliable since the measured standard deviation is not explained just by the noise but also from the geometry of $U$. The Percentile method tries to estimate the standard deviation $\sigma$ from $w \times w$ blocks of a high-passed version of $\tilde{U}$ by a small $p$-percentile of these standard deviations. The idea behind is that edges and textures in a block of the image increase the observed standard deviation but they never make it decrease. Therefore, a small percentile (0.5%, for example) in the list of standard deviations of the blocks is less likely to be affected by the edges and textures than a higher percentile (50%, for example). The 0.5%-percentile is empirically proven to be adequate for most natural, medical, and microscopy images. The Percentile method is adapted to deal with signal-dependent noise, which is realistic with the Poisson noise model obtained by a CCD or CMOS detector in a digital camera.

## 1. Introduction

The Percentile method [1] assumes the signal-dependent noise model. The output of the Percentile method is a noise curve, that is, a function that relates the intensity of the image with a noise standard deviation. The Percentile method tries to estimate the standard deviation $\sigma$ from $w \times w$ blocks of a high-passed version of $\tilde{U}$ by a small $p$-percentile of these standard deviations. The idea behind is that edges and textures in a block of the image increase the observed standard deviation but they never make it decrease. Therefore, a small percentile (0.5%, for example) in the list of standard deviations of the blocks is less likely to be affected by the edges and textures than a higher percentile (50%, for example). However, it might happen that for a certain intensity interval all the samples belong to textures and edges. In that case, the variance measured is explained by the geometry of the image and not by the noise. In order to minimize that effect, the noise curves can be filtered as explained in Section 1.2 of Chapter 9.

For a review of several noise estimation methods we refer the reader to the work of Lebrun et al. [1] and Olsen [55].
## 2. Noise Estimation Method

**2.1. Notation and Terminology.** This chapter prepares the detailed description of the noise estimation algorithm given in Chapter 2.2 by fixing its notation and terminology.

* **U**: the noiseless ideal image.

* **$\tilde{U}$**: the noisy input image.

* $N_x, N_y$: the size of $\tilde{U}$. $N_x$ and $N_y$ are always odd. If they are not, the leftmost column or the bottom row are first removed from $\tilde{U}$.

* $\tilde{U}_c$: the discrete noisy image $\tilde{U}$ after cropping $(s-1)/2$ columns and rows from each of the four sides of $\tilde{U}$ using the function $CROP_{(s-1)/2}$, where $s$ is odd. The details of this function are given in algorithm 28. The size of $\tilde{U}_c$ is therefore $(N_x - (s-1)) \times (N_y - (s-1)) = (N_x - s + 1) \times (N_y - s + 1)$.

* **R**: the operator used to obtain the discrete filter **F** with support $s \times s$, with $s$ odd.

* **F**: the discrete filter with support $s \times s$ used to high-pass the noisy image, obtained from operator **R**.

* **F(x, y)**: the value of the discrete filter **F** at position $(x, y)$, $x \in [0, s-1]$, $y \in [0, s-1]$.

* $\wedge$: logical conjunction (and operator).

* $\vee$: logical disjunction (or operator).

* $\neg$: logical negation (not operator).

* $\tilde{U}^f$: the cropped high-passed version of $\tilde{U}$: $\tilde{U}^f := CROP_{s-1}[(\tilde{U} * F)]$, where $*$ is the discrete convolution operator:

$$ \tilde{U}^f(x, y) := \mathrm{CROP}_{s-1}[(\tilde{U} * F)(x, y)] = \mathrm{CROP}_{s-1}\left[\sum_{i=-\infty}^{\infty} \sum_{j=-\infty}^{\infty} \mathbf{F}_z(i, j) \tilde{U}_z(x - i, y - j)\right] = \\ (30) \qquad = \mathrm{CROP}_{s-1}\left[\sum_{i=0}^{s-1} \sum_{j=0}^{s-1} \mathbf{F}(i, j) \tilde{U}_z(x - i, y - j)\right], $$

$$ \text{where } \mathbf{F}_z(x, y) := \begin{cases} \mathbf{F}(x, y), & x \in [0, s-1] \wedge y \in [0, s-1] \\ 0 & \text{otherwise} \end{cases} \quad \text{extends } \mathbf{F} \text{ with zeros,} $$

$$ \tilde{U}_z(x, y) := \begin{cases} \tilde{U}(x, y), & x \in [0, N_x - 1] \wedge y \in [0, N_y - 1] \\ 0 & \text{otherwise} \end{cases} \quad \text{extends } \tilde{U} \text{ with zeros and the} $$

$CROP_{s-1}$ function removes $s-1$ columns or rows at each of the four boundaries of the filtered image to avoid boundary effects. Since the convolution $(\tilde{U} * F)(x, y)$ is defined for $x \in [0, N_x + s - 2] \wedge y \in [0, N_x + s - 2]$, the cropped image $\tilde{U}^f(x, y)$ is defined for $x \in [0, N_x - s] \wedge y \in [0, N_y - s]$. The Fast Fourier Transform (FFT) algorithm [123] is used to speed-up the computation of the convolution; the details are given in Chapter 2.2.1.

* $\tilde{U}_c(x, y)$: the gray-level intensity of the pixel of $\tilde{U}_c$ at column $x$ and row $y$, $x \in [0, N_x - s]$, $y \in [0, N_y - s]$.
* $\tilde{\mathbf{U}}^f(x, y)$: the gray-level intensity of the pixel of $\tilde{\mathbf{U}}^f$ at column $x$ and row $y$, $x \in [0, N_x-s]$, $y \in [0, N_y-s]$.

* $w$: the side of the overlapping $w \times w$ pixels blocks $\mathbf{W}(x, y)$.

* $M$: the total number of overlapping blocks. $M = (N_x - w - s + 1)(N_y - w - s + 1)$.

* $\mathbf{W}(x, y)$: a $w \times w$ pixels block in $\tilde{\mathbf{U}}_c$, $\mathbf{W}(x, y) = \{\tilde{\mathbf{U}}_c(x + i, y + j) : i \in [0, w-1], j \in [0, w-1]\}$.

* $\mathbf{W}^f(x, y)$: a $w \times w$ pixels block in $\tilde{\mathbf{U}}^f$, $\mathbf{W}^f(x, y) = \{\tilde{\mathbf{U}}^f(x + i, y + j) : i \in [0, w-1], j \in [0, w-1]\}$.

* $\mathbf{W}_m$: a $w \times w$ pixels block in $\tilde{\mathbf{U}}_c$ according to its index $m$ in a list of overlapping blocks $\mathbf{W}_m = \mathbf{W}(x, y)$ where $y = \lfloor \frac{m}{N_x-s+1} \rfloor, x = m - y(N_x - s + 1), m \in [0, M-1]$.

* $\tilde{\mathbf{W}}_m^f$: a $w \times w$ pixels block in $\tilde{\mathbf{U}}^f$ according to its index $m$ in a list of overlapping blocks $\tilde{\mathbf{W}}_m^f = \mathbf{W}^f(x, y)$ where $y = \lfloor \frac{m}{N_x-s+1} \rfloor, x = m - y(N_x - s + 1), m \in [0, M-1]$.

* $\mathbf{V}$: the list of $M$ variances from the blocks $\{\mathbf{W}_m^f\}$. $\mathbf{V}[m]$ corresponds to the variance of the block $\mathbf{W}_m^f$.

* $\mathbf{T}$: the list of $M$ means from the blocks $\{\mathbf{W}_m\}$. $\mathbf{T}[m]$ corresponds to the mean of the block $\mathbf{W}_m$.

* $\mathbf{IV}_n$: the index of the element at the position $n$ in the list $\mathbf{V}$ after sorting the elements of $\mathbf{V}$ in ascending order.

* $p$: the (small) $p$-percentile.

* $\hat{\sigma}^2$: biased variance at the $p$-percentile of $\mathbf{V}$.

* $\tilde{\sigma}^2$: final unbiased variance at the $p$-percentile of $\mathbf{V}$.

## 2.2. The Algorithm.

### 2.2.1. Step 1: Pre-filtering the Input Noisy Image.
In order to get rid of deterministic tendencies due to signal structure, the image is first pre-filtered with a high-pass filter $\mathbf{F}$ implemented as a discrete stencil with support $s \times s$. This stencil corresponds to a discretization of a certain operator $\mathbf{R}$, typically a differential operator or a waveform. Convolving the image with such a filter removes smooth variations inside the blocks, which increases the number of blocks where noise dominates and on which the variance estimate will be reliable. Mastis proposed a similar approach [57], where operator $\mathbf{F}$ writes as a simple subtraction of the average or median to each $7 \times 7$ block. For the Percentile method a filter based on the DCT with support $7 \times 7$ is proposed. Given an $s \times s$ block in the image $\tilde{\mathbf{U}}$ at position $(x, y)$, its orthonormal 2D DCT-II is

$$ (31) \qquad \text{DCT}(\tilde{\mathbf{U}}(x, y))(i, j) := Q_s(i)Q_s(j) \sum_{n_x=0}^{s-1} \sum_{n_y=0}^{s-1} \tilde{\mathbf{U}}(x+n_x, y+n_y) \cos\left[\frac{\pi}{s}\left(n_x + \frac{1}{2}\right)i\right] \cos\left[\frac{\pi}{s}\left(n_y + \frac{1}{2}\right)j\right] $$

with $x \in [0, N_x - s - 1]$, $y \in [0, N_y - s - 1]$, $i \in [0, s - 1]$, $j \in [0, s - 1]$ and $Q_N(k)$ is the normalization factor
<figure><img src="image_104.png" /><figcaption>FIGURE 1. From left to right, the discrete stencils associated $\Delta$, $\Delta\Delta$ and $\Delta\Delta\Delta$ discrete operators, respectively.</figcaption></figure>

$$ (32) \qquad Q_N(k) := \begin{cases} \frac{1}{\sqrt{N}}, & k=0 \\ \sqrt{\frac{2}{N}}, & k \neq 0. \end{cases} $$

Filter $\mathbf{F}$ is made by taking the normalized product of cosines that correspond to the highest frequency $[s-1, s-1]$, that is,

$$ (33) \qquad \mathbf{F}(n_x, n_y) := \frac{2}{s} \cos \left[ \pi \left( n_x + \frac{1}{2} \right) \frac{s-1}{s} \right] \cos \left[ \pi \left( n_y + \frac{1}{2} \right) \frac{s-1}{s} \right], \quad (n_x, n_y) \in [0, s-1]^2. $$

The filter $\mathbf{F}$ presented here was empirically proven to give the best results. However, other typical differential operators can be used, like directional derivatives, the $\Delta$ (Laplace) operator, or its iterations $\Delta\Delta$, $\Delta\Delta\Delta$, all implemented as discrete stencils. Figure 1 shows the discrete stencils associated to these operators.

The filtered image $\tilde{\mathbf{U}}^f$ is obtained by cropping the discrete convolution, $\text{CROP}_{s-1}[(\tilde{\mathbf{U}} * \mathbf{F})(x, y)]$. Note that this cropping operation avoids the boundary effects of the convolution. To speed-up the computation, the Fast Fourier Transform (FFT) algorithm is used:

(1) Consider the signal

$$ (34) \qquad \tilde{\mathbf{U}}_z(x, y) := \begin{cases} \tilde{\mathbf{U}}(x, y), & x \in [0, N_x - s] \land y \in [0, N_y - s] \\ 0 & \text{otherwise.} \end{cases} $$
$\tilde{\mathbf{U}}_z(x, y)$ is defined for $x \in [0, N_x + s - 2] \land y \in [N_y + s - 2]$.

(2) Consider the signal

$$ (35) \qquad \mathbf{F}_z(x, y) := \begin{cases} \mathbf{F}(x, y), & (x, y) \in [0, s-1]^2 \\ 0 & \text{otherwise.} \end{cases} $$
$\mathbf{F}_z(x, y)$ is defined for $x \in [0, N_x + s - 2] \land y \in [N_y + s - 2]$.

(3) Compute the FFT of $\tilde{\mathbf{U}}_z(x, y)$: $\text{FFT}[\tilde{\mathbf{U}}_z(x, y)]$.

(4) Compute the FFT of $\mathbf{F}_z(x, y)$: $\text{FFT}[\mathbf{F}_z(x, y)]$.

(5) Compute the point-wise product of the FFTs:

$$ \text{FFT}[\tilde{\mathbf{U}}_z(x, y)] \times \text{FFT}[\mathbf{F}_z(x, y)]. $$
(6) Compute the inverse FFT of the point-wise product:

$$ \text{FFT}^{-1}(\text{FFT}[\tilde{\mathbf{U}}_z(x, y)] \times \text{FFT}[\mathbf{F}_z(x, y)]). $$

(7) Crop (see algorithm 28) the result to get $\tilde{\mathbf{U}}^f$, the cropped and low-pass filtered version of the noisy input image $\tilde{\mathbf{U}}$:

$$ (36) \qquad \tilde{\mathbf{U}}^f(x, y) := \text{CROP}_{s-1} \left[ \text{FFT}^{-1} \left( \text{FFT}[\tilde{\mathbf{U}}_z(x, y)] \times \text{FFT}[\mathbf{F}_z(x, y)] \right) \right]. $$

**Algorithm 28** Crops the boundary of the input image.

1: **CROP**<sub>b</sub> - Crops the boundary of the input image.

2: **Input** **I**: input image.

3: **Input** $N_x$: width of **I**.

4: **Input** $N_y$: height of **I**.

5: **Input** **b**: width of the boundary that will be removed at each of the four sides of **I**.

6: **Output** **V**: the cropped image

7: **V** = zeros($N_x - 2b, N_y - 2b$)

8: **for** $y = b...N_y - b - 1$ **do**

9:     **for** $x = b...N_x - b - 1$ **do**

10:         **V**[$x - b, y - b$] = **I**[$x, y$]

11:     **end for**

12: **end for**

2.2.2. Step 2: Computing the Sample Variances from $\tilde{\mathbf{U}}^f$ and the Means from $\tilde{\mathbf{U}}$. The size of the filtered noisy image $\tilde{\mathbf{U}}^f$ is $N_x - s + 1 \times N_y - s + 1$ pixels. Therefore, there are $M = (N_x - w - s + 1)(N_y - w - s + 1)$ overlapping blocks of size $w \times w$ pixels. Each overlapping block is referred to as $\{\mathbf{W}_m^f\}$, where $m$ is the index of the block, $m \in \{0, 1, ..., M - 1\}$. Many noise estimation algorithms [4, 60, 6, 5] compute local estimates of the noise variance in small blocks that are used for a final statistical estimation (median, average, percentile, ...).

The empirical variance of each block $\mathbf{W}_m^f$ is computed as

$$ (37) \qquad \text{Var}\mathbf{W}_m^f := \frac{1}{w^2 - 1} \sum_{x=0}^{w-1} \sum_{y=0}^{w-1} [\mathbf{W}_m^f(x, y) - \bar{\mathbf{W}}_m^f]^2 $$

where $\bar{\mathbf{W}}_m^f = \frac{1}{w^2} \sum_{x=0}^{w-1} \sum_{y=0}^{w-1} \mathbf{W}_m^f(x, y)$ is the mean of the block. Let **V** be the list of variances of the blocks $\{\mathbf{W}_m^f\}$, $\mathbf{V}[m] := \text{Var}\mathbf{W}_m^f$. The corresponding means from $\{\mathbf{W}_m\}$ are stored in the list $\mathbf{T}[m] = \bar{\mathbf{W}}_m = \frac{1}{w^2} \sum_{x=0}^{w-1} \sum_{y=0}^{w-1} \mathbf{W}_m(x, y)$, $m \in \{0, 1, ..., M - 1\}$. The mean of each patch will be needed when extending the method to signal-dependent noise (Section 1.1 of Chapter 9), and therefore is stored at this stage of the method.

2.2.3. Step 3: Obtaining a (biased) Noise Variance Estimation from **V** Using **p**. Once the list **V** is built a biased noise variance estimation is obtained by the p-percentile. Set $IVn := \text{SORTED}(\mathbf{V})[n]$ as the function that given a list of real numbers, sorts them in ascending order and returns the sorting indices. For example, if $n = 0$ then $IV0$ is the index of the minimum in
<figure><img src="image_105.png" /><figcaption>FIGURE 2. Left: Noise-free test image computer (see Section 1 of Chapter 13) after adding homoscedastic white Gaussian noise of $\sigma = 10$. Right: pure homoscedastic white Gaussian noise of $\sigma = 10$ image. Both images have the same 704 × 469 size.</figcaption></figure>

the list $V$ and therefore $V[IV0]$ is that minimum. In this step, a biased estimation of the variance is obtained by the small $p$-percentile and is given by

$$ (38) \qquad \hat{\sigma}^2 = V \left[ IV \left[ \frac{p}{100} M + \frac{1}{2} \right] \right], $$

where $M$ is the cardinal of $V$. We obtain a list of variances $IVn$ of $V$ sorted in ascending order, $V[IVn]$ with $n \in [0, M - 1]$. Since in general the variance of the signal (geometry of the image) is higher than the variance of the noise, small percentiles of $V$ are related more to the noise than to the signal. To illustrate it, figure 2 shows the noise-free test image computer after adding homoscedastic white Gaussian noise of $\sigma = 10$, and an image of pure homoscedastic Gaussian noise of $\sigma = 10$. Figure 3 shows the values of $V[IVn]$ depending on $n$ and using 21 × 21 blocks without any filtering. We refer to Section 1 of Chapter 13 for the details about how the noise-free images are obtained. Only for small percentiles the estimation of the variance on the computer image is close to the estimation on pure noise because of the effect of the image edges and textures. The Percentile method tries to avoid the effect of edges and textures by considering the variances under a very low percentile of the block variance histogram.

2.2.4. Step 4: Correcting the Biased Estimation $\hat{\sigma}^2$ to Obtain the Final $\hat{\sigma}^2$ Estimation. When a percentile different from the median ($p = 0.5$) is used, the estimation of the variance obtained is biased by the percentile. Figure 4 shows the values of $V[IVn]$ depending on $n$ using 21 × 21 blocks without any filtering in the image of pure white Gaussian noise (figure 2, right). If the percentile is under the median of the distribution an underestimation of the variance of the noise is obtained; on the other hand, if it is over the median, the result is an overestimation. The median is attained at position $n = \frac{(N_x - w + 1)(N_y - w + 1)}{2} = \frac{(704 - 21 + 1)(469 - 21 + 1)}{2} = 153558$. Since only a small percentile gives a reliable biased estimation of the noise (see Chapter 2.2.3), the estimation is lower than the real average block variance and it has to be corrected in order to get an estimation close to the median. The correction consists of multiplying the biased estimation $\hat{\sigma}$ by a factor. This correction only depends on the percentile, block size and on the chosen operator used to filter the image. Figure 5 shows the correction $C_{w,R}(p) = (\hat{\sigma} - \tilde{\sigma})_{w,R}$ vs. the direct estimation $\hat{\sigma}$ learned
<figure><img src="image_106.png" /><figcaption>FIGURE 3. Values of $V[IVn]$ depending on $n$ using 21 × 21 blocks without any filtering in the noise-free test image computer after adding homoscedastic white Gaussian noise of $\sigma = 10$ (red) and for pure homoscedastic white Gaussian noise of $\sigma = 10$ image (green). Only a part of the values is shown. The estimation in the computer image is only reliable when the percentile of $V$ is small.</figcaption></figure>

with pure-noise images. The correction is linear with the observed $\hat{\sigma}$. As a matter of fact it can be easily proven that there exists a constant $k_{w,R}$ such that $\hat{\sigma} = k_{w,R}C_{w,R}(p) = k_{w,R}(\hat{\sigma} - \tilde{\sigma})$ and then $(k_{w,R} - 1)\hat{\sigma} = k_{w,R}\tilde{\sigma}$. As a consequence,

$$ (39) \qquad \tilde{\sigma} = \frac{k_{w,R} - 1}{k_{w,R}} \hat{\sigma}. $$

Nevertheless, this constant $k_{w,R}$ is not easy to calculate explicitly, but can be learned from simulations. To obtain it, a large image of 4320 × 3232 pixels with all pixels set to zero is used. Homoscedastic Gaussian noise with standard deviation $\sigma$ is simulated and added to this image. Then, the noise is estimated from the noisy image using 200 bins, with a percentile $p \in \{0.01\%, 0.1\%, 0.5\%, 5\%, 10\%, 50\%\}$, a pre-filter operator $R$, which can be chosen between the following: Identity (no filtering), Directional derivative, Laplace, Laplace (2 iterations), Laplace (3 iterations), Laplace (4 iterations), DCT with support $7 \times 7$, DCT with support $5 \times 5$, DCT with support $3 \times 3$ or the filter of the article Fast Noise Variance Estimation [4]. The size of the block is $w \times w$ with $w \in \{3, 7, 8, 21\}$. No curve filtering iterations are used. The averaged estimation along all the bins gives $\hat{\sigma}$.

The Fast Noise Variance Estimation method tries to avoid the influence of image structures (edges and textures) on the image when estimating the noise. To do it, it detects these structures using an operator based on the Laplacian and cancels them. It therefore considers two $3 \times 3$ Laplacian stencils $L_1$ and $L_2$ and computes their difference to obtain the noise estimation operator $L = 2(L_2 - L_1)$.
<figure><img src="image_107.png" /><figcaption>FIGURE 4. values of $V[IVn]$ depending on $n$ using 21 × 21 blocks without any filtering in the image of pure white Gaussian noise (figure 2, right). If the percentile is under the median one obtains an underestimated value; if it is above the median, an overestimation. The median is attained at position $n = \frac{(N_x - w+1)(N_y - w+1)}{2} = \frac{(704-21+1)(469-21+1)}{2} = 153558$.</figcaption></figure>

<figure><img src="image_108.png" /><figcaption>FIGURE 5. Corrections $C_{w,R}(p) = \hat{\sigma} - \tilde{\sigma}$ for several different percentiles, with 21 × 21 blocks and the $R = \Delta\Delta\Delta$ operator, learnt on pure-noise patches.</figcaption></figure>

For example, with $p = 0.5\%$, $w = 21$ and $R = \Delta\Delta\Delta$, this empirical $\frac{k_{w,R}}{k_{w,R}-1}$ factor learned on pure noise is 1.208610869.

The complete algorithmic description of the Percentile method is summarized in algorithm 29.
**Algorithm 29** Percentile noise estimation algorithm.

1: **PERCENTILE** - Returns the standard deviation of the image noise.

2: **Input** $\tilde{U}$: discrete noisy image $U$ after cropping.

3: **Input** $N_x$: (odd) width of the image before cropping.

4: **Input** $N_y$: (odd) height of the image before cropping.

5: **Input** $\mathbf{F}$: discrete filter with support $s \times s$.

6: **Input** $s$: support parameter of the discrete filter $\mathbf{F}$.

7: **Input** $w$: side of the $w \times w$ pixels blocks.

8: **Output** $\tilde{\sigma}$: estimated standard deviation of the image noise.

9: $\tilde{U}^f = \text{CROP}_{s-1}[\tilde{U} * \mathbf{F}]$. ▷ Filter $\tilde{U}$ with the discrete filter $\mathbf{F}$ according to formula (2.1).

10: $M = (N_x - w - s + 1)(N_y - w - s + 1)$. ▷ Number of overlapping blocks.

11: **for** $m = 0... M-1$ **do**

12: $\tilde{W}^f(x, y) = \{\tilde{U}^f(x + i, y + j) : i \in [x, x + w - 1], j \in [y, y + w - 1]\}$

13: $\tilde{W}_m^f = \tilde{W}^f(x, y)$ where $y = \lfloor \frac{m}{N_x - s+1} \rfloor, x = m - y(N_x - s + 1), m \in [0, M-1]$.

14: $\mathbf{V}[m] = \text{Var}\tilde{W}_m^f = \frac{1}{s^2-1}\sum_{x=0}^{s-1}\sum_{y=0}^{s-1}[\tilde{W}_m^f(x, y) - \tilde{W}_m^f]^2$. ▷ Compute sample variances of the blocks.

15: **end for**

16: $\mathbf{IV}n = \text{SORTED}(\mathbf{V})(n) \forall n \in [0, M-1]$. ▷ Get ascending sorting indices.

17: $\hat{\sigma}^2 = \mathbf{V}[\mathbf{IV}[\lfloor \frac{p}{100}M + \frac{1}{2} \rfloor]]$. ▷ Get biased variance estimation.

18: $\hat{\sigma} = \sqrt{\hat{\sigma}^2}$. ▷ Get biased standard deviation.

19: $\tilde{\sigma} = \frac{k_{w,R}^{-1}}{k_{w,R}}\hat{\sigma}$. ▷ Obtain the final unbiased estimation by correcting $\hat{\sigma}$.

### 3. Optimal Parameters

The optimal parameter choice depends on the size of the image. Three possible sizes were fixed: $S_0 = 6M$, $S_1 = \frac{S_0}{4} = \frac{6M}{4}$, $S_2 = \frac{S_1}{4} = \frac{6M}{16}$, $S_3 = \frac{S_2}{4} = \frac{6M}{32}$ and $S_4 = \frac{S_3}{4} = \frac{6M}{128}$, where $M$ stands for megapixels. Table 1 shows the choice of the parameters values according to the size of the image.

<table><caption>TABLE 1. Best percentile $p$, block size $w \times w$, and pre-filter operator $\mathbf{R}$ for the Percentile method according to the size of the image.</caption><thead><tr><th>Image size</th><th>Percentile p</th><th>Block size w × w</th><th>Pre-filter operator R</th></tr></thead><tbody><tr><td>S<sub>0</sub></td><td>0.005</td><td>21 × 21</td><td>DCT supp. 7 × 7</td></tr><tr><td>S<sub>1</sub></td><td>0.005</td><td>15 × 15</td><td>DCT supp. 7 × 7</td></tr><tr><td>S<sub>2</sub></td><td>0.005</td><td>15 × 15</td><td>DCT supp. 7 × 7</td></tr><tr><td>S<sub>3</sub></td><td>0.005</td><td>15 × 15</td><td>DCT supp. 7 × 7</td></tr><tr><td>S<sub>4</sub></td><td>0.005</td><td>5 × 5</td><td>ΔΔΔ</td></tr></tbody></table>

### 4. Online demo

**4.1. Example: traffic Image.** The results of this example can be reproduced by adding noise with parameters $A = 0$ and $B = 0.5$ to the traffic image. The rest of the parameters are
the default parameters of the demo. Figure 1 shows the input noiseless image traffic before adding signal dependent noise with variance $\sigma^2 = 0.5$U.

Figure 6 shows the noise estimated for the three first scales of the signal-dependent noise with variance $\sigma^2 = 0.5$U added to the traffic image. Because the noise was added to a noise-free image, we can compute the RMSE for the different scales $S_0, S_1$ and $S_2$, and the corresponding errors are 0.15, 0.18 and 0.16, respectively. Note that, as expected, the noise standard deviation is divided by approximately two when down-scaling the image by the same ratio.

<figure><img src="image_109.png" /><figcaption>FIGURE 6. The noise estimated with the IPOL Percentile demo, for the three first scales of the signal-dependent noise with variance $\sigma^2 = 0.5$U added to the traffic image. From left to right: scales $S_0$ (original), $S_1$ and $S_2$.</figcaption></figure>
# The PCA method

In the article Image Noise Level Estimation by Principal Component Analysis [3], S. Pyatykh, J. Hesser, and L. Zheng propose a new method to estimate the variance of the noise in an image from the eigenvalues of the covariance matrix of the overlapping blocks of the noisy image. Instead of using all patches of the noisy image, the authors propose an iterative strategy to adaptively chose the optimal set containing the patches with lowest variance. Although the method measures uniform Gaussian noise, it can be easily adapted to deal with signal-dependent noise, which is realistic with the Poisson noise model obtained by a CCD or CMOS detector in a digital camera.

## 1. Noise estimation method

Noise estimation is a necessary preliminary step for most image processing and computer vision algorithms. Most noise estimation algorithms start by applying a high-pass filter to the image or image patches. The energies of these high-passed image patches are used in order to estimate the noise standard deviation by some statistical criterion. Instead of using any fixed high-pass filter, S. Pyatykh, J. Hesser, and L. Zheng [3] decided to apply Principal Component Analysis in order to select an adapted representation for image patches. The energy of the patches with fewer variation are used in order to estimate the noise variance since they are supposed to be less affected by edges and textures [124, 59]. Computing a local estimate of the noise variance from small blocks than then infering the variance of the noise by applying some statistic (median, average, percentile, ...) is technique shared by most noise estimation methods [4, 60, 6, 5]. Unlike other methods that pre-filter the image before extracting noise variance information from the blocks [57], the PCA method measures the variances directly from the noisy patches. For a review of several noise estimation methods we refer the reader to the "Secrets of image denoising cuisine" [1] and "Estimation of noise in images: an evaluation" [55] articles.

The algorithm proposed in [3] applies only for the estimation of white Gaussian noise. However, this assumptions is not realistic because of the way a CCD or a CMOS detector responds to light. The photon emission by a body follows a Poisson distribution which can be approximated by a Gaussian distribution when the number of photons is large enough. But the variance of this Gaussian is signal-dependent. For this reason, the PCA noise estimation algorithm is adapted (see Section 1 of Chapter 9) in order to estimate signal-dependent noise and get a noise curve that assigns an standard deviation of the noise for each gray level intensity of the input image.

1.1. Notation and terminology. This chapter prepares the detailed description of the noise estimation algorithm given in Section 3 by fixing its notation and terminology.
* **x**: the discrete noise-free image of size $S_1 \times S_2$ pixels.

* **n**: a discrete additive Gaussian noise of variance $\sigma^2$ and zero mean of size $S_1 \times S_2$ pixels.

* **y = x+n**: the corrupted image resulting from adding the noise **n** to the noise-free image **x**.

* $M = M_1 \times M_2$: the size in pixels of the overlapping blocks; $M_1 = M_2 = 5 \Rightarrow M = 25$.

* $Q(p)$: the p-quantile of the list of the variances of the overlapping blocks in **y**.

* $N = (S_1 - M_1 + 1) \times (S_2 - M_2 + 1)$: the number of overlapping blocks.

* **X**: matrix of samples of a random vector made from realizations $x_i$, $i \in [1, N]$. It can be written in matrix notation as

$$ \mathbf{X} = \begin{pmatrix} \mathbf{x}_1 \\ \mathbf{x}_2 \\ \vdots \\ \mathbf{x}_N \end{pmatrix} = \begin{pmatrix} x_{11} & x_{12} & \dots & x_{1M} \\ x_{21} & x_{22} & \dots & x_{2M} \\ \vdots & \vdots & \ddots & \vdots \\ x_{N1} & x_{N2} & \dots & x_{NM} \end{pmatrix} $$

where each row of the matrix contains a different blocks, that is, a realization $x_i$ with $i \in [0, N]$.

* **Y**: matrix of samples of a random vector made from realizations $y_i$, $i \in [1, N]$.

* **N**: matrix of samples of a random vector made from realizations $n_i$, $i \in [1, N]$ such that $N \sim N_M(0, \sigma^2 I)$ and $\text{cov}(X, Y) = 0$, where $N$ represents the normal distribution and $I$ is the identity matrix.

* $S_X$ and $S_Y$: the empirical covariance matrices obtained from $X$ and $Y$, respectively.

* $\lambda_{X,1} \ge \lambda_{X,2} \ge \dots \lambda_{X,M}$: the eigenvalues of $S_X$.

* $\lambda_{Y,1} \ge \lambda_{Y,2} \ge \dots \lambda_{Y,M}$: the eigenvalues of $S_Y$.

## 2. Principal components on natural images

The authors claim [3] that information in a natural image patch can be represented by fewer values than the number of pixels in the patch. Formally, the authors claim that any patch lies in a sub-space $V_{M-m} \subset \mathbb{R}^M$ whose dimension $M-m$ is smaller than the number of pixels of the patch, $M$.

**ASSUMPTION 1.** Let $\mathbf{X}$ be the matrix of samples made of $N$ rows, each describing a noise free patch of size $M$. Then we assume that all realizations $x_i$ lie in a sub-space $V_{M-m} \subset \mathbb{R}^M$ whose dimension $M-m$ is smaller than the dimension $M$ of the $x_i$. Therefore, $x_i$ have zero variance along any direction orthogonal to $V_{M-m}$

$$ (40) \qquad \lambda_{\mathbf{X},M-m-1} = \dots = \lambda_{\mathbf{X},M} = 0. $$

By assuming the previous result the authors show that for noisy patch observations **Y** the following results holds.

**THEOREM 2.** Let $\mathbf{Y} = \mathbf{X} + \mathbf{n}$, where $\mathbf{n}$ is uniform Gaussian noise of variance $\sigma^2$ then

$$ (41) \qquad \mathbb{E}|\lambda_{\mathbf{Y},i} - \sigma^2| \le T\sigma^2/\sqrt{N}, \quad i = M-m+1, \dots, M $$
where $\lambda_{\mathbf{Y},i}$ denote the eigenvalues of the covariance matrix $S_{\mathbf{Y}}$, $N$ the number of rows in $\mathbf{Y}$ and $T$ is a positive constant.

If the number of blocks is large enough, the variance of noise can be estimated from $\lambda_{\mathbf{Y},M}$. Indeed, taking $N \to \infty$ and the $M$-th eigenvector $\lambda_{\mathbf{Y},M}$ leads to

$$ (42) \qquad \lim_{N \to \infty} \mathbb{E} |\lambda_{\mathbf{Y},M} - \sigma^2| = 0. $$

This result provides the way to estimate the variance of the noise from the $M$-th eigenvector $\lambda_{\mathbf{Y},M}$ when the number of blocks $N$ is large enough. In order to check that this approximation is correct the authors use inequality (41) with $i = M - m + 1$ and approximating $\sigma^2$ by $\lambda_M$

$$ (43) \qquad \lambda_{\mathbf{Y},M-m+1} - \lambda_{\mathbf{Y},M} < T\sigma^2/\sqrt{N}. $$

The authors fix the parameter $m = 7$ for blocks of $M_1 \times M_2 = 5 \times 5$ pixels and $T = 49$.

In order to impose patches to actually belong to a $V_{M-m}$ subspace, the authors choose a subset of patches defined as those with smallest variance. The authors choose to use the variance of the blocks in $\mathbf{y}$ as a measure of the distance of a block $\mathbf{y}_i$ to $V_{M-m}, i \in [1, N]$. A simple strategy to select the subset of blocks for which equation (43) is verified consists in taking a small enough $p$-quantile $Q(p)$ of the list of variances of the noisy blocks and to set

$$ (44) \qquad B(p) := \{\mathbf{y}_i : s^2(\mathbf{y}_i) \le Q(p), i \in [1, N]\}. $$

In practice, the authors propose an iterative strategy which reduces the quantile $p$ while equation (43) is not satisfied.

### 3. Algorithm

The noise variance estimation algorithm (Algorithm 30) first obtains an upper bound of the variance of the noise. This upper bound is also used as initialization for the iterative procedure which refines the estimated variance, Function `GetNextEstimate` (in Algorithm 30). This function is iterated until the absolute difference between the new variance and the previous one is negligible.

The function `GetNextEstimate` uses a noise variance initialization, $\sigma_{est}$, in order to test inequality (43). The eigenvalues of the set of patches are computed and checked to see if they satisfy Equation (43) with $\sigma = \sigma_{est}$. The number of patches used for this estimation is reduced by decreasing $p$ for the set of blocks $B(p)$ while inequality (43) does not hold. The iteration begins setting $p = 1$ and it is decremented by $\Delta p = 0.05$ after each iteration.

### 4. Optimizing the PCA computation

Computing the eigenvectors and the corresponding eigenvalues of the sample covariance matrix is a time consuming operation which is called several times in Function `GetNextEstimate`. The PCA is computed on the set $B(p)$ for different values of $p$. It is possible to use the nested structure of these sets in order to avoid recomputing the covariance matrix at each step.
**Algorithm 30 PCA noise estimation.**

1: **PCA noise estimation** - Compute noise variance.

**Input** $\mathbf{y}$: list of blocks of the noisy image.

**Output** $\sigma_{est}^2$: variance of the noise.

2: $i_{max} = 10$

3: $\varepsilon = 0.001$

4: $\sigma_{ub}^2 = 3.1 Q(0.0005, \mathbf{y})$ ▷ Function $Q(p, \mathbf{y})$ gives the $p$-quantile of $\mathbf{y}$

5: $\sigma_{est}^2 = \sigma_{ub}^2$

6: **for** $i = 1 \dots i_{max}$ **do**

7: $\sigma_{next}^2 = \text{GetNextEstimate}(\mathbf{y}, \sigma_{est}^2, \sigma_{ub}^2)$

8: **if** $|\sigma_{est}^2 - \sigma_{next}^2| < \varepsilon$ **then**
   **return** $\sigma_{est}^2$

9: **end if**

10: $\sigma_{est}^2 = \sigma_{next}^2$

11: **end for**
   **return** $\sigma_{est}^2$

1: **GetNextEstimate** - Refine current variance estimation

**Input** $\mathbf{y}$: list of blocks of the noisy image.

**Input** $\sigma_{est}^2$: current value of the estimated noise variance.

**Input** $\sigma_{ub}^2$: upper bound of the noise variance.

**Output** $\sigma_{next}^2$: next approximation of the estimated noise variance.

2: $p_{min} = 0.06$

3: $\Delta p = 0.05$

4: $M = M_1 \times M_2 = 5 \times 5 = 25$

5: $m = 7$

6: $T = 49$

7: $p = 1$

8: $\sigma_{next}^2 = 0$

9: **while** $p \ge p_{min}$ **do**

10: $\lambda_{\mathbf{y},i} = \text{ApplyPCA}(B(p))$ ▷ $\lambda_{\mathbf{y},i}$ is the $i$-th eigenvalue of $S_{\mathbf{y}}$.

11: $\sigma_{next}^2 = \lambda_{\mathbf{y},M}$

12: **if** $(\lambda_{\mathbf{y},M-m+1} - \lambda_{\mathbf{y},M} < \frac{T\sigma_{est}^2}{\sqrt{|B(p)|}}) \land (\sigma_{next}^2 \le \sigma_{ub}^2)$ **then**
   **return** $\sigma_{next}^2$

13: **end if**

14: $p = p - \Delta p$

15: **end while**
   **return** $\sigma_{next}^2$

The sample covariance matrix can be written as

$$ (45) \qquad \frac{1}{|B(p)|-1} \left( \sum_{\mathbf{y}_i \in B(p)} \mathbf{y}_i \mathbf{y}_i^T - \frac{1}{|B(p)|} \sum_{\mathbf{y}_i \in B(p)} \mathbf{y}_i \sum_{\mathbf{y}_i \in B(p)} \mathbf{y}_i^T \right). $$
<figure><img src="image_110.png" /><figcaption>FIGURE 1. Raw image IMG_0177.</figcaption></figure>

The number of operations needed to compute directly this matrix is proportional to $|B(p)| M^2$.

Since the sets $B(p)$

$$B(p) = \{\mathbf{y}_i : s^2(\mathbf{y}_i) \le Q(p), i \in [1, N]\}.$$

satisfy

$$B(1) \supset B(1 - \Delta p) \supset B(1 - 2\Delta p) \supset B(1 - 3\Delta p) \supset \dots,$$

then the set $B(1 - j\Delta p)$ can be written as $B(1 - j\Delta p) = B(1 - (j + 1)\Delta p) \cup \mathbf{Y}_j$ being the set $\mathbf{Y}_j := \{\mathbf{y}_i : Q(1 - (j + 1)\Delta p) < s^2(\mathbf{y}_i) \le Q(1 - j\Delta p)\}$.

Since for disjoint sets $X_1$ and $X_2$ one has $C_{X_1 \cup X_2} = C_{X_1} + C_{X_2}$, then

$$C_{B(1-j\Delta p)} = C_{B(1-(j+1)\Delta p)} + C_Y.$$

Let $C_X := \sum_{\mathbf{y}_i \in X} \mathbf{y}_i \mathbf{y}_i^T$ and $c_X := \sum_{\mathbf{y}_i \in X} \mathbf{y}_i$. With this notation, Expression (45) can be written as

$$(46) \qquad \frac{1}{|B(p) - 1|} \left( C_{B(p)} - \frac{1}{|B(p)|} c_{B(p)} c_{B(p)}^T \right).$$

This strategy permits to use the covariance matrices of previous iteration, but not the eigen-value decomposition. With Expression (46) the number of operations needed is proportional to $M^2$.

## 5. Online demo

**5.1. Example: raw image IMG_0177 (ISO 1250, t=1/30s).** The results of this example can be reproduced by estimating the noise of the raw image without adding any noise (set $A = B = 0$) and choosing the raw image IMG_0177 (Figure 1) The rest of the parameters are the default parameters of the demo. Figure 1 shows the input noisy image.

Figure 2 shows the noise estimated for the three first scales of the signal-dependent noise with variance $\sigma^2 = 0.5$U added to the IMG_0177 image. From left to right: scales $S_0$ (original), $S_1$ and $S_2$. Note that, as expected, the noise standard deviation is divided by approximately two when down-scaling the image by the same ratio.
<figure><img src="image_111.png" /><figcaption>FIGURE 2. The noise estimated for the two first scales of the noisy raw image IMG_0177 (Figure 1) as shown in the IPOL PCA demo. Left: scale $S_0$. Right: scale $S_1$.</figcaption></figure>

## 6. Appendix: proof of Theorem 2

The PCA noise estimation method is based on the fact that the noise variance can be estimated from $\lambda_{Y,M}$ when the condition of Theorem 2 holds and converge is reached (Algorithm 15). To justify completely this result, the proof of Theorem 2 the authors used is reproduced here. To proof Theorem 2 the Gerschgorin circle theorem has to be reviewed, since a lemma from it will be used in the demonstration.

**THEOREM 3. (Gerschgorin) [125]** For a matrix $A \in \mathbb{C}^{M \times M}$ with entries $a_{ij}$ let

$$a_i := \sum_{j \neq i} |a_{i,j}|$$

and

$$G_i(A) := \{z \in \mathbb{C} : |z - a_{ii}| \le a_i\}$$

Then all eigenvalues of $A$ lie in $\bigcup_{i=1}^{M} G_i(A)$. Moreover, if $m$ of the Gerschgorin disks $G_i(A)$ are isolated from the other $M-m$, then there are exactly $m$ eigenvalues of $A$ in their union.

**LEMMA 1.** Let $A \in \mathbb{C}^{M \times M}$ be a random Hermitian matrix, $\lambda$ an eigenvalue of $A$ with multiplicity $m$ and normalized eigenvectors $v_1, \dots, v_m$ and $\tilde{A} = A + B$ a perturbed matrix. If there is a constant $\delta > 0$ such that the distance between $\lambda$ and the other eigenvalues of $A$ is always greater than $\delta$, then there are exactly $m$ eigenvalues of $\tilde{A}$, which satisfy

$$|\tilde{\lambda}_k - \lambda| \le \max_{i=1,\dots,m} \sum_{j=1}^{m} |v_i^T B v_j| + O(\|B\|^2)$$

**Proof:** let $\Lambda = Q^T A Q$ be the eigen-decomposition of $A$ and $q_i$ the columns of $Q$. Without loss of generality

$$\Lambda = \operatorname{diag}\{\underbrace{\lambda, \dots, \lambda}_{m}, \underbrace{\lambda', \dots, \lambda'}_{M-m}\}$$
where $\lambda'$ is an eigenvalue of A other than $\lambda$. Then the vectors $\mathbf{v}_1, \dots, \mathbf{v}_m$ form the first $m$ columns of $Q$:

$$ \mathbf{q}_i = \mathbf{v}_i, \quad i = 1, \dots, m. $$

Let us consider the matrix $\tilde{\Lambda}_a = D_a^{-1}Q^T \tilde{A} Q D_a$, where

$$ D_a = \operatorname{diag}\{\underbrace{1, \dots, 1}_{m}, \underbrace{a, \dots, a}_{M-m}\}, \quad a > 0. $$

The matrices $\tilde{A}$ and $\tilde{\Lambda}_a$ are similar and, hence, have the same eigenvalues.

$$ \begin{aligned} \tilde{\Lambda}_a &= D_a^{-1} Q^T (A+B) Q D_a = \\ &= \Lambda + D_a^{-1} Q^T B Q D_a = \\ &= \Lambda + \begin{bmatrix} b_{11} & \cdots & b_{1m} & ab_{1m+1} & \cdots & ab_{1M} \\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ b_{m1} & \cdots & b_{mm} & ab_{m+1} & \cdots & ab_{mM} \\ a^{-1}b_{m+1} & \cdots & a^{-1}b_{m+1m} & b_{m+1m+1} & \cdots & b_{m+1M} \\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ a^{-1}b_{M1} & \cdots & a^{-1}b_{Mm} & b_{Mm+1} & \cdots & b_{MM} \end{bmatrix} \end{aligned} $$

where $b_{ij} = \mathbf{q}_i^T B \mathbf{q}_j$. Note that $|b_{ij}| \le \|\mathbf{q}_i\| \|B\| \|\mathbf{q}_j\| \le \|B\|$. The first $m$ Gerschgorin disks of $\tilde{\Lambda}_a$ have centers $\lambda + b_{ii}$ and radii bounded by

$$ \sum_{j=1}^{m} |b_{ij}| - |b_{ii}| + (M-m)a \|B\|. $$

Hence, these disks lie entirely in the circle

$$ C = \{z \in \mathbb{C} : |z - \lambda| \le p + (M-m)a \|B\|\} $$

where $p = \max_{i=1,\dots,m} \sum_{j=1}^{m} |b_{ij}|$. The other $M-m$ Gerschgorin disks have centers $\lambda' + b_{ii}$ and radii bounded by

$$ ma^{-1} \|B\| + (M - m - 1) \|B\|. $$

Therefore, if

$$ (47) \qquad p + (M - m)a \|B\| + ma^{-1} \|B\| + (M - m) \|B\| < \delta, $$

the circle $C$ will be disjoint from the last $M-m$ Gerschgorin disks. Let $B$ be so small that

$$ p + (M - m) \|B\| < \frac{\delta}{2}. $$

In this case, (47) is satisfied if

$$ (48) \qquad (M - m)a \|B\| + ma^{-1} \|B\| < \frac{\delta}{2}, $$

and (48) is satisfied if

$$ (49) \qquad Ma \|B\| + Ma^{-1} \|B\| < \frac{\delta}{2}, $$
If $B$ is small enough so that $4M \|B\| < \delta$, (49) holds for

$$a = \frac{4M \|B\|}{\delta}.$$

Then, the circle C:

(1) has the radius bounded by

$$p + 4M(M - m) \|B\|^2 / \delta = p + O(\|B\|^2)$$

(2) contains the first $m$ Gerschgorin disks of $\tilde{\Lambda}_a$.

(3) is disjoint from the last $M - m$ Gerschgorin disks of $\tilde{\Lambda}_a$.

From Gerschgorin's theorem, C contains exactly $m$ eigenvalues of $\tilde{\Lambda}_a$ and, therefore, exactly $m$ eigenvalues of $\tilde{A}$.

**Proof of Theorem 2:** Consider the matrices

$$A := S_X + \frac{\sigma^2 N}{N - 1} I,$$

$$B := S_N - \frac{\sigma^2 N}{N - 1} I + S_{XN} + S_{NX}.$$

The matrix $S_Y$ can be represented in the following form: $S_Y = S_X + S_N + S_{XN} + S_{NX} = A + B$. Since

$$A \mathbf{v}_{\mathbf{X},i} = S_{\mathbf{X}} \mathbf{v}_{\mathbf{X},i} + \frac{\sigma^2 N}{N-1} \mathbf{v}_{\mathbf{X},i} = \left( \lambda_{\mathbf{X},i} + \frac{\sigma^2 N}{N-1} \right) \mathbf{v}_{\mathbf{X},i},$$

each eigenvalue of $A$ equals the sum of an eigenvalue of $S_X$ and $\frac{\sigma^2 N}{N-1}$, and the eigenvectors of $S_X$ and $A$ are the same. Under assumption 1, the last $m$ eigenvalues of $S_X$ are zeros, therefore the last $m$ eigenvalues of $A$ equal $\frac{\sigma^2 N}{N-1}$. Let $J = \{M - m + 1, \dots, M\}$ be the set of indices of zero eigenvalues of $S_X$. Using Lemma 1, for $j \in J$

$$(50) \qquad \left| \lambda_{\mathbf{Y},k} - \frac{\sigma^2 N}{N-1} \right| \le \max_{i \in J} \sum_{j \in J} \left| \mathbf{v}_{\mathbf{X},i}^T B \mathbf{v}_{\mathbf{X},j} \right| + O(\|B\|^2).$$

Consider the first part on the right side of (50). Denoting the sample covariance by $q$, for $i, j \in J$ we have

$$|\mathbf{v}_{\mathbf{X},i}^T S_{\mathbf{X}N} \mathbf{v}_{\mathbf{X},j}| = |q(\mathbf{v}_{\mathbf{X},i}^T \mathbf{X}, \mathbf{v}_{\mathbf{X},j}^T \mathbf{N})| \le \sqrt{s^2(\mathbf{v}_{\mathbf{X},i}^T \mathbf{X}) s^2(\mathbf{v}_{\mathbf{X},j}^T \mathbf{N})} = \sqrt{\lambda_{\mathbf{X},i} s^2(\mathbf{v}_{\mathbf{X},j}^T \mathbf{N})} = 0$$

and, similarly,

$$\mathbf{v}_{\mathbf{X},i}^T S_{\mathbf{N}X} \mathbf{v}_{\mathbf{X},j} = 0.$$

Hence, for $i \in J$

$$\begin{align*} \sum_{j \in J} |\mathbf{v}_{\mathbf{X},i}^T B \mathbf{v}_{\mathbf{X},j}| &= \sum_{j \in J} |\mathbf{v}_{\mathbf{X},i}^T (S_N - \frac{\sigma^2 N}{N-1} I) \mathbf{v}_{\mathbf{X},j}| \le \\ &\le \sum_{j \in J} \|\mathbf{v}_{\mathbf{X},i}\| \|S_N - \frac{\sigma^2 N}{N-1} I\| \|\mathbf{v}_{\mathbf{X},j}\| = \\ &= m \|S_N - \frac{\sigma^2 N}{N-1} I\| \end{align*}$$
and

$$ (51) \qquad \max_{i \in J} \sum_{j \in J} |\mathbf{v}_{\mathbf{X},i}^T B \mathbf{v}_{\mathbf{X},j}| \le m \left\| S_N - \frac{\sigma^2 N}{N-1} I \right\| \le m \sqrt{\sum_{i,j=1}^M \left( n_{ij} - \frac{\delta_{ij} \sigma^2 N}{N-1} \right)^2}, $$

where $n_{ij}$ are the entries of $S_N$ and $\delta_{ij}$ is the Kronecker delta. Since $N \sim N_M(0, \sigma^2 I)$, the matrix $(N-1)S_N$ has the Wishart distribution $W_M(\sigma^2 I, M)$. Therefore,

$$ \mathbb{E}(n_{ii}) = \frac{N\sigma^2}{N-1}, \quad \mathrm{Var}(n_{ii}) = \frac{2\sigma^4 N}{(N-1)^2}, $$

$$ (52) \qquad \mathbb{E}(n_{ij}) = 0, \quad \mathrm{Var}(n_{ij}) = \frac{\sigma^4 N}{(N-1)^2}, \quad i \neq j $$

and

$$ (53) \qquad \left[ \mathbb{E} \left( n_{ij} - \frac{\delta_{ij} \sigma^2 N}{N-1} \right) \right]^2 = \mathrm{Var} (n_{ij}) \le \frac{2 \sigma^4 N}{(N-1)^2}. $$

Since $\mathbb{E}(\sqrt{X}) \le \sqrt{\mathbb{E}(X)}$, (51) and (53) can be combined to give

$$ (54) \qquad \mathbb{E} \left( \max_{i \in J} \sum_{j \in J} |\mathbf{v}_{\mathbf{X},i}^T B \mathbf{v}_{\mathbf{X},j}| \right) \le m \sqrt{\sum_{i,j=1}^M \left( n_{ij} - \frac{\delta_{ij} \sigma^2 N}{N-1} \right)^2} \le \frac{m M \sigma^2 \sqrt{2N}}{N-1} = O(\sigma^2/\sqrt{N}). $$

For the second part on the right side of (50), let us construct an upper bound for $\mathbb{E}(\|B\|^2)$ now. Let $b_{ij}$ be the entries of $B$ and $c_{ij}$ the entries of $S_{\mathbf{X}N} + S_{\mathbf{N}X}$. Since

$$ \mathbb{E}(S_{\mathbf{X}N} + S_{\mathbf{N}X}) = \mathrm{cov}(\mathbf{X}, \mathbf{N}) + \mathrm{cov}(\mathbf{N}, \mathbf{X}) = 0, $$

$\mathbb{E}(c_{ij}) = 0$. Additionally, $\mathrm{Var}(c_{ij}) = \sigma^2 \frac{\mathrm{Var}(X)}{N-1}$. Combining it with (52), it gives

$$ \mathbb{E}(b_{ij}) = 0 $$

and

$$ \begin{aligned} \mathrm{Var}(b_{ij}) &= \mathrm{Var}(n_{ij}) + \mathrm{Var}(c_{ij}) + 2\mathrm{cov}(n_{ij}, c_{ij}) \le \\ &\le (\sqrt{\mathrm{Var}(n_{ij})} + \sqrt{\mathrm{Var}(c_{ij})})^2 \le \\ &\le \left( \frac{\sigma^2 \sqrt{2N}}{N-1} + \frac{\sigma \sqrt{\mathrm{Var}(\mathbf{X})}}{\sqrt{N-1}} \right)^2 = \\ &= O\left( \frac{(\sigma^2 + \sigma \sqrt{\mathrm{Var}(\mathbf{X})})^2}{N} \right). \end{aligned} $$

Therefore,

$$ (55) \qquad \mathbb{E}(\|B\|^2) \le \mathbb{E} \left( \sum_{i,j=1}^M b_{ij}^2 \right) = \sum_{i,j=1}^M \mathrm{Var}(b_{ij}) = O \left( \frac{(\sigma^2 + \sigma \sqrt{\mathrm{Var}(\mathbf{X})})^2}{N} \right). $$
Combining (50), (54), and (55) gives the final result. Since $N \to \infty$, $\frac{1}{N}$ is infinitesimal compared to $\frac{1}{\sqrt{N}}$,

$$ \mathbb{E}\left(\left|\lambda_{\mathbf{Y},i} - \frac{\sigma^2 N}{N-1}\right|\right) = O\left(\frac{\sigma^2}{N}\right) + O\left(\frac{\left(\sigma^2 + \sigma\sqrt{\mathrm{Var}(\mathbf{X})}\right)^2}{N}\right) = O\left(\frac{\sigma^2}{\sqrt{N}}\right) \square $$
# Evaluation of the adapted methods

In Chapters 10, 11, and 12, the Ponomarenko et al., Percentile and PCA signal-dependent noise estimation methods were described in detail. In this chapter we evaluate them using three different tests:

* A test with simulated white Gaussian noise using the noise-free images in Figure 1 (Section 1). The aim of this test is to verify that the methods give accurate results for the noise is homoscedastic (fixed variance, not depending on the intensity). As explained in Chapter 9, it is needed that the noise estimator is accurate before adapting it to signal-dependent noise, since the number of samples will be shared by total number of bins. We consider seven bins to classify the blocks according to their means (see Section 1.1 of Chapter 9).

* A test with raw images (Section 2), which compares the noise curve obtained by the methods to the ground-truth noise curve obtained by the camera that took the photograph, under the same ISO speed (see Chapter 2).

* A test on multiscale coherence (Section 3). The standard deviation of a Gaussian white noise is divided by two when the image is down-scaled. By down-scaling the image we mean a sub-sampling of the image where each block of four pixels is substituted by their mean. This test checks if the measured noise is divided by two at each image down-scaling.

## 1. Evaluation with simulated white Gaussian noise

In this test, white Gaussian noise was simulated and then added to a set of ten noise-free images. Since the noise is perfectly known a priori, it can be used as a ground truth. One can therefore compute the RMSE of the standard deviation estimations. Fig. 1 shows a set of 704 × 469 pixels noise-free images that were used in this test. In order to get the noise-free images, we have applied the following procedure. The pictures were taken with a Canon EOS 30D reflex camera of scenes under good lighting conditions and with a low ISO level. To reduce further the noise level, the average of each block of 5 × 5 pixels was computed, reducing therefore the noise by a factor of 5. Since the images are RGB, the mean of the three channels was computed, reducing the noise by a further $\sqrt{3}$ factor. Therefore the noise was reduced by a $5\sqrt{3} \approx 8.66$ factor. Finally, the images, which already had a good SNR before they were processed, can be considered noise-free.

To measure the error made when estimating the standard deviation $\sigma$ of the simulated noise in the bin $b$ in the image $i$, the RMSE along all the bins was used. This RMSE is denoted by $E_{i,\sigma}^{(1)}$
<figure><img src="image_112.png" /><figcaption>FIGURE 1. Set of noise-free images used to test the noise estimation algorithm with white Gaussian noise. From left to right and from top to bottom: bag, building1, computer, dice, flowers2, hose, lawn, leaves, stairs and traffic. Each image is 704 × 469 pixels.</figcaption></figure>

and it is defined by

$$ (56) \qquad E_{i,\sigma}^{(1)} := \sqrt{\frac{1}{|B|} \sum_{b=1}^{|B|} |\hat{\sigma}_{i,b} - \sigma|^2}, $$

where $|I|$ is the number of images, $i$ is the image index ($1 \le i \le |I|$), $|B|$ is the number of bins, $b$ is the index of the bin ($1 \le b \le |B|$), $\sigma$ is the standard deviation of the simulated noise and $\hat{\sigma}_{i,b}$ is the estimated noise for the image $i$ at the bin $b$. Tables 1, 2, and 3 show the obtained $E_{i,\sigma}^{(1)}$ for each image $i$ and each $\sigma$ of the simulated noise, for the Ponomarenko, Percentile, and PCA methods, namely. For the Ponomarenko et. al and the Percentile methods seven bins are used. For the PCA method three bins are used, since it requires more samples/bin. A new image is added to the set of noise-free images, the flat image, which is a constant image where all pixels have the value 127, useful to measure the performance of the noise estimators with an image of pure homogeneous Gaussian noise.

Note that the PCA method is able to give an accurate estimation of the noise even when $\sigma = 1$ when the image has large flat zones (building1, computer, dice, flowers2, traffic), but fails to give a good estimation for very textured images (bag, hose, leaves, lawn). For very high noise levels ($\sigma = 50, \sigma = 80$) the estimation is inaccurate. Surprisingly, the PCA method behaves better for small levels of noise.

For the Percentile method, it is apparent that the highly textured images create a significant error, particularly when little noise was added. Estimates of noise below $\sigma = 2$ are therefore obviously clearly unreliable. All in all, the estimate is nevertheless quite reliable for values $\sigma > 5$.
<table><caption>TABLE 1. This table shows, for the Ponomarenko method, the $E_{i,σ}^{(1)}$ RMSE after adding simulated noise to the set of noise-free images (Figure 1) with several values of standard deviation $σ$. The last row is the $E_{σ}^{(2)}$ RMSE using the estimated $\hat{\sigma}_{i,b}$ of all the images. The percentile $p = 0.005$ and seven bins are used.</caption><thead><tr><th>Image / E<sup>(1)</sup><sub>i,σ</sub></th><th>σ = 1</th><th>σ = 2</th><th>σ = 5</th><th>σ = 10</th><th>σ = 20</th><th>σ = 50</th><th>σ = 80</th></tr></thead><tbody><tr><td>bag</td><td>0.86</td><td>0.48</td><td>0.28</td><td>0.59</td><td>0.45</td><td>1.10</td><td>1.05</td></tr><tr><td>building1</td><td>0.19</td><td>0.16</td><td>0.06</td><td>0.29</td><td>0.52</td><td>1.16</td><td>1.66</td></tr><tr><td>computer</td><td>0.20</td><td>0.10</td><td>0.19</td><td>0.21</td><td>0.48</td><td>2.26</td><td>1.82</td></tr><tr><td>dice</td><td>0.12</td><td>0.05</td><td>0.11</td><td>0.18</td><td>0.43</td><td>0.94</td><td>2.17</td></tr><tr><td>flowers2</td><td>0.19</td><td>0.08</td><td>0.13</td><td>0.30</td><td>0.50</td><td>1.60</td><td>0.89</td></tr><tr><td>hose</td><td>0.86</td><td>0.60</td><td>0.39</td><td>0.42</td><td>0.58</td><td>1.68</td><td>1.18</td></tr><tr><td>lawn</td><td>1.46</td><td>1.25</td><td>0.68</td><td>0.47</td><td>0.51</td><td>1.40</td><td>1.92</td></tr><tr><td>leaves</td><td>1.47</td><td>1.09</td><td>0.65</td><td>0.56</td><td>0.44</td><td>1.43</td><td>2.17</td></tr><tr><td>stairs</td><td>0.59</td><td>0.32</td><td>0.34</td><td>0.28</td><td>0.49</td><td>0.80</td><td>1.30</td></tr><tr><td>traffic</td><td>0.13</td><td>0.09</td><td>0.21</td><td>0.22</td><td>0.64</td><td>1.44</td><td>2.21</td></tr><tr><td>Flat image</td><td>0.02</td><td>0.03</td><td>0.05</td><td>0.17</td><td>0.37</td><td>1.34</td><td>1.23</td></tr><tr><td>E<sup>(2)</sup><sub>σ</sub></td><td><strong>0.56</strong></td><td><strong>0.39</strong></td><td><strong>0.28</strong></td><td><strong>0.34</strong></td><td><strong>0.49</strong></td><td><strong>1.38</strong></td><td><strong>1.60</strong></td></tr></tbody></table>

<table><caption>TABLE 2. This table shows, for the Percentile method, the $E_{i,σ}^{(1)}$ RMSE after adding simulated noise to the set of noise-free images (Figure 1) with several values of standard deviation $σ$. The last row is the $E_{σ}^{(2)}$ RMSE using the estimated $\hat{\sigma}_{i,b}$ of all the images. The parameters used in the tests are percentile $p = 0.005$, block of size 15 × 15, DCT filter with support 7 × 7 and seven bins.</caption><thead><tr><th>Image / E<sup>(1)</sup><sub>i,σ</sub></th><th>σ = 1</th><th>σ = 2</th><th>σ = 5</th><th>σ = 10</th><th>σ = 20</th><th>σ = 50</th><th>σ = 80</th></tr></thead><tbody><tr><td>bag</td><td>0.74</td><td>0.62</td><td>0.47</td><td>0.37</td><td>0.81</td><td>0.73</td><td>2.26</td></tr><tr><td>building1</td><td>0.34</td><td>0.24</td><td>0.55</td><td>0.62</td><td>0.82</td><td>1.20</td><td>1.58</td></tr><tr><td>computer</td><td>0.35</td><td>0.36</td><td>0.55</td><td>0.64</td><td>0.86</td><td>1.25</td><td>3.10</td></tr><tr><td>dice</td><td>0.12</td><td>0.12</td><td>0.17</td><td>0.24</td><td>0.50</td><td>1.20</td><td>1.68</td></tr><tr><td>flowers2</td><td>0.15</td><td>0.13</td><td>0.15</td><td>0.27</td><td>0.81</td><td>1.54</td><td>2.82</td></tr><tr><td>hose</td><td>0.87</td><td>0.62</td><td>0.49</td><td>0.41</td><td>0.49</td><td>1.37</td><td>1.55</td></tr><tr><td>lawn</td><td>0.99</td><td>1.20</td><td>0.86</td><td>0.64</td><td>0.62</td><td>1.64</td><td>1.90</td></tr><tr><td>leaves</td><td>1.43</td><td>1.11</td><td>0.95</td><td>0.74</td><td>0.70</td><td>0.96</td><td>1.67</td></tr><tr><td>stairs</td><td>0.94</td><td>0.89</td><td>0.66</td><td>0.56</td><td>0.64</td><td>0.89</td><td>1.34</td></tr><tr><td>traffic</td><td>0.45</td><td>0.42</td><td>0.56</td><td>0.58</td><td>0.91</td><td>1.36</td><td>2.23</td></tr><tr><td>Flat image</td><td>0.02</td><td>0.04</td><td>0.15</td><td>0.10</td><td>0.20</td><td>1.25</td><td>2.63</td></tr><tr><td>E<sup>(2)</sup><sub>σ</sub></td><td><strong>0.58</strong></td><td><strong>0.52</strong></td><td><strong>0.51</strong></td><td><strong>0.47</strong></td><td><strong>0.67</strong></td><td><strong>1.22</strong></td><td><strong>2.07</strong></td></tr></tbody></table>

The last row is the RMSE obtained for a given $σ$ and all the images. It is denoted by $E_{σ}^{(2)}$ and defined as

$$ (57) \qquad E_{\sigma}^{(2)} := \sqrt{\frac{1}{|B||I|} \sum_{i=1}^{|I|} \sum_{b=1}^{|B|} |\hat{\sigma}_{i,b} - \sigma|^2} = \sqrt{\frac{1}{|I|} \sum_{i=1}^{|I|} (E_{i,\sigma}^{(1)})^2}. $$
<table><caption>TABLE 3. This table shows, for the PCA method, the $E_{i,σ}^{(1)}$ RMSE after adding simulated noise to the set of noise-free images (Figure 1) with several values of standard deviation $σ$. The last row is the $E_{σ}^{(2)}$ RMSE using the estimated $\hat{\sigma}_{i,b}$ of all the images. Note that the PCA method is able to give an accurate estimation of the noise even when $σ = 1$ when the image has large flat zones (building1, computer, dice, flowers2, traffic), but fails to give a good estimation for very textured images (bag, hose, leaves, lawn). For very high levels of noise ($σ = 50, σ = 80$) the estimation is inaccurate. Surprisingly, the method behaves better for small levels of noise. Three bins are used.</caption><thead><tr><th>Image / E<sup>(1)</sup><sub>i,σ</sub></th><th>σ = 1</th><th>σ = 2</th><th>σ = 5</th><th>σ = 10</th><th>σ = 20</th><th>σ = 50</th><th>σ = 80</th></tr></thead><tbody><tr><td>bag</td><td>0.65</td><td>0.49</td><td>0.29</td><td>0.27</td><td>0.40</td><td>6.04</td><td>22.98</td></tr><tr><td>building1</td><td>0.18</td><td>0.10</td><td>0.21</td><td>0.46</td><td>0.65</td><td>2.17</td><td>2.38</td></tr><tr><td>computer</td><td>0.14</td><td>0.11</td><td>0.17</td><td>0.37</td><td>0.67</td><td>1.10</td><td>3.43</td></tr><tr><td>dice</td><td>0.13</td><td>0.06</td><td>0.09</td><td>0.30</td><td>0.59</td><td>1.81</td><td>3.77</td></tr><tr><td>flowers2</td><td>0.17</td><td>0.10</td><td>0.07</td><td>0.22</td><td>0.81</td><td>2.78</td><td>9.28</td></tr><tr><td>hose</td><td>0.76</td><td>0.66</td><td>0.35</td><td>0.28</td><td>0.25</td><td>9.61</td><td>26.78</td></tr><tr><td>lawn</td><td>0.87</td><td>0.62</td><td>0.46</td><td>0.57</td><td>0.32</td><td>13.71</td><td>29.49</td></tr><tr><td>leaves</td><td>1.16</td><td>0.71</td><td>0.35</td><td>0.31</td><td>0.40</td><td>3.67</td><td>17.83</td></tr><tr><td>stairs</td><td>0.55</td><td>0.55</td><td>0.50</td><td>0.53</td><td>1.01</td><td>17.15</td><td>33.67</td></tr><tr><td>traffic</td><td>0.14</td><td>0.23</td><td>0.34</td><td>0.59</td><td>0.59</td><td>0.58</td><td>3.36</td></tr><tr><td>Flat image</td><td>0.58</td><td>1.18</td><td>2.90</td><td>5.85</td><td>11.59</td><td>28.92</td><td>46.80</td></tr><tr><td>E<sup>(2)</sup><sub>σ</sub></td><td><strong>0.49</strong></td><td><strong>0.44</strong></td><td><strong>0.52</strong></td><td><strong>0.89</strong></td><td><strong>1.57</strong></td><td><strong>7.96</strong></td><td><strong>18.16</strong></td></tr></tbody></table>

For completeness, the results corresponding to the estimation using just a single bin are shown in Tables 4 (Ponomarenko), 5 (Percentile), and 6 (PCA). Assuming homoscedastic noise is not a realistic model, but it is needed that the algorithms give accurate results with homoscedastic noise before attempting their adaptation to signal-dependent noise (see Chapter 2).

Tables 7 (Ponomarenko), 8 (Percentile), and 9 (PCA) show the obtained $E_{σ}^{(2)}$ RMSE depending on the number of the iterations of the noise curve filter (see Section 1.2). Using five filtering iterations seems to be safe for any image with independence of the standard deviation of the noise, the kind of textures or the number and type of the edges the image may contain. More than five iterations is useless.

## 2. Evaluation comparing the noise curve of the raw image with the ground truth

In this test, the noise curve obtained by the algorithm for the raw images in Figure 2 (12 bits/channel, ISO 1600, t=1/30s) was compared to the “ground truth” noise curve of the camera for that ISO and exposure times. The ground truth was obtained by computing for each pixel the standard deviation of a large burst [24] of fixed snapshots of the same calibration pattern (Figure 3).
<table><caption>TABLE 4. This table shows, for the Ponomarenko method, the $E_{1,σ}^{(1)}$ RMSE after adding simulated noise to the set of noise-free images (Figure 1) with several values of standard deviation $σ$. The last row is the $E_{σ}^{(2)}$ RMSE using the estimated $\hat{\sigma}_{1,b}$ of all the images. It corresponds to the original signal-independent method, using a single bin and the iterations show in Algorithm 26 to fix percentile K.</caption><thead><tr><th>Image / E<sup>(1)</sup><sub>i,σ</sub></th><th>σ = 1</th><th>σ = 2</th><th>σ = 5</th><th>σ = 10</th><th>σ = 20</th><th>σ = 50</th><th>σ = 80</th></tr></thead><tbody><tr><td>bag</td><td>1.78</td><td>1.57</td><td>0.73</td><td>0.32</td><td>0.26</td><td>0.13</td><td>0.16</td></tr><tr><td>building1</td><td>0.16</td><td>0.09</td><td>0.02</td><td>0.10</td><td>0.07</td><td>0.04</td><td>0.30</td></tr><tr><td>computer</td><td>0.15</td><td>0.05</td><td>0.08</td><td>0.05</td><td>0.05</td><td>0.15</td><td>0.23</td></tr><tr><td>dice</td><td>0.11</td><td>0.06</td><td>0.01</td><td>0.03</td><td>0.07</td><td>0.51</td><td>0.29</td></tr><tr><td>flowers2</td><td>0.09</td><td>0.03</td><td>0.02</td><td>0.05</td><td>0.11</td><td>0.17</td><td>0.07</td></tr><tr><td>hose</td><td>0.75</td><td>0.36</td><td>0.24</td><td>0.17</td><td>0.06</td><td>0.44</td><td>0.22</td></tr><tr><td>lawn</td><td>1.86</td><td>0.42</td><td>0.75</td><td>0.27</td><td>0.37</td><td>0.09</td><td>0.00</td></tr><tr><td>leaves</td><td>3.04</td><td>1.30</td><td>0.55</td><td>0.67</td><td>0.44</td><td>0.41</td><td>0.05</td></tr><tr><td>stairs</td><td>1.06</td><td>0.85</td><td>0.45</td><td>0.23</td><td>0.14</td><td>0.41</td><td>0.27</td></tr><tr><td>traffic</td><td>0.10</td><td>0.07</td><td>0.08</td><td>0.13</td><td>0.11</td><td>0.47</td><td>0.27</td></tr><tr><td>Flat image</td><td>0.00</td><td>0.01</td><td>0.01</td><td>0.02</td><td>0.05</td><td>0.23</td><td>0.22</td></tr><tr><td>E<sup>(2)</sup><sub>σ</sub></td><td><strong>0.83</strong></td><td><strong>0.44</strong></td><td><strong>0.27</strong></td><td><strong>0.18</strong></td><td><strong>0.16</strong></td><td><strong>0.28</strong></td><td><strong>0.19</strong></td></tr></tbody></table>

<table><caption>TABLE 5. This table shows, for the Percentile method, the $E_{1,σ}^{(1)}$ RMSE after adding simulated noise to the set of noise-free images (Figure 1) with several values of standard deviation $σ$. The last row is the $E_{σ}^{(2)}$ RMSE using the estimated $\hat{\sigma}_{1,b}$ of all the images. The parameters used in the tests are percentile $p = 0.005$, block of size $15 \times 15$, DCT filter with support $7 \times 7$ and a single bin.</caption><thead><tr><th>Image / E<sup>(1)</sup><sub>i,σ</sub></th><th>σ = 1</th><th>σ = 2</th><th>σ = 5</th><th>σ = 10</th><th>σ = 20</th><th>σ = 50</th><th>σ = 80</th></tr></thead><tbody><tr><td>bag</td><td>0.71</td><td>0.60</td><td>0.43</td><td>0.32</td><td>0.65</td><td>0.29</td><td>1.22</td></tr><tr><td>building1</td><td>0.21</td><td>0.22</td><td>0.29</td><td>0.21</td><td>0.51</td><td>0.55</td><td>1.24</td></tr><tr><td>computer</td><td>0.29</td><td>0.25</td><td>0.38</td><td>0.32</td><td>0.08</td><td>0.99</td><td>0.86</td></tr><tr><td>dice</td><td>0.11</td><td>0.11</td><td>0.05</td><td>0.08</td><td>0.29</td><td>1.06</td><td>1.24</td></tr><tr><td>flowers2</td><td>0.13</td><td>0.10</td><td>0.12</td><td>0.03</td><td>0.00</td><td>0.54</td><td>2.03</td></tr><tr><td>hose</td><td>0.60</td><td>0.44</td><td>0.28</td><td>0.19</td><td>0.11</td><td>0.29</td><td>0.48</td></tr><tr><td>lawn</td><td>0.68</td><td>0.81</td><td>0.65</td><td>0.50</td><td>0.28</td><td>0.74</td><td>0.39</td></tr><tr><td>leaves</td><td>1.34</td><td>1.10</td><td>0.91</td><td>0.64</td><td>0.51</td><td>0.76</td><td>0.05</td></tr><tr><td>stairs</td><td>0.82</td><td>0.80</td><td>0.59</td><td>0.40</td><td>0.56</td><td>0.38</td><td>0.52</td></tr><tr><td>traffic</td><td>0.32</td><td>0.33</td><td>0.42</td><td>0.21</td><td>0.58</td><td>0.50</td><td>0.34</td></tr><tr><td>Flat image</td><td>0.02</td><td>0.03</td><td>0.05</td><td>0.05</td><td>0.16</td><td>0.61</td><td>1.93</td></tr><tr><td>E<sup>(2)</sup><sub>σ</sub></td><td><strong>0.47</strong></td><td><strong>0.44</strong></td><td><strong>0.38</strong></td><td><strong>0.27</strong></td><td><strong>0.34</strong></td><td><strong>0.61</strong></td><td><strong>0.94</strong></td></tr></tbody></table>
<table><caption>TABLE 6. This table shows, for the PCA method, the $E_{1,σ}^{(1)}$ RMSE after adding simulated noise to the set of noise-free images (Figure 1) with several values of standard deviation $σ$. The last row is the $E_{σ}^{(2)}$ RMSE using the estimated $\hat{\sigma}_{1,b}$ of all the images. A single bin is used.</caption><thead><tr><th>Image / E<sup>(1)</sup><sub>i,σ</sub></th><th>σ = 1</th><th>σ = 2</th><th>σ = 5</th><th>σ = 10</th><th>σ = 20</th><th>σ = 50</th><th>σ = 80</th></tr></thead><tbody><tr><td>bag</td><td>0.55</td><td>0.40</td><td>0.21</td><td>0.08</td><td>0.35</td><td>0.15</td><td>1.27</td></tr><tr><td>building1</td><td>0.17</td><td>0.08</td><td>0.08</td><td>0.12</td><td>0.18</td><td>1.48</td><td>1.83</td></tr><tr><td>computer</td><td>0.11</td><td>0.05</td><td>0.12</td><td>0.28</td><td>0.25</td><td>0.69</td><td>3.03</td></tr><tr><td>dice</td><td>0.11</td><td>0.07</td><td>0.07</td><td>0.26</td><td>0.52</td><td>1.71</td><td>3.76</td></tr><tr><td>flowers2</td><td>0.14</td><td>0.09</td><td>0.02</td><td>0.08</td><td>0.54</td><td>1.87</td><td>2.76</td></tr><tr><td>hose</td><td>0.32</td><td>0.25</td><td>0.08</td><td>0.20</td><td>0.18</td><td>0.91</td><td>2.83</td></tr><tr><td>lawn</td><td>0.33</td><td>0.12</td><td>0.15</td><td>0.32</td><td>0.38</td><td>1.97</td><td>1.65</td></tr><tr><td>leaves</td><td>1.02</td><td>0.62</td><td>0.06</td><td>0.18</td><td>0.18</td><td>0.29</td><td>2.34</td></tr><tr><td>stairs</td><td>0.38</td><td>0.24</td><td>0.24</td><td>0.44</td><td>0.29</td><td>1.25</td><td>2.61</td></tr><tr><td>traffic</td><td>0.25</td><td>0.28</td><td>0.17</td><td>0.33</td><td>0.29</td><td>0.03</td><td>2.99</td></tr><tr><td>Flat image</td><td>0.04</td><td>0.08</td><td>0.20</td><td>0.40</td><td>0.85</td><td>2.47</td><td>3.22</td></tr><tr><td>E<sup>(2)</sup><sub>σ</sub></td><td><strong>0.31</strong></td><td><strong>0.21</strong></td><td><strong>0.13</strong></td><td><strong>0.25</strong></td><td><strong>0.37</strong></td><td><strong>1.17</strong></td><td><strong>2.57</strong></td></tr></tbody></table>

<table><caption>TABLE 7. This table shows, for the Ponomarenko method, the obtained $E_{σ}^{(2)}$ RMSE values depending on the number of iterations of the noise curve filtering and the standard deviation of the noise (see Section 1.2). The percentile $p = 0.005$ and seven bins are used. Five iterations is the recommend value, since using more iterations does not improve the result significantly and it could soften too much the noise curves for certain images.</caption><thead><tr><th>Image / E<sup>(2)</sup><sub>σ</sub></th><th>σ = 1</th><th>σ = 2</th><th>σ = 5</th><th>σ = 10</th><th>σ = 20</th><th>σ = 50</th><th>σ = 80</th></tr></thead><tbody><tr><td>No filtering</td><td>0.56</td><td>0.39</td><td>0.28</td><td>0.34</td><td>0.49</td><td>1.38</td><td>1.60</td></tr><tr><td>1 iteration</td><td>0.55</td><td>0.38</td><td>0.26</td><td>0.27</td><td>0.41</td><td>1.23</td><td>1.35</td></tr><tr><td>2 iterations</td><td>0.54</td><td>0.37</td><td>0.24</td><td>0.24</td><td>0.38</td><td>1.16</td><td>1.22</td></tr><tr><td>3 iterations</td><td>0.53</td><td>0.37</td><td>0.23</td><td>0.23</td><td>0.36</td><td>1.12</td><td>1.13</td></tr><tr><td>4 iterations</td><td>0.53</td><td>0.36</td><td>0.23</td><td>0.22</td><td>0.36</td><td>1.10</td><td>1.11</td></tr><tr><td>5 iterations</td><td>0.52</td><td>0.35</td><td>0.22</td><td>0.21</td><td>0.35</td><td>1.09</td><td>1.10</td></tr><tr><td>6 iterations</td><td>0.51</td><td>0.35</td><td>0.21</td><td>0.21</td><td>0.35</td><td>1.08</td><td>1.10</td></tr><tr><td>7 iterations</td><td>0.51</td><td>0.34</td><td>0.21</td><td>0.21</td><td>0.35</td><td>1.08</td><td>1.09</td></tr></tbody></table>

To get the ground truth of the camera, we fixed the ISO sensitivity (in this case at ISO 1600) and used four exposure times, $t \in \{1/30s, 1/250s, 1/400s, 1/640s\}$. For each exposure time about two hundred pictures of the pattern were taken (see Figure 3). After cropping the area of the image that does not contain the calibration pattern, the final size of the raw image was 1352 × 1952. Since
<table><caption>TABLE 8. This table shows, for the Percentile method, the obtained $E_{\sigma}^{(2)}$ RMSE values depending on the number of iterations of the noise curve filtering and the standard deviation of the noise (see Section 1.2). The parameters used are percentile $p = 0.005$, block of size $15 \times 15$, the DCT filter with support $7 \times 7$ and seven bins. Five iterations is the recommended value, since using more iterations does not improve the result significantly and it could soften too much the noise curves for certain images.</caption><thead><tr><th>Image / E<sub>σ</sub><sup>(2)</sup></th><th>σ = 1</th><th>σ = 2</th><th>σ = 5</th><th>σ = 10</th><th>σ = 20</th><th>σ = 50</th><th>σ = 80</th></tr></thead><tbody><tr><td>No filtering</td><td>0.58</td><td>0.52</td><td>0.51</td><td>0.47</td><td>0.67</td><td>1.22</td><td>2.07</td></tr><tr><td>1 iteration</td><td>0.57</td><td>0.52</td><td>0.49</td><td>0.44</td><td>0.60</td><td>1.12</td><td>1.84</td></tr><tr><td>2 iterations</td><td>0.56</td><td>0.52</td><td>0.49</td><td>0.42</td><td>0.57</td><td>1.07</td><td>1.75</td></tr><tr><td>3 iterations</td><td>0.56</td><td>0.51</td><td>0.48</td><td>0.41</td><td>0.55</td><td>1.04</td><td>1.70</td></tr><tr><td>4 iterations</td><td>0.55</td><td>0.50</td><td>0.47</td><td>0.39</td><td>0.53</td><td>1.02</td><td>1.68</td></tr><tr><td>5 iterations</td><td>0.55</td><td>0.50</td><td>0.46</td><td>0.38</td><td>0.52</td><td>1.00</td><td>1.67</td></tr><tr><td>6 iterations</td><td>0.54</td><td>0.49</td><td>0.46</td><td>0.38</td><td>0.51</td><td>0.99</td><td>1.66</td></tr><tr><td>7 iterations</td><td>0.54</td><td>0.48</td><td>0.45</td><td>0.37</td><td>0.50</td><td>0.98</td><td>1.65</td></tr></tbody></table>

<table><caption>TABLE 9. This table shows, for the PCA method, the obtained $E_{\sigma}^{(2)}$ RMSE values depending on the number of iterations of the noise curve filtering and the standard deviation of the noise (see Section 1.2). Three bins are used. Five iterations is the recommended value, since using more iterations does not improve the result significantly and it could smooth too much the noise curves for certain images.</caption><thead><tr><th>Image / E<sub>σ</sub><sup>(2)</sup></th><th>σ = 1</th><th>σ = 2</th><th>σ = 5</th><th>σ = 10</th><th>σ = 20</th><th>σ = 50</th><th>σ = 80</th></tr></thead><tbody><tr><td>No filtering</td><td>0.49</td><td>0.44</td><td>0.52</td><td>0.89</td><td>1.57</td><td>7.96</td><td>18.16</td></tr><tr><td>1 iteration</td><td>0.49</td><td>0.42</td><td>0.48</td><td>0.82</td><td>1.46</td><td>7.22</td><td>16.90</td></tr><tr><td>2 iterations</td><td>0.49</td><td>0.41</td><td>0.46</td><td>0.78</td><td>1.39</td><td>6.66</td><td>15.93</td></tr><tr><td>3 iterations</td><td>0.49</td><td>0.41</td><td>0.45</td><td>0.76</td><td>1.35</td><td>6.23</td><td>15.17</td></tr><tr><td>4 iterations</td><td>0.48</td><td>0.40</td><td>0.45</td><td>0.75</td><td>1.34</td><td>6.25</td><td>15.19</td></tr><tr><td>5 iterations</td><td>0.48</td><td>0.40</td><td>0.45</td><td>0.75</td><td>1.34</td><td>6.26</td><td>15.20</td></tr><tr><td>6 iterations</td><td>0.48</td><td>0.40</td><td>0.45</td><td>0.74</td><td>1.33</td><td>6.27</td><td>15.21</td></tr><tr><td>7 iterations</td><td>0.47</td><td>0.40</td><td>0.45</td><td>0.74</td><td>1.33</td><td>6.28</td><td>15.23</td></tr></tbody></table>

each $2 \times 2$ block of the CFA<fnref n="1" /> contains one sample of the red channel, two samples of the green channel and one sample of the blue channel, one of the green channels can be discarded to get a single color pixel of each $2 \times 2$ block of the CFA, given an effective size of the color raw image of $676 \times 976$ pixels. The position of the camera was fixed when taking the snapshots of the image. Assuming constant lighting, the variance along several samples of image coming from different

<footnote n="1">¹Color Filter Array.</footnote>
<figure><img src="image_113.png" /><figcaption>FIGURE 2. Set of raw images used to test the noise estimation algorithm using 25 bins and without any noise curve filtering. The images are raw 12 bits/channel, taken with a Canon EOS 30D camera, ISO 1600 and exposure time t=1/30s. From left to right and up to bottom: images 1, 2, 3, 4, 5, 6, 7, and 8.</figcaption></figure>

<figure><img src="image_114.png" /><figcaption>FIGURE 3. One of the pictures of the calibration pattern mire used to build the ground truth noise curve of the camera.</figcaption></figure>

images for the same pixel position can only be explained by the noise. Therefore, it was possible to measure the mean of a block and the temporal standard deviation along all the snapshots to create an association mean→standard deviation, that is, a ground truth for camera noise curve, given the ISO and exposure times. Moreover, since the exposure time only affects the photon count and not the noise model, it was possible to overlap the noise curves for the four exposure times tested in a single ground truth noise curve depending only on the ISO parameter (see Figure 4).

Given an estimated noise curve $A$ of a test image, its control points are the pairs $(\hat{\mu}_{A,i,b}, \hat{\sigma}_{A,i,b}) \in A$ where $\hat{\mu}_{A,i,b}$ is the mean intensity for bin $b$ and image $i$ in $A$ and $\hat{\sigma}_{A,i,b}$ is the corresponding standard deviation values for bin $b$ and image $i$ in $A$. In the same way, given a ground truth noise curve $G$, its control points are the pairs $(\mu_{G,v}, \sigma_{G,v}) \in G$. Unfortunately the means of the noise
<figure><img src="image_115.png" /><figcaption>FIGURE 4. Ground truth of the Canon EOS 30D camera with ISO=1600.</figcaption></figure>

<figure><img src="image_116.png" /><figcaption>FIGURE 5. Checking a noise curve A (red) against the ground truth G (green), where *i* is the index of the image, *b* is the index of the bin, ($\hat{\mu}_{A,i,b}$, $\hat{\sigma}_{A,i,b}$) are the control points of the noise curve A, ($\mu_{G,v}$, $\sigma_{G,v}$) are the control points of G, and $\tilde{\sigma}_{G,i,i,b}$ is the standard deviation value projected from A into G.</figcaption></figure>

curve A and those in G do not necessarily coincide; that is, $\hat{\mu}_{A,i,b} \neq \mu_{G,v}$ for most $(b,v)$ pairs. To solve this problem, instead of using G, a new ground truth curve $\tilde{G}_i$ is used. This $\tilde{G}_i$ curve has the same means $\hat{\mu}_{A,i,b}$ as A (and therefore the same number of bins), and its standard deviation values are obtained by a simple proportionality rule. Therefore, the control points in the new curve $\tilde{G}_i$ are $(\hat{\mu}_{A,i,b}, \tilde{\sigma}_{G_i,i,b}) = (\hat{\mu}_{A,i,b}, \frac{\sigma_{G,v+1}-\sigma_{G,v}}{\mu_{G,v+1}-\mu_{G,v}}(\hat{\mu}_{A,i,b} - \mu_{G,v+1}) + \sigma_{G,v})$ where *v* is the index of the bin in the curve G such that $\mu_{G,v} \le \hat{\mu}_{A,i,b} < \mu_{G,v+1}$ (see Figure 5).
The error between the ground truth noise curve G and the test noise curve A for the image i and bin b is defined as

$$E_{G,A,i,b}^{(3)} := |\tilde{\sigma}_{G,i,i,b} - \hat{\sigma}_{A,i,b}| = \left| \frac{(\sigma_{G,v+1} - \sigma_{G,v})(\hat{\mu}_{A,i,b} - \mu_{G,v+1})}{\mu_{G,v+1} - \mu_{G,v}} + \sigma_{G,v} - \hat{\sigma}_{A,i,b} \right|$$

The values of $E_{G,A,i,b}^{(3)}$ for each test image in Figure 2 and each method are shown in Table 10. The Ponomarenko method gives systematically the smallest error. The PCA method is significantly sensitive to textures: the error is worse for those images that are highly textured (for example, Img. 3) and better when they contain large flat zones (for example, Img. 6 and Img. 7).

<table><caption>TABLE 10. Values of $E_{G,A,i,b}^{(3)}$ measuring the error between the noise curve A obtained for each test image i (Figure 2) and the ground truth curve G for the Canon EOS 30D with ISO 1600. The Ponomarenko method gives systematically the smallest error. The PCA method is significatly sensitive to textures: the error is worse for those images that are highly textured (for example, Img. 3) and better when they contain large flat zones (for example, Img. 6 and Img. 7).</caption><thead><tr><th>Method</th><th>Img. 1</th><th>Img. 2</th><th>Img. 3</th><th>Img. 4</th><th>Img. 5</th><th>Img. 6</th><th>Img. 7</th><th>Img. 8</th></tr></thead><tbody><tr><td>Ponomarenko</td><td><b>0.802</b></td><td><b>0.381</b></td><td><b>0.372</b></td><td><b>0.307</b></td><td><b>0.437</b></td><td><b>0.694</b></td><td><b>0.436</b></td><td><b>0.580</b></td></tr><tr><td>Percentile</td><td>0.910</td><td>0.680</td><td>0.836</td><td>0.850</td><td>0.871</td><td>0.867</td><td>0.779</td><td>0.669</td></tr><tr><td>PCA</td><td>0.984</td><td>1.107</td><td>3.577</td><td>1.640</td><td>2.103</td><td>1.010</td><td>0.788</td><td>1.381</td></tr></tbody></table>

To test the average behavior of the algorithm in all the bins of any test image, we define a mean error function $E_{G,A,b}^{(4)}$ as the mean of the $E_{G,A,i,b}^{(3)}$ values over the $|I|$ images for each of the $|B|$ bins, that is,

$$E_{G,A,b}^{(4)} := \frac{1}{|I|} \sum_{i=0}^{|I|-1} E_{G,A,i,b}^{(3)} \quad (58)$$

Figure 6 shows the $E_{G,A,b}^{(4)}$ error for all the 49 bins (25 bins in the PCA method) of the test images in Figure 2 for the first green channel, in the three methods.

### 3. Evaluation of the multiscale coherence of the result

Consider the down-scaling operator $S$ that tessellates the image into 2 × 2 pixels blocks, and replaces each block by a pixel having the mean of the four previous pixels as new value. If $\tilde{U}$ is a discrete pure Gaussian noise image with standard deviation $\sigma$, then $S(\tilde{U})$ has standard deviation $\frac{\sigma}{2}$. Indeed, if a block $W$ contains the pixels $\{u_1, u_2, u_3, u_4\}$ each one with variance $\sigma^2$, the variance of the mean of $W$ is $\text{Var}\tilde{W} = \text{Var}\frac{u_1+u_2+u_3+u_4}{4} = \frac{1}{16}[\text{Var}u_1 + \text{Var}u_2 + \text{Var}u_3 + \text{Var}u_4] = \frac{1}{16}[4\sigma^2] = \frac{\sigma^2}{4}$. Therefore, the standard deviation is $\text{Std}(\tilde{W}) = \frac{\sigma}{2}$; the noise has been divided by two. The objective of this test is to check if the noise estimation algorithm indeed divides the noise by two when the
<figure><img src="image_117.png" /><figcaption>FIGURE 6. Mean error $E_{G,A,b}^{(4)}$ for the 49 bins (25 in the case of the PCA method) of all the eight test images in Figure 2, for the Ponomarenko (a), Percentile (b) and PCA (c) methods.</figcaption></figure>

image is down-scaled several times.

$$ (59) \qquad \text{Set: } E_{A_0, A_k, i, b}^{(5)} = \left| \frac{\hat{\sigma}_{A_0, i, b}}{\hat{\sigma}_{A_k, i, b}} - 2^k \right| = \left| \frac{(\hat{\sigma}_{A_0, i, v+1} - \hat{\sigma}_{A_0, i, v})(\hat{\mu}_{A_k, i, b} - \hat{\mu}_{A_0, i, v+1})}{\hat{\sigma}_{A_k, i, b}(\hat{\mu}_{A_0, i, v+1} - \hat{\mu}_{A_0, i, v})} + \frac{\hat{\sigma}_{A_0, i, v}}{\hat{\sigma}_{A_k, i, b}} - 2^k \right|, $$

where

* $A_k$ is the noise curve corresponding to the input image $i$ after applying the down-scaling operator $k$ times. For example, if $k = 2$ then $A$ corresponds to the curve of the noise estimation of $SS(\tilde{U})$.

* $i$ is the image index, for the raw images in Fig. 2, $1 \le i \le |I|$, where $|I|$ is the number of images. $|I| = 8$ images were used.

* $b$ is the bin index, $1 \le b \le |B_k|$ where $|B_k|$ is the number of bins of the noise curve at scale $k$. For the test images $|B_0| = 49$, $|B_1| = 12$ and $|B_2| = 3$ bins are used.

* $v$ is the index of the bin in the curve $A_0$ such that $\hat{\mu}_{A_0, i, v} \le \hat{\mu}_{A_k, i, b} < \hat{\mu}_{A_0, i, v+1}$ (see Figure 7).

Remark: since the noise should be divided by two when operator $S$ is applied and an ideal noise estimator is used, the relation $\frac{\hat{\sigma}_{A_0, i, b}}{\hat{\sigma}_{A_k, i, b}}$ between the standard deviation estimations at scale 0 and scale $k$ (applying $k$ times $D$) should be equal to $2^k$. The error $E_{A_0, A_k, i, b}^{(5)}$ measures, for the tested noise estimator, the absolute deviation from the ideal value $2^k$ at each bin. To get a mean estimation of the error $E_{A_0, A_k, i, b}^{(5)}$ along all the test images and bins in Figure 2, we define another error function as

$$ (60) \qquad E_{A_0, A_k, b}^{(6)} := \frac{1}{|I| |B_k|} \sum_{i=1}^{|I|} \sum_{b=1}^{|B_k|} E_{A_0, A_k, i, b}^{(5)} $$

Table 11 shows the obtained mean down-scale error $E_{A_0, A_k, b}^{(6)}$ for the raw images in Figure 2 depending on the scale $k$, for the three methods. The measurements are done for one of the green channels of the raw images. The results of the PCA method are significantly worse than those given by the Ponomarenko et al. [15] and the Percentile [17] methods. The reason is that the
<figure><img src="image_118.png" /><figcaption>FIGURE 7. Checking a noise curve $A_k$ at scale $k$ (red) against the noise curve $A_0$ of the same image at scale 0 (green), where $i$ is the index of the image, $b$ is the index of the bin, ($\hat{\mu}_{A_k, i, b}$, $\hat{\sigma}_{A_k, i, b}$) are the control points of the noise curve of the sub-scaled image, ($\hat{\mu}_{A_0, i, v}$, $\hat{\sigma}_{A_0, i, v}$) are the control points of the noise curve of the image at scale 0, and $\tilde{\sigma}_{A_0, i, b}$ is the standard deviation value projected from $A_k$ into $A_0$.</figcaption></figure>

<table><caption>TABLE 11. Evaluation $E_{A_0, A_k, b}^{(6)}$ for the raw images in Figure 2 depending on the scale $k$. The measurements are done for one of the green channels of the raw images.</caption><thead><tr><th>Method</th><th>$k=1$</th><th>$k=2$</th><th>$k=3$</th></tr></thead><tbody><tr><td>Ponomarenko</td><td>0.195</td><td><strong>0.738</strong></td><td><strong>1.684</strong></td></tr><tr><td>Percentile</td><td><strong>0.158</strong></td><td>0.922</td><td>3.416</td></tr><tr><td>PCA</td><td>0.627</td><td>2.407</td><td>4.997</td></tr></tbody></table>

PCA method needs 112000 samples/bin while Ponomarenko et al. and the Percentile need 42000 samples/bin (see Section 1.1 of Chapter 9). When the image is down-scaled it becomes too small for the PCA method to get an accurate estimation, since the minimum number of samples/bin is not attained and only one bin is used for the third and fourth scales.

## 4. Online demo

An online demo is available for all three methods in IPOL. The users can upload any image to measure its noise. The demo also offers several types of pre-uploaded images to test the algorithm:

* Raw images obtained by splitting the raw channels $R, G_1, G_2, B$ and leaving out the $G_2$ channel. Then, an RGB image is formed by using the $R, G_1, B$ channels. Since in the
raw image no gamma correction has been done yet, the values of the image are multiplied by 32 to increase their dynamics and screen visibility. The colors of these images are not quite adapted to human visualization, because no white balance has been applied to them.

* The JPEG versions of the same raw images, as they are encoded by the camera.

* Various JPEG images.

* High SNR raw images, down-scaled by eight with their color channels averaged, so that they are nearly noiseless. In the demo they are referred to as "no noise" images.

Once an image has been chosen the following parameters can be configured:

* **Treatment of groups (2 × 2) of equal pixels.** It allows to choose between ignoring the blocks that contain a group of four equal pixels in any channel (default), or using all the blocks (see Section 1.3 of Chapter 9) unconditionally.

* **Curve filter iterations.** It indicates the number of filtering iterations that are applied to filter the noise curve (see Section 1.2). Default: five iterations.

* **Number of bins.** It is the number of bins in the noise curve (see Section 1.1 of Chapter 9). If it is set to *automatic selection*, each bin will contain approximately 112000 samples/bin.

* **A and B noise parameters.** Add a simulated noise with variance $A + B\tilde{U}$ is added to the input image. If $A = B = 0$ no noise will be added. If $B = 0$ white Gaussian noise with variance $A$ will be added. Default: $A = B = 0$.

For the Ponomarenko method, the following parameters can be configured:

**Percentile.:** The possible values are 0.01%, 0.1%, 0.5% (default), 5%, 10%, 50%. It can be also configured to use the iterations of the original method in order to find the percentile $K$ (see Algorithm 26).

**Block size.:** The possible choices are 3 × 3, 5 × 5, 7 × 7, 8 × 8 (default), 11 × 11, 15 × 15 and 21 × 21.

**Mean of blocks computation.:** permits to choose how the intensity associated to each bin is calculated. The possible choices are the average of the mean value of the pixels of the blocks that belong to the bin, or the median of the pixels of the blocks that belong to the bin (default).

**Curve filter iterations.:** It indicates the number of filtering iterations that are applied to filter the noise curve (see Section 1.2). Default: five iterations.

**Treatment of groups (2 × 2) of equal pixels.:** It allows to choose between ignoring the blocks that contain a group of four equal pixels in any channel (default), or using all the blocks unconditionally (see Section 1.3 of Chapter 9).

**Number of bins.:** It is the number of bins in the noise curve (see Section 1.1 of Chapter 9). The number of bins that are used depends on the size of the image when "automatic selection" is chosen. First, a nearest compatible size of the image is considered. For images whose
size is $S_0$, $S_1$ or $S_2$, the number of bins is given by dividing the total number of pixels of the image by 42000. Default: automatic selection.

**A and B noise parameters.:** Add a simulated noise with variance $A+B\tilde{U}$ is added to the input image. If $A = B = 0$ no noise will be added. If $B = 0$ white Gaussian noise with variance $A$ will be added. Default: $A = B = 0$.

For the Percentile method, the following parameters can be configured:

* Percentile. The possible values are 0.01%, 0.1%, 0.5% (default), 5%, 10% and 50%.

* Pre-filter operator. Operator **R** whose stencil **F** is used to convolve the image with. The possible operators are: Identity (no filter), Directional derivative, Laplace, Laplace (2 iterations), Laplace (3 iterations), Laplace (4 iterations), DCT with support $7 \times 7$ (default), DCT with support $5 \times 5$, DCT with support $3 \times 3$ and the filter of the article Fast Noise Variance Estimation [4].

* Block size. Possible choices: $3 \times 3$, $7 \times 7$, $8 \times 8$, $15 \times 21$ and $21 \times 21$ (default).

* Curve filter iterations.

* Treatment of groups ($2 \times 2$) of equal pixels.

* Number of bins.

* A and B noise parameters.

For the PCA method, the following parameters can be configured:

* Treatment of groups ($2 \times 2$) of equal pixels.

* Curve filter iterations.

* Number of bins.

* A and B noise parameters.

**4.1. Subtraction of the quantization noise.** In the online demo, all the images are encoded using 8 bits/pixel/channel. This adds a quantization error over the noise being estimated that must be subtracted. Indeed, the variance of a uniform random variable is $\sigma_q^2 = \int_{-\frac{1}{2}}^{\frac{1}{2}} (x-\bar{x})^2 dx = \int_{-\frac{1}{2}}^{\frac{1}{2}} x^2 dx = \left[ \frac{x^3}{3} \right]_{-\frac{1}{2}}^{\frac{1}{2}} = \frac{1}{12}$. This is the variance of the quantization error that must be subtracted at each scale. The standard deviation of the noise is computed at each bin as the square root of the noise variance computed directly by the algorithm minus the variance of the variance of the (independent) quantization error. At each scale $k$ the variance is divided by $4^k$ and thus the corrected standard deviation of the noise given by the demo is

$$ \tilde{\sigma}_k = \sqrt{\hat{\sigma}_k^2 - \frac{\sigma_q^2}{4^k}} = \sqrt{\hat{\sigma}_k^2 - \frac{1}{4^k 12}} $$

## 5. Complexity analysis of the algorithms

**5.1. Common subalgorithms.** This chapter discusses the algorithmic complexity of the sub-algorithms which are common to any method adapted to signal-dependent noise with the techniques explained in Chapter 9.
Algorithm 21 first executes the *argsort* (implemented with the Quicksort algorithm) operation with complexity $O(N \log N)$. The loop that iterates idx = 0...N just copies data in linear time $O(N)$. Therefore, Algorithm 21 is executed with complexity $O(N \log N)$. Algorithm 22, simply checks the condition at the *if* statement and computes simple arithmetic operations. Therefore, Algorithm 22 is executed in constant time $O(1)$. Algorithm 24 loops over the number of bins of the noise curve and inside the loop Algorithm 22 is called. Since Algorithm 22 is executed in constant time $O(1)$, Algorithm 24 has a linear complexity $O(B)$, being B the number of bins. Algorithm 25 loops over all possible pixels in the image (with the exception of the boundary of the image). The loop iterates through all the channels of the image and looks for groups of four connected pixels. Therefore, the inner loop is executed in linear time with the number of channels, $O(\text{num\_channels})$. It can actually be considered executed in constant time $O(1)$ once the number of channels have been fixed. The complexity of Algorithm 25 is given by its main loop, that is executed in linear time with complexity $O(N)$.

Algorithm 25 loops over all possible pixels in the image (with the exception of the boundary of the image). The loop iterates through all the channels of the image and looks for groups of four connected pixels. Therefore, the inner loop is executed in linear time with the number of channels, $O(\text{num\_channels})$. Since the number of channels is fixed and smaller than $M$, the complexity of Algorithm 25 is given by its main loop, that is executed in linear time with complexity $O(M)$.

**5.2. Analysis specific to Ponomarenko et al.** This complexity analysis is specific to the Ponomarenko et al. method.

The noise estimation procedure (Algorithm 26) first computes the DCT of all the $w \times w$ blocks in the image. Since the DCT is computed using the FFT algorithm, that has a complexity $O(M \log M)$. The loop that iterates seven times to get the optimal $K$ executes the *argmin* operation that involves the Quick-sort algorithm and therefore it can be done with $O(n \log n)$ with $n = |V_m^L|$.

The computation of $V_m^L = \frac{1}{\theta} \sum_{i=0}^{w-1} \sum_{j=0}^{w-1} [\mathbf{D}_m(i, j)]^2 \delta(i, j)$ and $V^H(i, j) = \frac{1}{K} \sum_{k=0}^{K-1} [\mathbf{D}_{(k)}(i, j)]^2$ have linear complexity $O(M)$ and $O(K)$, namely. Therefore, Algorithm 26 has a linear complexity $O(M)$ since $M > K$. Algorithm 27 performs the same operations with the exception of substituting the seven iterations to get $K$ with the product $k = pM$. Therefore, Algorithm 27 has also linear complexity $O(M)$.

**5.3. Analysis specific to Percentile.** This complexity analysis is specific to the Percentile method.

Algorithm 28 creates a matrix of $(N_x - 2b \times N_y - 2b)$ elements and then fills matrix $\mathbf{V}$ with $N_y - b - 1 \times N_x - b - 1$ values. Therefore, the complexity is linear, $O(N)$ with $N \sim N_x \times N_y$. The noise estimation procedure (Algorithm 29) computes the convolution $\tilde{\mathbf{U}} * \mathbf{F}$ in the Fourier domain using the FFT algorithm and then crops the result using Algorithm 28, that has linear complexity. The complexity of computing the convolution is $O(N \log N)$, where $N \sim (N_x - s + 1) \times (N_y - s + 1)$. Therefore, the complexity of cropping the convolution is $O(N \log N)$. Then, Algorithm 28 iterates $M$ times through a loop that reads the blocks $\mathbf{W}_m^f$ and computes its variance. Since $s$ is fixed,
the complexity of the loop is linear with the number of iterations, $O(M)$, where $M = (N_x - w - s + 1) \times (N_y - w - s + 1)$ is the number of overlapping blocks. After the execution of this loop, the SORTED function (implemented with the Quicksort algorithm) is called. Thus, the sorting operation has complexity $O(M \log M)$. Therefore, Algorithm 29 has complexity $O(M \log M)$.

### 5.4. Analysis specific to PCA. This complexity analysis is specific to the PCA method.

In Function `GetNextEstimate` (in Algorithm 30) the most computationally expensive part is the `ApplyPCA` function that is iterated several times until convergence is reached. The `ApplyPCA` function implies two steps:

* Computing the sample covariance matrix (45). If computed directly, it comprises about $|B(p)|M^2$ operations.

* Computing of the eigenvalues of the sample covariance matrix. Without any optimization, the number of operations is proportional to $M^3$. Again, since $M$ is fixed, the cost of this step can be considered constant, $O(1)$.

With the optimization explained in Section 4 of Chapter 12, the sample covariance matrix can be expressed as

$$ \frac{1}{|B(p)| - 1} \left( C_{B(p)} - \frac{1}{|B(p)|} \mathbf{c}_{B(p)} \mathbf{c}_{B(p)}^{T} \right) $$

where $C_X := \sum_{\mathbf{y}_i \in X} \mathbf{y}_i \mathbf{y}_i^T$ and $\mathbf{c}_X := \sum_{\mathbf{y}_i \in X} \mathbf{y}_i$.

Since the optimized algorithm first pre-computes the matrix $C_{B(1-j\Delta p)}$ and the vector $\mathbf{c}_{B(1-j\Delta p)}$ for $j = n-1, \dots, 0$, if the worst case is considered ($j = 0$), then $B(1-j\Delta p) = B(1)$, that is, all the $N$ blocks in the images would be used. The matrices and vectors are computed at most $\frac{1}{\Delta p} = \frac{1}{0.05} = 20$ times. To get the subsets of $BB(1-j\Delta p)$ according to $j$ it suffices to call just one the `argsort` function, implemented with the Quicksort algorithm with complexity $O(N \log N)$. Assuming that the matrices $C_{B(p)}$ and vectors $\mathbf{c}_{B(p)}$ and computed 20 times with the worst case $p = 1$, we have that $C_{B(p)} = C_{B(1)} = \sum_{\mathbf{y}_i \in B(1)} \mathbf{y}_i \mathbf{y}_i^T$ is executed with complexity $O(|B(1)|) = O(N)$. The same for $\mathbf{c}_{B(1)} := \sum_{\mathbf{y}_i \in B(1)} \mathbf{y}_i$, also with complexity $O(N)$. In summary, the computation of the sample covariance matrix consists on the execution of the `argsort` operation with complexity $N \log N$ and then looping at most 20 times through a loop that executes operations with complexity $O(N)$. Therefore, the complexity of Function `GetNextEstimate` (in Algorithm 30) is $N \log N$, being $N$ the number of overlapping $M \times M$ blocks in the image.

## 6. Conclusion

The Ponomarenko et al. method gave the best results according to its RMSE, when adapted to the signal-dependent noise. The Percentile method achieves a similar accuracy, but is slightly worse.

Although the PCA method gave the best results when the noise is a priori known to be homoscedastic, when it is signal-dependent the method loses accuracy. It is easily explained if we take into account that to get a reliable estimation of the noise variance it needs at least 112000
samples/bin, whereas the Ponomarenko et al. and the Percentile methods only require 42000 samples/bin.

The strategy follow by the Ponomarenko et al. method (measuring the variance of the noise at high-frequency coefficients in a low quantile of the set of blocks whose variance measured at the low-frequencies is minimal) gives the best results and can be considered the state of the art in signal-dependent noise estimation.
<figure><img src="image_119.png" /></figure>
Final conclusion

In the first part of this thesis, Noise Estimation, we have discussed several strategies to estimate
the noise, from the simplest techniques, that simply consist on simulating white Gaussian noise
and adding it to a noise-free image, signal-dependent noise estimators, to the most complete noise
estimators, that permit to estimate a complex model on which the noise depends on both the
intensity and the frequency.

For the homoscedastic white Gaussian noise estimators, Chapter 1 introduces the importance
of the noise estimation for denoising and reviews the most representative classic methods. However,
the simply homoscedastic noise model is only useful for raw images, where the noise is known to
depend on the intensity, as explained in Chapter 2. This chapter also proposes a way to obtain a
ground-truth noise curve for a particular camera configured with some ISO speed.

Chapter 3 showed the noise curves of the noise throughout the complete camera processing
pipeline. This chapter made it clear that the signal-dependent noise model is not enough to measure
the noise in most of the pictures, because after demosaicing the noise becomes highly correlated
and signal-dependent noise estimators give underestimations in these conditions.

Therefore, a new noise estimation algorithm was developed and presented in Chapter 8. This
new multiscale method is able to accurately measure the variance of highly correlated noise, even in
the case of JPEG-compressed images, which have undergone demosaicing, white balance, gamma
correction, and lossy JPEG-encoding. We introduced a new method to compute the difference
between to noise blocks, the *sparse distance*, which is robust to noise and compares only the
geometry of the blocks, without the interference of the noise. The comparison between the spatial
and the temporal noise curves demonstrates that indeed the method is valid. Other tests also
proved the method valid: a comparison of the noise estimates with the values expected after
filtering the image with a low-pass filter and the observation of denoising results.

In the second part of this thesis, Patch Denoising, we propose in Chapter 8 a new denoising
algorithm. It uses the noise estimation of the method presented in Chapter 4 to denoise each patch
of the noisy image. With this fully automatic tool, users input a noisy image with unknown noise
model, the noise is characterized, the noise is removed at each scale, and finally a denoised image
is returned to the users without any other interaction. We called this tool the "Noise Clinic",
as a metaphor of a real clinic, where the patients enter with some unknown problem, which is
diagnosed, and sometimes cured.

Definitely, this thesis studied in detail the problem of how to estimate noise to afterwards
denoise an image, using only the information of the noisy image (blind noise estimation and blind
denoising) and assuming only SFD noise. This minimal assumption of SFD permits to denoise
most of the images obtained with digital cameras, even when they are JPEG-compressed.

The third part of the thesis, Reproducible Research Contributions, presents three of the meth-
ods that were adapted to measure signal-dependent noise, along with some proposed general tech-
niques that permit to reliably adapt almost any block-based noise estimator to signal-dependent
noise. The articles presented in this third part were published in the Image Processing On Line
(IPOL) journal.

Some problems are still open and invite to further research, as the case of images with position-
dependent noise, adapting the size of the scanning window to the characteristics of the image, or
the search of optimal functions to compute the similarity between two noisy patches.
Bibliography

[1] M. Lebrun, M. Colom, A. Buades, and J.M. Morel, "Secrets of image denoising cuisine", *Acta Numerica*, vol. 21, pp. 475-576, 2012, doi: 10.1017/S0962492912000062.

[2] S.I. Olsen, "Estimation of noise in images: An evaluation", *CVGIP: Graphical Models and Image Processing*, vol. 55, no. 4, pp. 319-323, 1993.

[3] S. Pyatykh, Hesser J., and Zheng L., "Image noise level estimation by principal component analysis", *IEEE Transactions on Image Processing*, 2012, doi: 10.1109/TIP.2012.2221728.

[4] J. Immerkaer, "Fast noise variance estimation", *Computer Vision and Image Understanding*, vol. 64, no. 2, pp. 300-302, 1996, doi: 10.1006/cviu.1996.0060.

[5] K. Rank, M. Lendl, and R. Unbehauen, "Estimation of image noise variance", in *Vision, Image and Signal Processing*. IET, 1999, doi: 10.1049/ip-vis:19990238, vol. 146, pp. 80-84.

[6] P. Meer, J.M. Jolion, and A. Rosenfeld, "A fast parallel algorithm for blind estimation of noise variance", *IEEE Transactions on Pattern Analysis and Machine Intelligence*, vol. 12, no. 2, pp. 216-223, 1990, doi: 10.1109/34.44408.

[7] A. Foi, M. Trimeche, V. Katkovnik, and K. Egiazarian, "Practical poissonian-gaussian noise modeling and fitting for single-image raw-data", *IEEE Transactions on Image Processing*, vol. 17, no. 10, pp. 1737-1754, 2008, doi: 10.1109/TIP.2008.2001399.

[8] R. A. Boie and I. J. Cox, "An analysis of camera noise", *IEEE Transactions on Pattern Analysis and Machine Intelligence*, vol. 14, no. 6, pp. 671-674, 1992.

[9] M. Colom, A. Buades, and J.M. Morel, "Nonparametric noise estimation method for raw images", *Journal of the Optical Society of America A*, vol. 31, 2014. doi: 10.1364/JOSAA.31.000863.

[10] M. Colom, M. Lebrun, A. Buades, and J.M. Morel, "Multiscale estimation of intensity and frequency dependent noise", *IEEE Transactions On Image Processing*, 2014, Submitted.

[11] M. Raphan and E.P. Simoncelli, "An empirical bayesian interpretation and generalization of nl-means", Tech. Rep., Technical Report TR2010-934, Computer Science Technical Report, Courant Inst. of Mathematical Sciences, New York University, 2010.

[12] D. Van De Ville and M. Kocher, "Sure-based non-local means", *Signal Processing Letters, IEEE*, vol. 16, no. 11, pp. 973-976, 2009.

[13] C.A. Deledalle, V. Duval, and J. Salmon, "Non-local methods with shape-adaptive patches (NLM-SAP)", *Journal of Mathematical Imaging and Vision*, pp. 1-18, 2010.

[14] J. Portilla, "Blind non-white noise removal in images using Gaussian scale mixtures in the wavelet domain", *Benelux Signal Processing Symposium*, 2004.

[15] M. Colom and A. Buades, "Analysis and Extension of the Ponomarenko et al. Method, Estimating a Noise Curve from a Single Image", *Image Processing On Line*, vol. 3, pp. 173-197, 2013. doi: 10.5201/ipol.2013.45.

[16] N.N. Ponomarenko, V.V. Lukin, M.S. Zriakhov, A. Kaarna, and J.T. Astola, "An automatic approach to lossy compression of AVIRIS images", in *IEEE International Geoscience and Remote Sensing Symposium*, 2007, doi: 10.1109/IGARSS.2007.4422833, pp. 472-475.

[17] M. Colom and A. Buades, "Analysis and Extension of the Percentile Method, Estimating a Noise Curve from a Single Image", *Image Processing On Line*, vol. 2013, pp. 322-349, 2013. doi: 10.5201/ipol.2013.90.

[18] M. Colom and A. Buades, "Analysis and Extension of the PCA Method, Estimating a Noise Curve from a Single Image", *Image Processing On Line*, 2014, Submitted.
[19] M. Colom, G. Facciolo, M. Lebrun, N. Pierazzo, M. Rais, Y. Wang, and J.M. Morel, "A mathematical perspective of image denoising", *International Congress of Mathematicians*, 2014, To appear.

[20] M. Colom, M. Lebrun, A. Buades, and J.M. Morel, "A non-parametric approach for the estimation of intensity-frequency dependent noise", *IEEE International Conference on Image Processing*, 2014, Submitted.

[21] M. Lebrun, M. Colom, and J.M. Morel, "The Noise Clinic: a universal blind denoising algorithm", *IEEE International Conference on Image Processing*, 2014, Submitted.

[22] M. Lebrun, M. Colom, and J.M. Morel, "Multiscale image blind denoising", *IEEE Transactions On Image Processing*, 2014, Submitted.

[23] C. Chevalier, G. Roman, and J.N. Niepce, *Guide du photographe*, C. Chevalier, 1854.

[24] A. Buades, Y. Lou, J.M. Morel, and Z. Tang, "A note on multi-image denoising", in *Proceedings of the International Workshop on Local and Non-Local Approximation in Image Processing*. IEEE, 2009, doi: 10.1109/LNLA.2009.5278408, pp. 1–15.

[25] C.E. Shannon, "A mathematical theory of communication", *ACM SIGMOBILE Mobile Computing and Communications Review*, vol. 5, no. 1, pp. 3–55, 2001.

[26] S.R.J.L. HARRIS, "Image evaluation and restoration", JOSA, vol. 56, no. 5, pp. 569–570, 1966.

[27] W.H. Richardson, "Bayesian-based iterative method of image restoration", JOSA, vol. 62, no. 1, pp. 55–59, 1972.

[28] L. Yaroslavsky and M. Eden, "Fundamentals of Digital Optics", 2003.

[29] Y. Meyer, "Wavelets-algorithms and applications", *Wavelets-Algorithms and applications Society for Industrial and Applied Mathematics Translation*, 142 p., vol. 1, 1993.

[30] S. Geman and D. Geman, "Stochastic relaxation, gibbs distributions and the bayesian restoration of images", *IEEE Pat. Anal. Mach. Intell.*, vol. 6, pp. 721–741, 1984.

[31] A. Efros and T. Leung, "Texture synthesis by non parametric sampling", in *Proc. Int. Conf. Computer Vision*, 1999, vol. 2, pp. 1033–1038.

[32] A. Buades, B. Coll, J.M. Morel, et al., "A review of image denoising algorithms, with a new one", *Multiscale Modeling and Simulation*, vol. 4, no. 2, pp. 490–530, 2006.

[33] J. Portilla, V. Strela, M.J. Wainwright, and E.P. Simoncelli, "Image denoising using scale mixtures of gaussians in the wavelet domain", *Image Processing, IEEE Transactions on*, vol. 12, no. 11, pp. 1338–1351, 2003.

[34] A. Levin and B. Nadler, "Natural image denoising: Optimality and inherent bounds", in *Computer Vision and Pattern Recognition (CVPR), 2011 IEEE Conference on*. IEEE, 2011, pp. 2833–2840.

[35] P. Chatterjee and P. Milanfar, "Is denoising dead?", *Image Processing, IEEE Transactions on*, vol. 19, no. 4, pp. 895–911, 2010.

[36] K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian, "Image denoising by sparse 3D transform-domain collaborative filtering", *IEEE Transactions on image processing*, vol. 16, pp. 2007, 2007.

[37] J. Portilla, "Full blind denoising through noise covariance estimation using Gaussian scale mixtures in the wavelet domain", *Image Processing, 2004. ICIP'04. 2004 International Conference on*, vol. 2, pp. 1217–1220, 2004.

[38] B. Rajaei, "An analysis and improvement of the BLS-GSM denoising method", *Image Processing On Line*, vol. 2013, 2013, Preprint.

[39] A.B. Lee, K.S. Pedersen, and D. Mumford, "The nonlinear statistics of high-contrast patches in natural images", *International journal of Computer Vision*, vol. 54, no. 1, pp. 83–103, 2003.

[40] A. Foi, S. Alenius, and K. Katkovnik, V. and Egiazarian, "Noise measurement for raw-data of digital imaging sensors by automatic segmentation of non-uniform targets", *IEEE Sensors journal*, vol. 7, no. 10, pp. 1456–1461, 2007, doi: 10.1109/JSEN.2007.904864.

[41] Siméon-Denis Poisson, "Recherches sur la probabilité des jugements en criminelle et en matière civile, précédées des règles générales du calcul des probabilités", *Bachelier, Paris*, 1837.

[42] F. J. Anscombe, "The transformation of Poisson, binomial and negative-binomial data", *Biometrika*, vol. 35, no. 3, pp. 246–254, 1948. doi: 10.2307/2332343.
[43] M. Makitalo and A. Foi, "Optimal inversion of the Anscombe transformation in low-count Poisson image denoising", Image Processing, IEEE Transactions on, vol. 20, no. 1, pp. 99-109, 2011. doi: 10.1109/TIP.2010.2056693.

[44] Robert D. Nowak and Richard G. Baraniuk, "Wavelet-domain filtering for photon imaging systems", IEEE Transactions on Image Processing, vol. 8, no. 5, pp. 666-678, 1997. doi: 10.1109/83.760334.

[45] E. D. Kolaczyk, "Wavelet shrinkage estimation of certain Poisson intensity signals using corrected thresholds", Statist. Sin., vol. 9, pp. 119-135, 1999.

[46] Stamatios Lefkimmiatis, Petros Maragos, and George Papandreou, "Bayesian inference on multiscale models for poisson intensity estimation: Application to photo-limited image denoising", IEEE Transactions on Image Processing, vol. 18, no. 8, pp. 1724-1741, 2009.

[47] C.A. Deledalle, L. Denis, and F. Tupin, "Nl-insar: Nonlocal interferogram estimation", Geoscience and Remote Sensing, IEEE Transactions on, vol. 49, no. 4, pp. 1441-1452, 2011.

[48] C.A. Deledalle, F. Tupin, and L. Denis, "Poisson nl means: Unsupervised non local means for poisson noise", in Image Processing (ICIP), 2010 17th IEEE International Conference on. IEEE, 2010, pp. 801-804.

[49] C.A. Deledalle, F. Tupin, and L. Denis, "Polarimetric sar estimation based on non-local means", in Geoscience and Remote Sensing Symposium (IGARSS), 2010 IEEE International. IEEE, 2010, pp. 2515-2518.

[50] A. Foi, "Noise estimation and removal in mr imaging: The variance-stabilization approach", in Biomedical Imaging: From Nano to Macro, 2011 IEEE International Symposium on. IEEE, 2011, pp. 1809-1814.

[51] D. Zoran and Y. Weiss, "Scale invariance and noise in natural images", in Computer Vision, 2009 IEEE 12th International Conference on. IEEE, 2009, pp. 2209-2216.

[52] C. Liu, W. Freeman, R. Szeliski, and S. Kang, "Automatic estimation and removal of noise from a single image", IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 30, no. 2, pp. 299-314, February 2008.

[53] C. Liu, W.T. Freeman, R. Szeliski, and S.B. Kang, "Noise estimation from a single image", in Computer Vision and Pattern Recognition, 2006 IEEE Computer Society Conference on. Ieee, 2006, vol. 1, pp. 901-908.

[54] Aram Danielyan and Alessandro Foi, "Noise variance estimation in nonlocal transform domain", in Local and Non-Local Approximation in Image Processing, 2009. LNLA 2009. International Workshop on. IEEE, 2009, pp. 41-45.

[55] S.I. Olsen, "Estimation of noise in images: an evaluation", Computer Vision, Graphics, and Image Processing: Graphical Models and Image Processing, vol. 55, no. 4, pp. 319-323, July 1993. doi: 10.1006/cgip.1993.1022.

[56] J.S. Lee, "Refined filtering of image noise using local statistics", Computer graphics and image processing, vol. 15, no. 4, pp. 380-389, 1981.

[57] G.A. Mastin, "Adaptive filters for digital image noise smoothing: An evaluation", Computer Vision, Graphics, and Image Processing, vol. 31, no. 1, pp. 103-121, 1985, doi: 10.1016/S0734-189X(85)80078-5.

[58] R. Bracho and AC Sanderson, "Segmentation of images based on intensity gradient information", in Proc. IEEE Computer Soc. Conf. on Computer Vision and Pattern Recognition, 1985, pp. 19-23.

[59] H. Voorhees and T. Poggio, "Detecting textons and texture boundaries in natural images", in Proceedings of the First International Conference on Computer Vision, London. IEEE, Washington, DC, 1987, pp. 250-258.

[60] J. S Lee and K. Hoppel, "Noise modelling and estimation of remotely-sensed images", Proceedings of the International Geoscience and Remote Sensing Symposium, vol. 2, pp. 1005-1008, 1989. doi: 10.1109/IGARSS.1989.579061.

[61] D.L. Donoho and I.M. Johnstone, "Adapting to unknown smoothness via wavelet shrinkage", journal of the American Statistical Association, pp. 1200-1224, 1995, doi: 10.2307/2291512.

[62] P. Milanfar, "A tour of modern image filtering: New insights and methods, both practical and theoretical", Signal Processing Magazine, vol. 30, no. 1, pp. 106-128, 2013.

[63] Hossein Rabbani, Milan Sonka, and Michael D Abramoff, "Optical coherence tomography noise reduction using anisotropic local bivariate", International journal of Biomedical Imaging, vol. 3, no. 225, 2011.
[64] J. Schmitt, J.L. Starck, J.M. Casandjian, J. Fadili, and I. Grenier, "Multichannel poisson denoising and deconvolution on the sphere: Application to the fermi gamma ray space telescope", *Astronomy and Astrophysics*, vol. 546, pp. 10, 2012.

[65] Florian Luisier, Thierry Blu, and Michael Unser, "Image denoising in mixed poisson-gaussian noise", *Image Processing, IEEE Transactions on*, vol. 20, no. 3, pp. 696-708, 2011.

[66] F-X Dupé, Jalal M Fadili, and J-L Starck, "A proximal iteration for deconvolving poisson noisy images using sparse representations", *Image Processing, IEEE Transactions on*, vol. 18, no. 2, pp. 310-321, 2009.

[67] Hossein Rabbani, Reza Nezafat, and Saeed Gazor, "Wavelet-domain medical image denoising using bivariate laplacian mixture model", *Biomedical Engineering, IEEE Transactions on*, vol. 56, no. 12, pp. 2826-2837, 2009.

[68] Bo Zhang, Jalal M Fadili, and Jean-Luc Starck, "Wavelets, ridgelets, and curvelets for poisson noise removal", *Image Processing, IEEE Transactions on*, vol. 17, no. 7, pp. 1093-1108, 2008.

[69] Hossein Rabbani and Saeed Gazor, "Local probability distribution of natural signals in sparse domains", *International journal of Adaptive Control and Signal Processing*, 2013.

[70] Joseph Salmon Charles-Alban Deledalle and Arnak Dalalyan, "Image denoising with patch based pca: local versus global", in *Proceedings of the British Machine Vision Conference*. 2011. doi: 10.5244/C.25.25, pp. 25.1-25.10, BMVA Press.

[71] V.I.A. Katkovnik, V. Katkovnik, K. Egiazarian, and J. Astola, *Local approximation techniques in signal and image processing*, Society of Photo Optical, 2006.

[72] N.N. Ponomarenko, V.V. Lukin, S.K. Abramov, K.O. Egiazarian, and J.T. Astola, "Blind evaluation of adaptive noise variance in textured images by nonlinear processing of block DCT coefficients", in *Proceedings of the International Society for Optics and Photonics. Electronic Imaging, Image Processing: Algorithms and Systems II*, 2003, doi: 10.1117/12.477717, vol. 5014, pp. 178-189.

[73] D.L. Donoho and J.M. JOHNSTONE, "Ideal spatial adaptation by wavelet shrinkage", *Biometrika*, vol. 81, no. 3, pp. 425-455, 1994.

[74] D.L. Donoho, I.M. Johnstone, et al., "Adapting to unknown smoothness via wavelet shrinkage", *journal of the American Statistical Association*, vol. 90, no. 432, pp. 1200-1224, 1995.

[75] M.L. Uss, B. Vozel, V. V. Lukin, and K. Chehdi, "Image informative maps for component-wise estimating parameters of signal-dependent noise", *journal of Electronic Imaging*, vol. 22, no. 1, pp. 013019-013019, 2013.

[76] Mikhail Uss, Benoit Vozel, V Lukin, S Abramov, Igor Baryshev, and Kacem Chehdi, "Image informative maps for estimating noise standard deviation and texture parameters", *EURASIP journal on Advances in Signal Processing*, vol. 2011, no. 2011, pp. 806516, 2011.

[77] Poorvi L Vora, Joyce E Farrell, Jerome D Tietz, and David H Brainard, "Linear models for digital cameras", in *IS & T Annual Conference*. The Society for Imaging Science and Technology, 1997, pp. 377-382.

[78] M. Lebrun, A. Buades, and J.M. Morel, "A nonlocal bayesian image denoising algorithm", *SIAM journal on Imaging Sciences*, vol. 6, no. 3, pp. 1665-1688, 2013.

[79] G. K. Wallace, "The JPEG still picture compression standard", *Communications of the ACM*, vol. 34, no. 4, pp. 30-44, 1991.

[80] A. Buades, B. Coll, J. M. Morel, and C. Sbert, "Self-similarity driven demosaicking", *Image Processing On Line*, vol. 2011, 2011. doi: 10.5201/ipol.2011.bcms-ssdd.

[81] J. F. Hamilton Jr. and J. E. Adams Jr., "Adaptive color plan interpolation in single sensor color electronic camera", May 13 1997, US Patent 5,629,734.

[82] R. Neelamani, R. De Queiroz, Z. Fan, S. Dash, and R. G. Baraniuk, "JPEG compression history estimation for color images", *Image Processing, IEEE Transactions on*, vol. 15, no. 6, pp. 1365-1378, 2006.

[83] N. N. Ponomarenko, V. V. Lukin, K. O. Egiazarian, and J. T. Astola, "A method for blind estimation of spatially correlated noise characteristics", in *IS&T/SPIE Electronic Imaging*. International Society for Optics and Photonics, 2010, pp. 753208-753208.
[84] M. Lebrun, A. Buades, and J.M. Morel, "Implementation of the "non-local Bayes" (NL-Bayes) image denoising algorithm", Image Processing On Line, vol. 2013, pp. 1-42, 2013.

[85] F. J. Anscombe, "The transformation of Poisson, binomial and negative-binomial data", *Biometrika*, vol. 35, no. 3, pp. 246–254, 1948.

[86] A. Buades, B. Coll, and J.M. Morel, "A non local algorithm for image denoising", *IEEE Computer Vision and Pattern Recognition*, vol. 2, pp. 60–65, 2005.

[87] M. Aharon, M. Elad, and A. Bruckstein, "K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation", *Signal Processing, IEEE Transactions on*, vol. 54, no. 11, pp. 4311–4322, 2006.

[88] Pierre Brémaud, *Markov chains: Gibbs fields, Monte Carlo simulation, and queues*, vol. 31, springer verlag, 1999.

[89] P. Chatterjee and P. Milanfar, "Patch-based near-optimal image denoising.", *IEEE transactions on image processing: a publication of the IEEE Signal Processing Society*, 2011.

[90] Daniel Zoran and Yair Weiss, "From learning models of natural image patches to whole image restoration", in *Computer Vision (ICCV), 2011 IEEE International Conference on*. IEEE, 2011, pp. 479–486.

[91] L. Zhang, W. Dong, D. Zhang, and G. Shi, "Two-stage image denoising by principal component analysis with local pixel grouping", *Pattern Recognition*, vol. 43, no. 4, pp. 1531–1549, 2010.

[92] R.R. Coifman and D.L. Donoho, "Translation-invariant de-noising", *Lecture Notes In Statistics*, pp. 125–125, 1995.

[93] A. Nemirovski, "Topics in non-parametric statistics", *Lectures on probability theory and statistics (Saint-Flour, 1998)*, vol. 1738, pp. 85–277, 2000.

[94] O.G. Guleryuz, "Weighted averaging for denoising with overcomplete dictionaries", *Image Processing, IEEE Transactions on*, vol. 16, no. 12, pp. 3020–3034, 2007.

[95] Guillermo Sapiro Guoshen Yu, "DCT image denoising: a simple and effective image denoising algorithm", *Image Processing On Line*, 2011. doi: 10.5201/ipol.2011.ys-dct.

[96] C. Kervrann and J. Boulanger, "Local Adaptivity to Variable Smoothness for Exemplar-Based Image Regularization and Representation", *International Journal of Computer Vision*, vol. 79, no. 1, pp. 45–69, 2008.

[97] M. Raphan and E.P. Simoncelli, "Learning to be bayesian without supervision", *Advances in neural information processing systems*, vol. 19, pp. 1145, 2007.

[98] V. Duval, J.F. Aujol, and Y. Gousseau, "A bias-variance approach for the nonlocal means", *SIAM Journal on Imaging Sciences*, vol. 4, pp. 760, 2011.

[99] P. Milanfar, "A tour of modern image filtering", *Invited feature article to IEEE Signal Processing Magazine (preprint at http://users. soe. ucsc. edu/ milanfar/publications/)*, 2011.

[100] J.W. Tukey, "Exploratory data analysis. 1977", *Massachusetts: Addison-Wesley*, 1976.

[101] S. Osher, M. Burger, D. Goldfarb, J. Xu, and W. Yin, "Using geometry and iterated refinement for inverse problems (1): Total variation based image restoration", *Department of Mathematics, UCLA, LA, CA*, vol. 90095, pp. 04–13, 2004.

[102] Alexei A Efros and Thomas K Leung, "Texture synthesis by non-parametric sampling", in *Computer Vision, 1999. The Proceedings of the Seventh IEEE International Conference on*. IEEE, 1999, vol. 2, pp. 1033–1038.

[103] J.M. Morel, B. Coll, and A. Buades, "Image data processing method by reducing image noise, and camera integrating means for implementing said method", June 10 2010, US Patent App. 11/579,380.

[104] J.S. De Bonet, "Noise reduction through detection of signal redundancy", *Rethinking artificial intelligence*, 1997.

[105] T. Weissman, E. Ordentlich, G. Seroussi, S. Verdu, and MJ Weinberger, "Universal discrete denoising: Known channel", *IEEE Transactions on Information Theory*, vol. 51, no. 1, pp. 5–28, 2005.

[106] E. Ordentlich, G. Seroussi, S. Verdu, M. Weinberger, and T. Weissman, "A discrete universal denoiser and its application to binary images", in *International Conference on Image Processing*, 2003, vol. 1.

[107] S.P. Awate and R.T. Whitaker, "Unsupervised, information-theoretic, adaptive image filtering for image restoration", *IEEE Trans. PAMI*, vol. 28, no. 3, pp. 364–376, 2006.
[108] G. Yu and G. Sapiro, “DCT image denoising: a simple and effective image denoising algorithm”, Image Processing On Line, vol. 2011, 2011.

[109] P. Getreuer, “Rudin-osher-fatemi total variation denoising using split bregman”, Image Processing On Line, vol. 2012, 2012.

[110] M. Lebrun and A. Leclaire, “An implementation and detailed analysis of the k-svd image denoising algorithm”, Image Processing On Line, vol. 2012, 2012.

[111] A. Buades, B. Coll, and J.M. Morel, “Non-local means denoising”, Image Processing On Line, vol. 2011, 2011.

[112] M. Lebrun, “An analysis and implementation of the BM3D image denoising method”, Image Processing On Line, vol. 2012, 2012.

[113] T. Rabie, “Robust estimation approach for blind denoising”, IEEE Transactions on Image Processing, vol. 14, no. 11, pp. 1755–1765, November 2005.

[114] S. Gabarda and G. Cristóbal, “The generalized Rényi image entropy as a noise indicator”, Noise and Fluctuations in Photonics, Quantum Optics, and Communications, vol. 6603, 2007. doi: 10.1117/12.725086.

[115] J. Mairal, M. Elad, and G. Sapiro, “Sparse representation for color image restoration”, IEEE Transactions on Image Processing, vol. 17, no. 1, pp. 53–69, 2008.

[116] J. Mairal, G. Sapiro, and M. Elad, “Learning multiscale sparse representations for image and video restoration”, SIAM Multiscale Modeling and Simulation, vol. 7, no. 1, pp. 214–241, 2008.

[117] L.P. Yaroslavsky, K.O. Egiazarian, and J.T. Astola, “Transform domain image restoration methods: review, comparison, and interpretation”, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, May 2001, vol. 4304, pp. 155–169.

[118] L.P. Yaroslavsky, “Local adaptive image restoration and enhancement with the use of DFT and DCT in a running window”, in Proceedings of SPIE, 1996, vol. 2825, pp. 2–13.

[119] J.L. Starck, E.J. Candès, and D.L. Donoho, “The curvelet transform for image denoising”, IEEE Transactions on Image Processing, vol. 11, pp. 670–684, 2002.

[120] L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms”, Phys. D, vol. 60, pp. 259–268, 1992.

[121] R. R. Coifman and D. L. Donoho, Translation-invariant de-noising, vol. 103, Springer New York, 1995. doi: 10.1007/978-1-4612-2544-7_9.

[122] A. Buades, B. Coll, J.M. Morel, and C. Sbert, “Self-similarity Driven Demosaicking”, Image Processing On Line, vol. 1, 2011. doi: 10.5201/ipol.2011.bcms-ssdd.

[123] J.W. Cooley and J.W. Tukey, “An algorithm for the machine calculation of complex fourier series”, Mathematics of Computation, pp. 297–301, 1965, doi: 10.2307/2003354.

[124] M. Black and G. Sapiro, “Edges as outliers: Anisotropic smoothing using local image statistics”, in Scale-Space Theories in Computer Vision, Mads Nielsen, Peter Johansen, Ole Olsen, and Joachim Weickert, Eds., vol. 1682 of Lecture Notes in Computer Science, pp. 259–270. Springer Berlin / Heidelberg, 1999, doi: 10.1007/3-540-48236-9_23.

[125] S. A. Gershgorin, “Uber die abgrenzung der eigenwerte einer matrix”, Izvestija Akademii Nauk SSSR, pp. 749–754, 1931.