<figure><img src="image_36.png" /><figcaption>FIGURE 9. Comparison between the temporal series standard deviation and the noise curve obtained by the Ponomarenko et al. algorithm after gamma correction. The images were taken with ISO 1600 and t=1/250s. The solid lines correspond to the temporal series and the dashed to the single image estimation with the Ponomarenko et al. algorithm. The dynamic range is clearly increased specially in the dark areas. The noise increases significantly because of the power law function. When the intensity is over some threshold the noise level begins to decrease instead of increasing because the signal gets saturated. Left: Canon camera. Right: Nikon camera.</figcaption></figure>

color space but in the $Y'C_BC_R$, where $Y'$ is the luma component and $C_B$ and $C_R$ the blue-difference and red-difference chroma components, namely. Also, the number of bits is reduced to 8 bits/channel/pixel. Since human perception is much more sensitive to changes in luma than in chroma components, the $C_BC_R$ components are subsampled. Usually the 4 : 2 : 0 subsampling is performed, which means that both the horizontal and vertical resolutions are halved. This step implies loss of information. After the chroma subsampling, the image is tessellated into blocks of 8 × 8 pixels and the 2D DCT-II of each is computed. Human perception is not adapted to distinguish accurately differences in luma and chroma when the signal varies rapidly, that is, at its high frequencies. Therefore, it is possible to quantize the information at the high-frequencies without the notice of human observers. Also, it is well known that the variance of the high frequency coefficients in the 2D-DCT transformation of a block is mainly explained by the noise. Since JPEG encoding implies quantizing these high-frequency coefficients, most of the methods will not be able to detect the noise by using these high frequencies. For example, the Ponomarenko et al. method [16] which tries to estimate the noise using these high-frequency coefficients, gives an underestimation if the image is JPEG encoded.

For each coefficient $C_{i,j}$ of the DCT block of the image $(i, j \in [0, 7])$, JPEG applies the operation

$$J_{i,j} = \text{round} \left( \frac{C_{i,j}}{R_{i,j}} \right) R_{i,j},$$
<figure><img src="image_37.png" /><figcaption>FIGURE 10. Example of quantization matrix for JPEG encoding, as given in the standard [79].</figcaption></figure>

where $J_{i,j}$ is the resulting quantized coefficient and $R_{i,j}$ is the value of the quantization matrix for the coefficient at frequency $[i, j]$. Figure 10 shows the quantization matrix proposed by the JPEG standard [79].

Matrix $R$ contains low values for the low-frequencies and higher values for the high-frequencies, to quantize them more. It is designed such a way that it causes a great energy loss in the high frequencies (the energy at the high-frequencies comes mainly from edges and textures) whereas the medium and low frequencies (where most of the visual information is located, with the exception of edges) are mostly respected. The actual values of the matrix depend on the quality at which the JPEG encoder is configured, where the less quality, the higher the quantization. The quality of the final JPEG image is determined by a quality parameter that goes from $Q = 0$ (worst quality) to $Q = 100$ (best quality). Most of the cameras use quality $Q = 92$ as the standard quality by default.

As an example, let us consider an image made of samples of pure white Gaussian noise of mean zero and standard deviation $\sigma = 10$. Then, compute the 2D DCT-II of all non-overlapping $8 \times 8$ blocks in the image. If we denote by $D_k[i, j]$ the coefficient of frequency $[i, j] \in [0, 7]^2$ of block $k \in [0, M-1]$, then the empirically estimated variance of the noise at that frequency along all the blocks is $\text{Var}_k(D_k[i, j]) = \frac{1}{M-1} \sum_{k=0}^{M-1} (D_k[i, j])^2$, where $M$ is the number of non-overlapping blocks. Table 11 shows the averaged STD of all blocks according to the frequency. On the left, result from an image of pure white Gaussian noise of STD=10. On the right, the same image after JPEG encoding with quality $Q = 70$. The [0, 0] frequency (DC) is at the leftmost top corner. Note that after JPEG compression, the STD decreases as the frequency increases, because of the quantization matrix. The slightly increased energy overall the matrix is due to the blocking artifacts created by JPEG compression.

After quantization, the number of possible values for each DCT coefficient is reduced. In this situation the lossless Huffman encoder can achieve good compression ratios. This is the last step of the compression, that encodes the quantized coefficients found following a zig-zag scan of the DCT block. Since the energy of the signal in natural images decreases as the frequency increases,
<figure><img src="image_38.png" /><figcaption>FIGURE 11. Averaged variance along all the blocks at all frequencies. Up: result from an image of pure white Gaussian noise of STD=10. Down: the result using same image after JPEG encoding with quality $Q = 70$. The [0,0] frequency (DC) is at the leftmost top corner. Note that after JPEG compression, the STD decreases as the frequency increases, because of the quantization matrix. The slightly increased energy overall the matrix is due to the blocking artifacts created by JPEG compression.</figcaption></figure>

the encoder uses a special code word EOB (end-of-block) which indicates that next coefficients in the zig-zag scan are all zero. This improves even more the compression ratio.

To study the effect of JPEG encoding, the parameter $Q = 92$ was used to encode the images for the experiments of this chapter. Figure 12 shows the noise curves for both the Canon and Nikon cameras after JPEG encoding.

It can be observed that

* the dynamic range of the image has not changed after JPEG encoding.

* the noise is reduced after JPEG encoding, because of the quantization of the coefficients, particularly those corresponding to the high frequencies.

Figure 13 shows the autocorrelation function after JPEG encoding for both the Canon and Nikon cameras. Each pixel is strongly correlated with the four pixels at distance 1, and since JPEG removes most of the energy at the high-frequencies, it can also be observed some small correlations
<figure><img src="image_39.png" /><figcaption>FIGURE 12. Comparison between the temporal series standard deviation and the noise curve obtained by the Ponomarenko et al. algorithm after JPEG encoding. The images were taken with ISO 1600 and t=1/250s. The solid lines correspond to the temporal series and the dashed to the single image estimation with the Ponomarenko et al. algorithm. The dynamic range is not modified because of JPEG encoding. The curves are shown using a 12 bits intensity range to compare with the previous steps, although in reality the image is encoded with 8 bits instead of 12 bits/color channel. Left: Canon camera. Right: Nikon camera.</figcaption></figure>

with pixels with distances up to 5 pixels. Also, the Nikon camera shows a higher autocorrelation because the noise is higher than the noise of the Canon camera.

<figure><img src="image_40.png" /><figcaption>FIGURE 13. Autocorrelation function for a rectangle in the calibration pattern for the Canon (left) and Nikon (right) cameras with ISO 1600 and t=1/250s in the JPEG-encoded image. Each pixel is strongly correlated with the four pixels at distance 1, and since JPEG removes most of the energy at the high-frequencies, it can also be observed some small correlations with pixels with distances up to 5 pixels. Also, the Nikon camera shows a higher autocorrelation because the noise is higher than the noise of the Canon camera.</figcaption></figure>

## 1.3. Synthesis: effect of the complete image processing pipeline on the noise curve.

In this chapter we show the evolution of the noise curves throughout the complete processing
<figure><img src="image_41.png" /><figcaption>FIGURE 14. Complete pipeline for ISO 1250, t=1/30s, Canon: raw image, demosaicing, white balance, gamma (tone curve) correction and JPEG compression. In the first step (raw image), all four color channels show share the same noise curve. After demosaicing, each color channel have a different noise curve, since the Adams-Hamilton algorithm treats each channel in a different way. Finally, the gamma correction saturates the noise curve, which starts decrease from a certain intensity. The final JPEG noise curve exhibits the combination of all these effects along the processing chain.</figcaption></figure>

pipeline of the cameras, from the initial raw to the final JPEG image. We fixed an exposure time t=1/30s and two ISO speeds, ISO 1250 (Figure 14) and ISO 1600 Figure 15 for both the Canon and Nikon cameras. In the first step (raw image), all four color channels show share the same noise curve. After demosaicing, each color channel have a different noise curve, since the Adams-Hamilton algorithm treats each channel in a different way. Finally, the gamma correction saturates the noise curve, which starts decrease from a certain intensity. The final JPEG noise curve exhibits the combination of all these effects along the processing chain. Comparing the raw noise curves from both pipelines, it can be observed that the noise of the images with ISO 1600 is higher than the noise of the images with ISO 1250, as expected. Also, the dynamic range is larger with ISO 1600, since the sensitivity is higher. Apart of the different dynamic ranges, the curves from ISO 1250 and ISO 1600 show the same behavior along the processing chain. Figure 16 and Figure 17 show the curves corresponding to ISO 1250 and ISO 1600 for the Nikon camera.

The main difference between the Canon and Nikon cameras that can be observed is that when using the same exposure time and ISO, the dynamic range of the Nikon camera is larger in the Nikon compared to the Canon camera. In other terms, the Nikon camera is more sensitive to light than the Canon, under the same ISO configuration. This explains why in all measurements the
<figure><img src="image_42.png" /><figcaption>FIGURE 15. Complete pipeline for ISO 1600, t=1/30s, Canon: raw image, demosaicing, white balance, gamma (tone curve) correction and JPEG compression. In the first step (raw image), all four color channels show share the same noise curve. After demosaicing, each color channel have a different noise curve, since the Adams-Hamilton algorithm treats each channel in a different way. Finally, the gamma correction saturates the noise curve, which starts to decrease from a certain intensity. The final JPEG noise curve exhibits the combination of all these effects along the processing chain.</figcaption></figure>

Nikon camera always show more noise. Since the Nikon camera is more sensitive to light than the Canon, the gamma correction functions start to saturate the signal earlier. This effect can be seen looking at the curves after the gamma correction step. For example, the noise at the red channel begins to decrease at intensity 2000 approximately with both ISO 1250 and ISO 1600 in the Canon camera. In the case of the Nikon camera it begins at intensity 700, approximately.

## 2. Overlapping of noise curves with different exposure times

Historically, the ISO speed (standard ISO 5800:2001) was related to the sensitivity of the photographic film used to take a picture with analogic cameras. With modern digital cameras, the ISO level follows the ISO 12232:2006 norm, that sets the reference values for the calibration of digital still cameras. In general, the ISO speed is simulated by multiplying the readouts at the raw image by some fixed factors. As a consequence, the standard deviation of the noise varies linearly with the ISO factors. However, changing the exposure time using the same ISO does not alter the shape of the noise curve, but only increases the dynamic range with the exposure time. In Section 1.2, several noise curves corresponding to different exposure times were shown. Since the shape of the noise curve just depends on the ISO, it makes sense to combine the curves obtained at different
<figure><img src="image_43.png" /><figcaption>FIGURE 16. Complete pipeline for ISO 1250, t=1/30s, Nikon: raw image, demosaicing, white balance, gamma (tone curve) correction and JPEG compression. In the first step (raw image), all four color channels show share the same noise curve. After demosaicing, each color channel have a different noise curve, since the Adams-Hamilton algorithm treats each channel in a different way. Finally, the gamma correction saturates the noise curve, which starts de decrease from a certain intensity. The final JPEG noise curve exhibits the combination of all these effects along the processing chain.</figcaption></figure>

exposure times and same ISO to get a unique noise curve. If this procedure is applied to the raw image, the resulting curve can be used as a ground-truth for some ISO of a given camera. Figures 18 and 19 show the overlapping of the noise curves corresponding to 1/30s, 1/250s, 1/400s and 1/640s for the Canon and Nikon cameras at the raw, demosaicing, white balance, gamma correction and JPEG encoding steps.

## 3. Mean ground-truth curves

Section 2 showed that noise curves coming from different exposure times overlap at the raw image step. After the demosaicing step, the curves do not overlap anymore because the demosaicing algorithm computes different interpolations depending on the channel. Since the image at the raw step has not yet been processed by the demosaicing and further algorithms, it is only a function of the exposure time and ISO speed and can be used as a ground-truth curve. This is particularly useful to get an estimation of the accuracy of noise estimation algorithms. The noise curves corresponding to different exposure times should overlap. In practice, there is some small error due to the random nature of noise itself and the limited accuracy of the noise estimation algorithms. Figures 18 and 19 show that the noise curves overlap almost exactly in the raw image step, but
<figure><img src="image_44.png" /><figcaption>FIGURE 17. Complete pipeline for ISO 1600, t=1/30s, Nikon: raw image, demosaicing, white balance, gamma (tone curve) correction and JPEG compression. In the first step (raw image), all four color channels show share the same noise curve. After demosaicing, each color channel have a different noise curve, since the Adams-Hamilton algorithm treats each channel in a different way. Finally, the gamma correction saturates the noise curve, which starts de decrease from a certain intensity. The final JPEG noise curve exhibits the combination of all these effects along the processing chain.</figcaption></figure>

with some negligible error. To get a unique noise curve, the average of the values of the curves corresponding to different exposure times is computed (see Algorithm 13). Figures 20 and 21 show the computed ground-truth for the Canon camera for ISO 1250 and ISO 1600. Figures 22 and 23 show the same for the Nikon camera. On the right, a detailed view of the noise curve at darkest intensities.

The darkest zone of the noise curve has a particular interest, since the gamma correction increases the energy of the image specially in these dark zones, making the noise really noticeable afterwards. After the white balance step, the noise has been multiplied by a different factor that depends on the channel. If the noise is convolved by some kernel, this visible noise will appear to the observer as colored stains on the image. The DCT coefficient quantization step in JPEG is similar to convolution with a kernel that low-passes the image and therefore these kind of artifacts are still visible after JPEG compression.

4. Comparison of the autocorrelation functions at different scales

Consider the operator **S** that tessellates the image into sets of 2 × 2 pixels blocks, and replaces each block by a pixel having the average of the four pixels. We define the *n*-th *subscale* of an
<figure><img src="image_45.png" /></figure>

FIGURE 18. Noise curves overlapping for 1/30 s, 1/250 s, 1/400 s and 1/640 s exposure times for the Canon camera with ISO 1250. Up: raw image, after demosaicing and after white balance. Down: after the application of the tone curve, and after JPEG compression.

<figure><img src="image_46.png" /></figure>

FIGURE 19. Noise curves overlapping for 1/30 s, 1/250 s, 1/400 s and 1/640 s exposure times for the Nikon camera with ISO 1250. Up: raw image, after demosaicing and after white balance. Down: after the application of the tone curve, and after JPEG compression.
**Algorithm 13 Compute the average of overlapping noise curves.**

1: **MEAN_OVERLAPPING** - Computes the average of the standard deviations coming from noise curves of different exposure times with the same ISO speed. **Input** $X_t, X_y$: noise curves vectors (intensity and standard deviation) from different exposure times and the same ISO speed. **Output** **g**: **u** with one quarter of its medium/high frequencies set to zero.

2: $x\_range = \text{arange}(0.0, 4096, 0.1)$ ▷ The whole intensity range with 12 bits and step 0.1

3: $Y\_mean = \text{zeros}(\text{num\_channels}, \text{len}(x\_range))$ ▷ Output array with the averaged stds.

4: **for** ch = 0 ... num\_channels **do**

5: $x\_idx = 0$

6: **for** $x \in x\_range$ **do**

7: $\text{num} = 0$

8: $\text{mean} = 0.0$

9: **for** $t \in \text{range}(4)$ **do** ▷ 4 is the number of captors at the CFA

10: **if** $\min(X_t[t, \text{ch}, :]) \le x \le \max(X_t[t, \text{ch}, :])$ **then** ▷ If intensity belongs to any curve

11: $\text{idx} = \text{argmin}_i (|X_t[t, \text{ch}, i] - x|)$

12: **while** $X_t[t, \text{ch}, \text{idx}] > x$ **do**

13: $\text{idx} = \text{idx} - 1$

14: **end while**

15: **if** $\text{idx} == \text{num\_bins} - 1$ **then**

16: $\text{idx} = \text{idx} - 1$

17: **end if**

▷ The coordinates of the control points in between

18: $x_1, x_2 = X_t[t, \text{ch}, \text{idx}]$, $X_t[t, \text{ch}, \text{idx} + 1]$

19: $y_1, y_2 = Y_t[t, \text{ch}, \text{idx}]$, $Y_t[t, \text{ch}, \text{idx} + 1]$

20: $\text{mean} = \text{mean} + y_1 + (x - x_1) \times (y_2 - y_1)/(x_2 - x_1)$ ▷ Interpolation of the std.

21: $\text{num} = \text{num} + 1$

22: **end if**

23: **end for**

24: **if** $\text{num} > 0$ **then** ▷ If at least one overlapping noise curve

25: $\text{mean} = \text{mean} / \text{num}$

26: $Y\_mean[\text{ch}, x\_idx] = \text{mean}$ ▷ Store final averaged std.

27: **end if**

28: $x\_idx = x\_idx + 1$ ▷ Next intensity

29: **end for**

30: **end for**

image **u** as the result of applying *n* times operator **S** to **u**. If **u** is a discrete pure Gaussian noise image with standard deviation $\sigma$, then **S**(**u**) has standard deviation $\frac{\sigma}{2}$. Indeed, if a block **W** contains pixels $\{p_1, p_2, p_3, p_4\}$ each one with variance $\sigma^2$, the variance of the mean of **W** is $\text{Var}(\bar{\text{W}}) = \text{Var}(\frac{p_1+p_2+p_3+p_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}$. The standard deviation is $\text{Std}(\bar{\text{W}}) = \frac{\sigma}{2}$, that is, the noise is divided by two at each subscale.
<figure><img src="image_47.png" /></figure>

FIGURE 20. Mean noise curve for the Canon camera with ISO 1250. Since this noise curve was obtained from raw images, it can be used as the ground-truth of the noise of the camera for that ISO. On the right, a detail of the curve at the darkest pixels.

<figure><img src="image_48.png" /></figure>

FIGURE 21. Mean noise curve for the Canon camera with ISO 1600. Since this noise curve was obtained from raw images, it can be used as the ground-truth of the noise of the camera for that ISO. On the right, a detail of the curve at the darkest pixels.

In this chapter we compare the autocorrelation matrix of a 10 × 10 pixels block of an rectangle in one snapshot of the calibration pattern with ISO 1600 and exposure time 1/250s for both the Canon and Nikon cameras. The autocorrelation function has support 11 × 11 pixels.

Figure 24 shows the autocorrelation matrices with support 11 × 11 with the calibration pattern obtained using the Canon camera. From up to bottom, the autocorrelation functions obtained at the raw, demosaiced, white balance, gamma correction and JPEG compression accumulated steps of the chain. From left to right, the measurements at the first, second, third and fourth subscales of the image. Note that the scale of the color bar changes at each processing step. Figure 25 shows the same for the Nikon camera. The following discussing is applicable to both cameras. At
<figure><img src="image_49.png" /><figcaption>FIGURE 22. Mean noise curve for the Nikon camera with ISO 1250. Since this noise curve was obtained from raw images, it can be used as the ground-truth of the noise of the camera for that ISO. On the right, a detail of the curve at the darkest pixels.</figcaption></figure>

<figure><img src="image_50.png" /><figcaption>FIGURE 23. Mean noise curve for the Nikon camera with ISO 1600. Since this noise curve was obtained from raw images, it can be used as the ground-truth of the noise of the camera for that ISO. On the right, a detail of the curve at the darkest pixels.</figcaption></figure>

the raw image step, it can be observed that the noise is not correlated at any scale. The central value at [0, 0] is the only significant value and corresponds to the correlation of the signal with itself without shifting. Since a subscale just consists of computing the average of each disjoint set of 2 × 2, the autocorrelation values decrease as the number of scales increase. This can be observed no matter the step of the camera processing chain. The effect of demosaicing (by the Adams-Hamilton algorithm) is clearly noticeable, as it correlates strongly the pixels at distance 1 and slightly the pixels at distance $\sqrt{2}$. As expected, this effect is minimized at the second scale and beyond, due to the pixel averaging. Finally, the JPEG step result is similar to that obtained at the gamma correction. The value at [0, 0] is lower, since the JPEG algorithm averages the red,
green and blue color channels in order to turn the *RGB* color space into $Y'CBC_R$. Note that the autocorrelation pattern observed after the demosaicing step will remain along all the rest of the camera processing chain steps. The origin of the correlation observed in the JPEG image is the demosaicing step.

<figure><img src="image_51.png" /><figcaption>FIGURE 24. Autocorrelation functions for the Canon camera, ISO 1600 and exposure time 1/250s. From up to bottom: results of the steps raw, demosaicing, white balance, gamma correction and JPEG correction. From left to right: first, second, third and fourth scales. Note that the scale of the color bar changes at each processing step.</figcaption></figure>

5. Conclusions

Under the same ISO speed and exposure time conditions, the noise curves obtained by different cameras differ. Some cameras do not pre-process at all the data acquired at the CCD or CMOS detector and therefore it would be possible to assume a Poisson model for the noise. However, in other cameras the data at the raw image has been already altered in an unknown way, thus
<figure><img src="image_52.png" /><figcaption>FIGURE 25. Autocorrelation functions for the Nikon camera, ISO 1600 and exposure time 1/250s. From up to bottom: results of the steps raw, demosaicing, white balance, gamma correction and JPEG correction. From left to right: first, second, third and fourth scales. Note that the scale of the color bar changes at each processing step.</figcaption></figure>

making impossible to assume any model. Therefore, it is preferable to use non-parametric models that directly estimate a noise curve from the image itself with assuming a predefined model. The only necessary assumption of such estimation methods is that the noise is signal and frequency dependent. The frequency dependency is justified because the final JPEG encoding is equivalent to low-pass filtering the image (besides adding block artifacts and shot noise). Since the quantization matrix of the JPEG encoder depends on the unknown quality parameter, the only way to find a proper model for the noise is again to estimate the noise at each frequency, avoiding any predefined model for the JPEG filter. In general, denoising algorithms need an accurate estimation of the noise to properly denoise an image. Since the exact transformations that have been applied to the
noisy image are unknown, assuming a model is too risky and it is preferable to get a profile of the
noise depending of the intensity and the frequency from the noisy image itself.

In raw images both the spatial and temporal estimations matched accurately, but after the
next step (demosaicing), they did not. The explanation for this is that demosaicing correlates the
noise (thus, it becomes frequency-dependent) and the signal-dependent noise model is not enough
to measure the noise in these kind of real images. Most of the signal-dependent noise estimators
assume that the variance at any frequency is the same, and thus they compute an average of the
variances measured at the high-frequencies (where the noise dominates over the geometry of the
image). After demosaicing, this does not hold. In Chapter 4 a method for estimating both intensity
and frequency dependent noise is presented.
# Multiscale estimation of intensity and frequency dependent noise

The camera calibration and the image processing chain that generated a given image are generally no more accessible to the receiver. This happens for example with scanned photographs and for most JPEG images. These images have undergone various nonlinear contrast changes and linear and nonlinear filters. To deal with remnant noise in such images, we introduce a general non parametric intensity and frequency dependent noise model. We demonstrate by simulated and real experiments that this model, which requires the estimation of more than 1000 parameters, performs an efficient noise estimation.

The proposed noise model is a patch model. Its estimation can therefore be used as a preliminary step to any patch-based denoising method. Our noise estimation method introduces several new tools for performing this complex estimation. One of them is a new sparse patch distance function permitting to detect efficiently pure noise patches. We also show a new way to avoid the bias of noise estimation methods based on the use of low variance patches.

A validation of the noise model and of its estimation method is obtained by comparing its results to ground-truth noise curves for both raw and JPEG-encoded images, and by visual inspection of the denoising results of real images. A fair comparison to the state of the art is also performed.

This chapter presents a complex noise model that makes it possible to measure the noise even in images where the noise is highly correlated, as in the case of JPEG-encoded images. It is the final step of the study on noise estimation presented in this thesis, that begins with the simple but absolutely unrealistic simulated homoscedastic white Gaussian noise (Chapter 1), later the signal-dependent noise, which is only useful for raw images (Chapter 2), and finally the intensity-frequency dependent noise that is presented in this chapter. Chapter 8 describes a new denoising algorithm able to denoise even JPEG-encoded images using the noise estimation algorithm presented here, where medium and specially low-frequency noise remains after compression.

## 1. Introduction

The noise initially present in a raw digital image is transformed at each step of the processing chain of the camera. When acquired at the focal plane in a color filter array, the noise is Poisson distributed, intensity-dependent and frequency-independent. Yet the image at the Color Filter Array (CFA) is possibly saturated, which infringes the simple linear dependency of the noise variance with the intensity [7]. Even without saturation, the variance of the noise may not follow the linear
model, depending on the characteristics of the detector [8]. At the very first step of the camera pro-
cessing chain a demosaicing algorithm [81, 80] must be applied to get a color image from the raw
mosaic acquired at the CFA. This causes the noise to be spatially correlated. It therefore becomes
frequency-dependent. The colored noise caused by the demosaicing undergoes further nonlinear
transforms such as the white balance and a gamma-correction. Finally, JPEG-encoding [79] accen-
tuates the frequency dependence, as JPEG encoding applies a frequency-dependent quantification
matrix to the coefficients of the 8 × 8 DCT-II blocks of the image. Therefore, the remaining noise in
a JPEG image is signal dependent and highly correlated. It generally contains little high frequency
noise, as the quantization removes the image high frequencies. But it too often still contains strong
noise at the low and medium DCT block frequencies. This annoying noise is hard to evaluate and
to remove.

Such noise characteristics are observed in modern digital images, but also in scans of old
photographs, which contain chemical noise. The assumption that their final observed noise is both
signal and frequency dependent (SFD) is clearly a minimal model. The purpose of this chapter is
to develop a method for estimating such SFD noise, and to validate it by comparing the estimated
results to the appropriate ground-truth.

Little has been written on SFD noise estimation from a single digital image. A method estimating a “JPEG compression history” from a single image can be found in [82]. The noise estimation method for JPEG images proposed in [52] estimates a signal dependent noise level which is not frequency dependent and therefore only gives a “noise level”. One of the most complete attempt to estimate a general noise model is contained in the blind denoising method [37], which estimates multiscale noise covariances for noise wavelet coefficients. This model is nevertheless not signal dependent. The recent method for estimating frequency dependent noise on patches in [83] is probably the closest to our endeavor. We will detail the points in common and the proposed extensions and improvements of this method. To the best of our knowledge, no method has proposed so far to estimate a general SFD noise patch model. The situation is nonetheless favorable, as most homoscedastic noise estimation algorithms are actually patch based [2, 3, 4, 5, 6, 1], and can therefore be adapted to measure SFD noise models on patches.

Our plan is this Chapter follows. Section 2 develops the principles of blind noise estimation,
defines the signal and frequency dependent (SFD) model, and explains progressively how to esti-
mate it. Section 3 details the proposed algorithm. Section 4 performs a comparison of the method
with the current state of the art. Section 5 is the core of the chapter. It validates the sufficiency of
the SFD model for JPEG images by calculating a ground truth SFD model and checking that it is
indeed obtained by the algorithm. This section also performs a final consistency check by display-
ing denoising results obtained with a multiscale version of the NL-Bayes algorithm [84, 22]. To
that aim, it compares the estimated noise model before and after denoising. Section 5.1 contains
the conclusions.
## 2. Blind noise estimation principles

This chapter is a progressive presentation of patch based noise modeling. It gives a definition of the SFD noise model and discusses how to estimate it with minimal bias from well chosen blocks extracted from the image itself.

2.0.0.1. The search for blocks with low variance. The first principle adopted by most recent methods is to look for small image blocks that contain only noise. Indeed, all image blocks contain a sum of signal and noise. Being independent, the variances of noise and signal add. This leads to the conclusion that blocks with minimal variance extracted from the image are likely to contain no signal, and therefore only noise. In the algorithms adopting this strategy, all image blocks are ordered by their variance, and the noise estimation is performed on their lowest quantile (typically taking the 0.5% blocks with lowest variance).

2.0.0.2. Compensating the bias. Computing a median of these variances gives an estimate of the noise variance. Yet, this estimation is biased. Indeed, if the image were a pure white noise realization, the above method would estimate the noise variance on the blocks with lowest variance, leading to an obvious underestimation. This bias can nevertheless be compensated by a multiplicative correction factor learnt on a white noise image of the same size.

2.0.0.3. Using instead block differences. A recent method [83] has proposed a clever way to extend the low variance block method by involving the image self similarity. The idea is to associate with each block its most similar block in a neighborhood. Then, assuming that this similarity is essentially caused by the signal, the difference of both blocks is assumed to give a pure noise block, with twice the variance of the original block. Then, again, only a small quantile of those "noise blocks" should be retained for the final noise estimation. In practice, most of the selected blocks correspond to flat zones (as we shall see in the comments of Figure 1). If the image lacks flat zones, this selection of block differences may nevertheless also contain differences of non flat patches with similar geometry.

2.0.0.4. Dealing with signal dependency. The emission of photons by a physical body is a random process that can be modeled with a Poisson distribution, for which the mean is equal to the variance. Thus, the image formed at the CFA contains noise that depends on the intensity of the underlying image. This intensity dependence of the noise model begins at the very first step (acquisition at the detector: the raw image) and remains until the last step of the camera processing chain (JPEG encoding). The value of the tabulated gamma correction function is generally unknown. Even when this information is available, the CCD or CMOS detector do not necessarily follow a simple linear relation intensity/variance [8] when acquired at the detector.

Therefore, the noise estimation algorithm must estimate intensity-dependent noise. A famous alternative is to transform the data to turn it into homoscedastic noise via an Anscombe transform [85, 43]. Yet an Anscombe transform is only possible if we deal with raw images. In the general setting of a signal dependency that can be different for every frequency, there is no other way to estimate the signal dependent noise model than dividing the set of blocks into disjoint bins,
each for a different intensity and to estimate a separate frequency-dependent noise model on each intensity bin.

Fortunately, as recalled in [1, 9], it is possible to adapt most patch-based homoscedastic noise estimation methods [2] [3] [4] [5] [17] [15] to measure intensity-dependent noise, by simply splitting the list of input blocks into sets of blocks disjoint in mean intensity (bins).

2.0.0.5. Dealing with frequency dependency: the DCT diagonal assumption. According to the above considerations, the SFD noise estimation should start by finding a set of pure noise patches for each intensity bin. Our main assumption will be that the unknown linear and nonlinear transforms that have been applied to the image can be approximated by a diagonal operator on the DCT patch coefficients. There are two arguments in favor of this diagonality. First, it is easily checked that every linear real symmetric filter applied to an image is a diagonal operator on the DCT transform. This observation is actually exact for a global DCT transform, but remains approximately true for the (local) block DCT. Second, we mentioned that JPEG 1985 also is a diagonal (nonlinear) operator on the DCT coefficients. The demosaicing operation itself is more complex, but a good demosaicing algorithm must avoid creating structure in noise, and must therefore be close to a linear isotropic filter, which again is a diagonal operator. Nevertheless, these arguments are no end proofs. The DCT diagonal assumption can only be conjectured and must be verified empirically, as we shall do in the experimental chapters. Under this assumption, a patch noise model is described by the variances of its DCT coefficients.

2.0.0.6. Definition of the SFD noise model and its evaluation. The proposed signal and frequency dependent (SFD) noise model follows from the above observations, that also hint at an estimation algorithm. It is enough to find sufficiently many noise patches in a given image, and to apply to them a DCT before measuring the variance of their DCT coefficients. In that way the SFD model is defined (and can be empirically estimated) for each patch size w as a map:

$$ (i, j, b) \in [[0, w-1]]^2 \times [[0, B-1]] \rightarrow \sigma(i, j, b), $$

where (i, j) is the DCT frequency, w the block size, B the number of color level bins, and σ(i, j, b) the observed noise standard deviation for this particular frequency and bin. This model must be estimated independently for each color channel. The SFD model has therefore many parameters. For example for a (minimal) block size w = 4, B = 25 bins and a three channels, 16 × 25 × 3 = 1200 parameters must be estimated. This also explains why a more sophisticated model can hardly be envisaged. It might for example seem natural to estimate a whole covariance matrix for the noise patches, but this would require many more samples than those available in a single image. (Nevertheless, this extension might be considered for video noise estimation.)

The method proposed in [83] is the closest to the above considerations. It computes a frequency dependent noise model for the purpose of estimating homoscedastic frequency dependent noise. We shall now propose several extensions and improvements to make it fit for accurate SFD estimation, when many more free parameters are being estimated. Our goal is indeed to extend the frequency dependent model to a general enough model permitting to cope with most images.
We saw that introducing a signal dependence in the noise model is a necessary extension. The next improvement will be to define an unbiased block distance as explained in the next paragraph.

2.0.0.7. A sparse semi-distance between blocks. The frequency dependent noise estimation of [83] proposed to generate pure noise blocks by pairing similar blocks and making their difference to eliminate the signal. Their distance is evaluated in this method on a random choice of 50% of the block DCT coefficients. We shall explain the bias caused by this method and try to alleviate it. We propose to compare blocks with a sparse semi-distance which we introduce now. Given an image DCT block $\tilde{m}_1$, the idea it to only use its $w^2/4$ coefficients with largest absolute value for comparing them with the coefficients at the same frequencies in another block $\tilde{m}_2$. According to the transform thresholding principle [73], these high-energy coefficients belong with high probability to the signal and not to the noise. Indeed, the $w \times w$ DCT-II matrices of the blocks in a natural and not noisy image are generally sparse, and the signal is mainly contained in some low and medium frequency coefficients. Thus, we can assume that the semi-distance<fnref n="1" /> computed with that set of coefficients is not guided by the noise in the patch, but by the image geometry.

To compute the sparse semi-distance between a reference block $\tilde{m}_1$ and a candidate $\tilde{m}_2$, the absolute values of the DCT coefficients $A[i, j], [i, j] \in [0, w-1]^2$ of $\tilde{m}_1$ are sorted. The decreasing sorting indices are stored in $Q$. The sparse semi-distance is then computed based on the $w^2/4$ DCT

**Algorithm 14** Sparse distance algorithm

1: **Input** : $\tilde{m}_1$: patch of size $w \times w$

2: **Input** : $\tilde{m}_2$: patch of size $w \times w$

3: **Output** : sparse distance between $\tilde{m}_1$ and $\tilde{m}_2$.

4: Compute the absolute value of the DCT-II coefficients of $\tilde{m}_1$, $A[i, j] = \{|\tilde{m}_1[i, j]| : [i, j] \in [0, w-1]^2\}$.

5: Sort $A$ from the highest to the lowest value and put the sorting indices in $Q$; $Q = argsort(A)$.

6: Compute the SD on the first $w^2/4$ sorting indices in $Q$: $SD_{\tilde{m}_1, \tilde{m}_2} = \sum_{q=0}^{w^2/4-1} (\tilde{m}_1[Q_q] - \tilde{m}_2[Q_q])^2$.

coefficients in $\tilde{m}_1$ with the largest absolute value. This computation is summarized in Algorithm 14.

2.0.0.8. Should noise be estimated on block differences? The proposal made in [83] to choose the differences $\frac{1}{\sqrt{2}}(\tilde{m}_1 - \tilde{m}_2)$ as noise block samples is quite tempting. Indeed, the sparse semi-distance ensures that in this difference the part of the signal that was contained jointly in $\tilde{m}_1$ and $\tilde{m}_2$ has been canceled, thus giving a pure noise sample. Nevertheless, we found that it was better not to operate this subtraction, and that the noise estimation based on the blocks $\tilde{m}_1$ such that $\tilde{m}_1 - \tilde{m}_2$ is minimal was more accurate than the estimation based on the normalized differences $\frac{1}{\sqrt{2}}(\tilde{m}_1 - \tilde{m}_2)$. Experimental evidence confirms that the subtraction is advised only in the case where the image is known a priori to be made mainly of textures and to lack flat or smooth zones. We can anticipate the following explanation. Consider an image of pure Gaussian noise of

<footnote n="1">¹We call it semi-distance and not distance, as it is not symmetric and does not satisfy the triangle inequality.</footnote>
variance $\sigma^2$. Given a block $\tilde{m}_1$, and its most resembling block $\tilde{m}_2$, the similar DCT coefficients in both $\tilde{m}_1$ and $\tilde{m}_2$ are anyway due to noise. Thus the (small) block $\tilde{m}_1 - \tilde{m}_2$ generally loses low frequency coefficients belonging to the noise. In this situation, the estimated noise can be drastically distorted and underestimated. On the other hand, it turns out that the blocks $\tilde{m}_1$ such that the distance $\tilde{m}_1 - \tilde{m}_2$ is minimal are a better choice than blocks with minimal variance. Indeed, the minimal difference criterion is less biased than the minimal variance criterion. This empirical statement will be checked experimentally in Section 4. This explains why we shall select the best noise blocks $\tilde{m}_1$ based on their difference with the most resembling block, but keep $\tilde{m}_1$ for the final noise estimation.

### 2.0.0.9. The multiscale approach to estimate low frequency noise.
The selection of noise blocks faces another dilemma: on the one hand it is much easier to find small (typically $4 \times 4$) blocks containing only noise, than larger (e.g. $8 \times 8$) blocks. Yet, small blocks do not permit to estimate the noise low frequencies. Such low frequencies are prominent in JPEG images because of the demosaicing, which creates sometimes long range correlation, and because of the JPEG compression itself. So noise can appear in large spots, as shown in Figure 9. This image is the result of convolving an image of pure Gaussian noise with mean 127 and $\sigma = 50$ with the kernel $G$ in Equation 4. A multiscale approach solves the dilemma. Defining the input noisy image as the image at the first scale, a second scale image can be obtained by a 2-subsampling. Estimating again noise in this subsampled image permits to catch the noise low frequencies.

### 2.0.0.10. Denoising.
Once the STD of the noise at each intensity $b$ and each frequency $[i, j] \in [0, w-1]^2$, $[i, j] \neq [0, 0]$ is known, it is possible to fully characterize the noise by its covariance matrix. Thus, it is possible to denoise the image by obtaining the covariance matrix of the noise at each scale, and then denoising each scale using the obtained noise profile. Since the number of samples is divided by 4 after each subscaling, the number of pixels of the input image is a limiting factor (highly correlated noise cannot be measured in small images). The Nonlocal Bayes [84, 22] algorithm will be used for that purpose, as it only requires a knowledge of the noise covariance matrix.

## 3. Noise estimation algorithm

Our proposed SFD noise estimator (Algorithm 15) follows from the considerations of the preceding chapter.

In the first step, $w \times w$ (typically, $w = 8$ or $w = 4$) overlapping blocks are extracted from the input noisy image, and their 2D orthonormal DCT-II is computed. For each of these DCT blocks $\tilde{m}_1$, the most similar block $\tilde{m}_2$ in a fixed neighborhood is found. The candidate blocks are at a spatial range belonging to the interval $[r1, r2] = [4, 14]$ (the spatial distance is defined as the maximum between the horizontal and vertical distances of the block centers). Given a block $\tilde{m}_1$ at spatial position $[x_1, y_1] \in [0, N_x - w] \times [0, N_y - w]$, the candidate blocks are searched only at positions $[x_2, y_2] : x_2 > x_1 \land y_2 > y_1$ to avoid repeating matching pairs. The similarity between $\tilde{m}_1$ and any other block in its neighborhood is evaluated with the Sparse Distance (SD) function
(Algorithm 14), designed to avoid the interference of the noise in the similarity evaluation. For each $\tilde{m}_1$ in the image, a corresponding most similar block $\tilde{m}_2$ is found, and $\tilde{m}_1$ and the sparse distance $\text{SD}_{\tilde{m}_1, \tilde{m}_2}$ are stored in list $\mathbf{L}$.

<caption>**Algorithm 15** Noise estimation algorithm</caption>

1: **Input** : Noisy image $\mathbf{u}$ of size $N_x \times N_y$ pixels.
2: **Input** : $w \times w$ size of the block in pixels.
3: **Output** : SFD noise curve $\tilde{\sigma}$.
4: Extract from the input image $\mathbf{u}$ of size $N_x \times N_y$ all possible $M = (N_x - w + 1)(N_y - w + 1)$ overlapping $w \times w$ blocks $\mathbf{B}_k$ and compute their 2D orthonormal DCT-II, $\tilde{\mathbf{B}}_k$, $k \in [0, M - 1]$.
5: Set $\mathbf{L} = \emptyset$
6: **for** each DCT block $\tilde{m}_1 \in \tilde{\mathbf{B}}$, **do**
7: Find the block $\tilde{m}_2$ that minimizes the sparse distance $\text{SD}_{\tilde{m}_1, \tilde{m}_2}$ (see Algorithm 14). Consider only $\tilde{m}_2$ blocks whose horizontal and vertical distance with respect to $\tilde{m}_1$ belongs to the interval $[r_1, r_2] = [4, 14]$.
8: Extract from $\tilde{m}_1$ the mean of $m_1$.
9: Add block $\tilde{m}_1$ and the associated sparse distance ($[\tilde{m}_1, \text{SD}_{\tilde{m}_1, \tilde{m}_2}]$) to list $\mathbf{L}$.
10: **end for**
11: Classify the elements of list $\mathbf{L}$ into disjoint bins according to the mean intensity of the blocks $[1, 75]$. Each bin contains (with the exception of the last) 42000 DCT blocks.
12: **for** each bin, **do**
13: Obtain the set $S_p$ made by those DCT blocks inside the current bin whose SD is below the $p$-quantile, with $p = 0.005$.
14: Assign to the current bin the intensity $I$ (Equation 3).
15: **for** each frequency $[i, j]$ with $[i, j] \in [0, w - 1]^2$, $[i, j] \neq [0, 0]$, **do**
16: Compute the (biased) STD of the noise at the current bin and frequency $[i, j]$ using the MAD estimator (Equation 1).
17: Correct the biased $\hat{\sigma}[I][i, j]$ and obtain the final STD estimate (Equation 2).
18: **end for**
19: **end for**

In order to estimate intensity-dependent noise, the set of blocks in $\mathbf{L}$ is split according to the mean intensity of the blocks. Each division of $\mathbf{L}$ is called a *bin* and contains only the blocks (and their associated SD) whose mean intensity belongs to a certain range. The bins are disjoint in intensity and their union gives the whole intensity range of the image.

At each bin the list $S_p$ is filled in with the blocks whose associated SD is below the $p$-quantile, with $p = 0.005$. The value of $p$ must be kept small, to ensure that $S_p$ contains blocks for which it was possible to find another block with minimal SD between them.

In step 9, block $\tilde{m}_1$ and the sparse distance $\text{SD}_{\tilde{m}_1, \tilde{m}_2}$ are stored in list $\mathbf{L}$.
Finally, for the DCT blocks in $S_p$, the standard deviation (STD) according to each frequency $[i, j] \in [0, w-1]^2$, $[i, j] \neq [0, 0]$ is computed using the MAD estimator (Equation 1).

$$ (1) \qquad \begin{aligned} \hat{\sigma}[I][i, j] &= \text{MAD}(S_p) \\ &= \text{median} \left( \left| \tilde{n}[i, j] - \text{median} (\tilde{m}[i, j]) \right| \right). \end{aligned} $$

Equation 2 gives the correction factor for the STD depending on the size of the blocks, for $p = 0.005$. A correction of the STD is needed because MAD is a biased estimator of the STD and also because the available number of coefficients to compute the (sample) variance is finite and thus biased. To obtain the correction factors, we added simulated homoscedastic noise of STD $\sigma = 5$ to a synthetic image of a calibration pattern with large flat zones of several different grayscale intensities. The biased STD $\hat{\sigma}$ was estimated with our algorithm and compared with $\sigma = 5$. The ratio $\sigma/\hat{\sigma}$ gives the correction factor. In Equation (2) we only give corrections for $w = 4$ and $w = 8$, but of course more correction factors may be tabulated if needed.

$$ (2) \qquad \tilde{\sigma}[I][i, j] = \begin{cases} 1.775 \times \hat{\sigma}[I][i, j] & \text{if } w = 4; \\ 1.677 \times \hat{\sigma}[I][i, j] & \text{if } w = 8. \end{cases} $$

The corresponding intensity $I$ is computed with Equation (3), as the median of the mean intensities under the $p$-quantile of blocks SD.

$$ (3) \qquad I = \text{median} \left( \frac{\tilde{m}[0, 0]}{w} \right) $$

The size of the block depends on the application of the noise estimation. A small block, say $4 \times 4$ ($w = 4$) is preferable for highly textured images, as the probability of capturing flat zones decreases with the size of the block. However, a larger block gives a more accurate frequency information about the noise and is more consistent with JPEG compression. As we mentioned the solution to this dilemma is to perform a two scale estimation with small blocks.

## 4. Comparison

In this chapter we compare our proposed noise estimation method with the best state-of-the-art paper on estimation of correlated noise [83]. This method is designed to estimate homoscedastic noise, but since it is patch-based, it is relatively easy to adapt it to intensity-dependent noise. As explained in Section 2, being able to estimate intensity-dependent noise is a mandatory step for any noise estimation methods that intends to be used with real images, where the noise is highly correlated after the demosaicing step.

We introduced the sparse distance function to avoid the interference of noise in the similarity measure between two noisy blocks. Figure 1 shows the blocks selected by Algorithm 15 using the sparse similarity function of Algorithm 14. If the image contains flat or smooth zones, the algorithm
<figure><img src="image_53.png" /><figcaption>FIGURE 1. Blocks selected by Algorithm 15 using the sparse similarity function of Algorithm 14, using the traffic test image shown in Figure 3 after adding homogeneous noise of STD = 10, for the bins #0 (a), #2 (b), #4 (c), and #6 (d) (7 bins are used). In a highly textured image, where is difficult to find smooth zones, the sparse similarity function will choose patches with similar geometry.</figcaption></figure>

will prefer them to estimate the noise, since they give the maximum similarity. Finally, in [83] the difference between similar blocks is stored for subsequent noise estimation, but as explained in Section 2, it is better not to compute the difference, to prevent the possibility of choosing blocks that are similar because of the noise and thus giving underestimations.

Now we will pass to compare and discuss the influence of two decisions taken in the design of our proposed SFD noise estimation: the subtraction or not of similar blocks in the list of DCT blocks under the 0.005-quantile and the performance of the new similarity function (Algorithm 14). In order to measure the influence of textures in the performance of the compared noise estimation methods, we used a synthetic noise-free calibration pattern (see Fig.2, left). Since calibration patterns lack texture, it is easy to find flat zones for which any variation of the intensity is due to the noise. Thus most noise estimation methods are expected to perform optimally on such images. To simulate the effect of textures, we considered an image that combines both the calibration pattern and a noise-free image. For example Figure 2 (right) shows the noise-free image obtained by adding 80% of the intensity of the calibration pattern and 20% of a noise-free textured image traffic. As both combined images are noise-free, the result is still noise-free, but textured.

To show the influence of textures in the noise estimation depending on the method, we added simulated intensity-dependent noise of variance $\sigma^2 = 100+7u$ to the combination of the calibration pattern image with several noise-free images and then estimated the Root Mean Squared Error (RMSE) along all frequencies and intensity bins ($u$ is the pixel intensity of the combined image). The level of texture of the image analyzed was controlled by a parameter $\alpha$, the combination being
<figure><img src="image_54.png" /><figcaption>FIGURE 2. On the left, synthetic noise-free calibration pattern. On the right, the weighted sum of the calibration pattern with the noise-free test image *traffic* (with weights 0.8 for the calibration pattern and 0.2 for the *traffic* image. Since both combined images are noise-free, the result is still noise-free, but textured.</figcaption></figure>

<figure><img src="image_55.png" /><figcaption>FIGURE 3. Noise-free images used to measure the robustness of the methods to the presence of texture. Each image is 704 × 469 pixels.</figcaption></figure>

$αP + (1 − α)T$, where *P* is the calibration pattern and *T* the noise-free image that brings the texture. We used the four noise-free images shown in Figure 3.

Figure 4 shows the RMSEs obtained for the test images in Figure 3 using 8 × 8 blocks. In the horizontal axis, the value of $α ∈ [0, 1]$ (the texture level) and in the vertical axis, the RMSE along all frequencies and intensity bins.

We compared the adaptation to intensity-dependent noise of the Ponomarenko et al. method [83] (state-of-the-art in frequency-dependent noise estimation using patches) and our method, with two variants for each method: subtracting the blocks under the MSE quantile and not subtracting them (see the discussion of Section 2). It can be observed that, regardless of the texture level, it is almost systematically better not to subtract similar blocks before estimation. Compare the Ponomarenko method (labeled *Ponomarenko sub*) with the variant which does not subtract similar blocks (labeled *Ponomarenko no-sub*): the estimation obtained with subtraction has a lower RMSE
<figure><img src="image_56.png" /><figcaption>FIGURE 4. RMSEs obtained for the test images in Figure 3 using 8 × 8 blocks. In the horizontal axis, the value of $\alpha \in [0, 1]$ (the texture level) and in the vertical axis, the RMSE along all frequencies and intensity bins. We compared the adaptation to intensity-dependent noise of the Ponomarenko et al. method [83] and our method, and two variants for each method: subtracting the blocks under the low quantile and not subtracting them (see the discussion in Section 2). In general, avoiding the subtraction is the best option and it is what we propose in the presented method.</figcaption></figure>

than when performing subtraction. Only in the case of extremely highly-textured images (as in the case of the *bag* image), the subtraction brings to a lower RMSE. The plots also compare the proposed method (labeled *Proposed no-sub*) with the variant that subtracts the blocks (labeled *Proposed sub*). It can be clearly seen that the proposed method gives a lower RMSE thanks to the use of a better similarity function (Algorithm 14. The new similarity function is less affected by textures, since it measures the difference between blocks using only the coefficients with most energy. This coefficients are related to the geometry of the image and are not biased by other coefficients that carry information from the noise (as it happens in [83]). In the plot of the *Bag* image it can be observed that for high values of $\alpha$, the RMSE of the variant that does not subtract the blocks is similar to the variant subtracting the blocks.

Depending on the image (*computer*, *flowers*), the variant subtracting the blocks is slightly better than the variant without subtraction, for high values of $\alpha$, but only for highly-textured images, as shown in Figure 4. But the proposed method has a better overall RMSE.
5. Validation with ground truth JPEG noise

The SFD noise estimation method gives an estimation of the standard deviation (STD) of the noise that depends both on the intensity and frequency in a single image. It uses the observation of blocks at many spatial locations and will therefore be called in this chapter the *spatial* estimation, to match it with the ground truth *temporal* estimation. We validated the spatial estimation method by taking raw and JPEG photographs with a given camera. The value of the spatially estimated STD on a single image should match the ground-truth STD for that camera for the configured ISO speed [7], obtained from multiple frames. Note that with JPEG images we do not refer to the result of a mere JPEG compression, but to the result of the whole camera processing chain applied the raw image acquired at the focal plane of the camera on the CCD (or CMOS), including demosaicing, white balance, gamma correction, and the final JPEG encoding.

For that purpose, we took a sequence of pictures of the same scene taken with fixed camera position and constant lighting. Under these conditions, any variation of the intensity in any pixel through the sequence is only attributable to the effect of the noise. It is therefore possible to build a ground-truth noise curves for both raw and JPEG-encoded images, associating with each observed mean signal value the corresponding STD of its observed samples. Similarly, by *frequency noise curve* we mean a numerical function associating with each value of the block mean the STD of the DCT coefficient of the noise at that frequency. *Thus, there are as many noise curves as DCT coefficients.* To obtain such curves, instead of measuring the variation of the intensity of the pixels in a fixed position along the sequence, we consider all *M* overlapping *w* × *w* blocks in the image, compute their orthonormal DCT-II, and measure the variance at the intensity of the bin and frequency [*i*, *j*] ∈ [0, *w* − 1]², [*i*, *j*] ≠ [0, 0] along the coefficients of the blocks at the same spatial position and with varying image index.

The noise curve obtained this way for each *DCT frequency* will be called the *temporal estimation* and can be used as a ground-truth to compare with the spatial estimation. Even if a noise model for JPEG images has never been proposed in the literature, it is therefore possible to obtain reliable empirical curves for JPEG images. To obtain them, it suffices to JPEG-encode each image of the snapshot with the same quality parameter, and to apply the described procedure.

The objective of this chapter is to verify that the spatial STD measured at any frequency
[i, j] ∈ [0, w - 1]², [i, j] ≠ [0, 0] using the algorithm in Section 3 coincides with the STD of the
temporal series measured only at that frequency for all intensities. To build the temporal STD
noise curve we used 100 snapshots of the same calibration pattern, for both raw and JPEG-encoded
images. In principle, any image might be used to get the temporal STD of the noise, but it is
preferable to use an object with large flat regions of different gray levels, in order to avoid the
effect of textures in the temporal estimation. To be robust to outliers (the edges between the large
flat zones), we considered only the 0.05-quantile [17] of the STD estimations that was corrected
afterwards to obtain an unbiased estimate.

The procedure to compute the ground-truth curve for JPEG-encoded images for frequency
[i, j] ∈ [0, w − 1]², [i, j] ≠ [0, 0] from a set of H images is detailed in Algorithm 16.
**Algorithm 16** Algorithm to obtain the SFD ground-truth noise curve from a sequence of images

1: **Input**: Sequence of JPEG (or raw) images.

2: **Input**: *H* number of images in the sequence.

3: **Input**: $N_x \times N_y$ the size of each image in the sequence, in pixels.

4: **Output**: Ground-truth noise curve $\tilde{\sigma}$.

5: Set $M = (N_x - w + 1)(N_y - w + 1)$ the number of overlapping blocks.

6: Set $E1 = E2 = E3 = \text{array}(M)$. ▷ Array of size *M*

7: **for** each JPEG (or raw) image of the series of *H* images, **do**
8: Extract from the input image **u** of size $N_x \times N_y$ all possible *M* overlapping $w \times w$ blocks
$\tilde{\mathbf{B}}_k$ and compute their 2D orthonormal DCT-II, $\tilde{\mathbf{B}}_k, k \in [0, M-1]$.

9: **for** $k \in [0, M-1]$ **do**

10: $E1[k] = E1[k] + (\tilde{\mathbf{B}}_k[i, j])^2$.

11: $E2[k] = E2[k] + \tilde{\mathbf{B}}_k[i, j]$.

12: $E3[k] = E3[k] + \tilde{\mathbf{B}}_k[0, 0]/w$. ▷ The mean of $\tilde{\mathbf{B}}_k$

13: **end for**

14: $E1[k] = E1[k]/H$.

15: $E2[k] = E2[k]/H$.

16: $E3[k] = E3[k]/H$. ▷ Normalization

17: **end for**

18: Set $\mathbf{L} = \text{array}(M)$. ▷ Array of size *M*

19: **for** $k \in [0, M-1]$ **do**

20: Set $\mathbf{L}[k] = \left[ \frac{k}{k-1} (E1[k] - (E2[k]))^2 \right]^{1/2}$. ▷ STD

21: **end for**

22: Classify the elements of list **L** into disjoint bins [1, 75] according to the mean intensity $E3[k]$ of the blocks. Each bin contains (with the exception of the last) 42000 sample variance estimations.

23: **for** each bin *b*, **do**

24: Set **X** the means of the blocks in bin *b*.

25: Set **Y** the STDs of the blocks in bin *b*.

26: Get the 0.05-quantile of **Y** and set $\hat{\mu}$ the mean in **X** associated with it.

27: Assign the 0.05-quantile of **Y** to $\hat{\sigma}[\hat{\mu}][i, j]$.

28: Set $\hat{\sigma}[\hat{\mu}][i, j]$ the 0.05-quantile of **X**. Set $\hat{\mu}$ the mean at the quantile position in **X**.

29: Correct $\hat{\sigma}$ biased by the quantile and obtain the final control point of the ground-truth for intensity $\hat{\mu}$ and frequency $[i, j]$:

$$ \hat{\sigma}[\hat{\mu}][i, j] = 1.22 \times \hat{\sigma}[\hat{\mu}][i, j]. $$

30: **end for**
<figure><img src="image_57.png" /><figcaption>FIGURE 5. Comparison of the temporal (ground-truth, in green) and spatial STD (in red) for the Canon EOS 30D in raw images for ISO speed 1600 using blocks of 4 × 4 DCT coefficients. The temporal and spatial STD match despite some oscillation in the spatial estimation. The curve at the bottom right is the comparison between the averaged mean temporal STDs and the averaged mean spatial STDs (along all frequencies except DC), showing that on average both estimations match accurately.</figcaption></figure>

In the sequel, we compare the results of the spatial estimation to the ground-truth, for both raw and JPEG-encoded images taken with a Canon EOS 30D camera with exposure time $t = 1/30$s, ISO speed 1600, and blocks of $w × w$ DCT coefficients with $w = 4$ and $w = 8$. Figure 5 compares the temporal and the spatial STDs for raw images and Figure 7 shows the same for JPEG-encoded images with compression factor $Q = 92$ for $w = 4$. Only coefficients [1, 1], [2, 2], and [3, 3] are shown, but equivalent results were obtained with all 15 coefficients. Respectively, Figure 6 and Figure 8 for $w = 8$. Only coefficients [2, 2], [5, 5], and [7, 7] are shown for $w = 8$. The average of the estimations along all coefficients $[i, j] ∈ [0, w-1]^2$, $[i, j] ≠ [0, 0]$ is also given in both cases.

Despite small oscillations in the spatial estimation, there is an accurate match between both the spatial and temporal estimations in the case of raw and JPEG images. It can be concluded that the method is able to estimate reliably SFD noise.

This test was performed with snapshots of the calibration pattern, which is not textured and contains large flat areas whose spatial variations are caused mainly by the noise. Thus, the final validation must use real natural images compressed with JPEG. Since a proper noise model for JPEG encoding has not been already described, a visual comparison of the quality of the images before and after denoising using the frequency-by-frequency estimation given by the proposed method is needed. This comparison is performed in Section 5.1.
<figure><img src="image_58.png" /><figcaption>FIGURE 6. Comparison of the temporal (ground-truth, in green) and spatial STD (in red) for the Canon EOS 30D in raw images for ISO speed 1600 using blocks of 8 × 8 DCT coefficients. The temporal and spatial STD match despite some oscillation in the spatial estimation. The curve at the bottom right is the comparison between the averaged mean temporal STDs and the averaged mean spatial STDs (along all frequencies except DC), showing that on average both estimations match accurately.</figcaption></figure>

We also compared the accuracy of the proposed method by simulating colored noise and comparing the temporal STD (ground-truth) with the spatial estimation given by our algorithm for images of pure noise, frequency by frequency. To obtain the temporal STD, we created a list of 210 blocks of size 8 × 8 pixels made of simulated Gaussian noise of mean 127 and σ = 10 after applying a convolution with the discrete Gaussian kernel $G$ in Equation 4. Figure 9 shows a crop of the convolved noise image, where it can be observed that it contains spatial structure, as the noise is correlated because of the Gaussian convolution.

$$ (4) \qquad \mathcal{G} = \frac{1}{273} \begin{bmatrix} 1 & 4 & 7 & 4 & 1 \\ 4 & 16 & 26 & 16 & 4 \\ 7 & 26 & 41 & 26 & 7 \\ 4 & 16 & 26 & 16 & 4 \\ 1 & 4 & 7 & 4 & 1 \end{bmatrix}. $$

Table 1 compares for several frequencies (first column) the temporal STD obtained for the 210 8 × 8 blocks of pure noise (second column), the spatial STD estimation obtained by our method for pure noise after convolution with the discrete Gaussian kernel $G$ in Equation 4 (third column), and the spatial STD estimation given by our algorithm after adding homoscedastic Gaussian noise
<figure><img src="image_59.png" /><figcaption>FIGURE 7. Comparison of the temporal (ground-truth, in green) and spatial STD (in red) for the Canon EOS 30D in JPEG-encoded images with quality factor $Q = 92$ for ISO speed 1600 using blocks of 4 × 4 DCT coefficients. The curve at the bottom right is the comparison between the averaged temporal STDs and the averaged mean spatial STDs (along all frequencies except DC), showing that in average both estimations match.</figcaption></figure>

of $\sigma = 50$ to the noise-free test image *computer* (in Figure 3, top right) and then convolving it with *G* (fourth column). With a small margin of error, the proposed method is able to measure the STD of the noise for both pure noise and textured natural images. If nevertheless the STD of the noise is below 0.4, the method is unable to estimate it accurately and in some cases the MAD estimator gives negative values that the algorithm sets to zero afterwards.

**5.1. Denoising results.** Old photographs are particularly well adapted to be evaluated with a SFD noise model. Indeed, as such images involve two different successive acquisition systems, one chemical and one digital, the noise model is fully unknown and could not be learnt but from the image itself. And there is, of course, no ground-truth. Yet the visual inspection of the noise gives a very good hint at its independence from the (recovered) signal. To denoise JPEG digital images of old photographs, we used a modified version of the NL-Bayes algorithm [84] using the noise DCT coefficients estimated by our algorithm in Section 3. The details of the denoiser can be found in [22]. Of course, other patch-based denoisers [37, 36, 86, 87] might be used instead. For the denoising tests shown in this chapter, we estimated the noise at two scales to go deeper in low frequencies: a noise patch model was estimated at the finer scale and a second noise model was also computed after a Gaussian image zoom in.

Figure 10 shows denoising results for images with unknown noise model. In the first row, details of the noisy images; in the second row, details of the denoised images; in the third row,
<figure><img src="image_60.png" /><figcaption>FIGURE 8. Comparison of the temporal (ground-truth, in green) and spatial STD (in red) for the Canon EOS 30D in JPEG-encoded images with quality factor $Q = 92$ for ISO speed 1600 using blocks of 8 × 8 DCT coefficients. The curve at the bottom right is the comparison between the averaged temporal STDs and the averaged mean spatial STDs (along all frequencies except DC), showing that on average both estimations match.</figcaption></figure>

<figure><img src="image_61.png" /><figcaption>FIGURE 9. Crop of the image of pure Gaussian noise with mean 127 and $\sigma = 50$ after convolution with the kernel $G$ in Equation 4. The noise has spatial structure, as it gets correlated after Gaussian convolution.</figcaption></figure>

enhanced difference image (removed noise) between the noisy and denoised image. The color spots in the difference image and their random aspect at zones with the same intensity indicate that the denoising algorithm removed colored noise and, since details are kept at the denoised image, it can be concluded that the noise estimation was successful.

Figs. 11 and 12 show the noise curves corresponding to the low and high frequencies of the *Apollo* and *Kleiner* images for which a detail is shown at the bottom of Figure 10 using DCT blocks of 4 × 4 coefficients. A coefficient at frequency [*i*, *j*] ∈ [0, 3]² is assumed to belong to a “low-frequency” if *i* + *j* ≤ 2 and to a “high-frequency” otherwise. Image *Apollo* was taken in 18 May 1969 during the prelaunch tasks at the Launch Control Center’s Firing Room 3 at the Kennedy
<table><caption>TABLE 1. This table compares for several frequencies (first column), the temporal STD obtained for the 210 8 × 8 blocks of Gaussian noise of σ = 50 (second column), the spatial STD estimation obtained by our method for pure noise after convolution with the Gaussian kernel G in Equation 4 (third column), and the spatial STD estimation given by our algorithm after adding homoscedastic Gaussian noise of σ = 50 to the noise-free test image computer and then convolving it with G (fourth column). Both STD estimation in pure noise and in a textured natural image match with small error the temporal STD. For the image of pure noise a single bin is used and 7 bins for the computer image.</caption><thead><tr><th>Frequency</th><th>Temporal STD</th><th>Spatial (pure noise)</th><th>Spatial (computer)</th></tr></thead><tbody><tr><td>[1, 1]</td><td>31.45</td><td>29.20</td><td>30.72</td></tr><tr><td>[2, 2]</td><td>21.77</td><td>19.93</td><td>21.06</td></tr><tr><td>[3, 3]</td><td>10.03</td><td>11.16</td><td>10.73</td></tr><tr><td>[4, 4]</td><td>3.44</td><td>3.72</td><td>3.76</td></tr><tr><td>[5, 5]</td><td>0.73</td><td>0.62</td><td>0.61</td></tr><tr><td>[6, 6]</td><td>0.15</td><td>0</td><td>0</td></tr><tr><td>[7, 7]</td><td>0.13</td><td>0</td><td>0</td></tr></tbody></table>

Space Center<fnref n="2" /> and image Kleiner is a picture of an old tramway called “Kleiner Hecht” taken in 1998 in Dresden<fnref n="3" />. Both images contain large low-frequency noise and JPEG compression artifacts. We show the mean of the noise curves from the low-frequencies before (a) and after (b) denoising, where it can be observed that most of the noise remains at the low-frequencies of the image and that is strongly reduced after denoising. We also show the means for the high-frequencies before (c) and after (d) denoising. Since JPEG quantizes the value of the DCT coefficients at the high-frequencies (thus canceling most of them), the noise is clearly lower that what is observed at the low-frequencies, but nevertheless the noise could also be removed.

This chapter presented a non-parametric noise estimation method for SFD noise. It can be applied to images where the noise model is not available [8], as in the case of JPEG images. In general, a non-parametric estimation of SFD noise is needed for almost any kind of image, since the very first step in the camera processing chain is to demosaic the raw image, which correlates the noise. In general, the information of what algorithm used to demosaic the image is not made available by camera makers. However, the goal of blind noise estimation and denoising is to retrieve the noise model from the image itself, without relying at additional information such as metadata

<footnote n="2">²This file is in the public domain because it was solely created by NASA. NASA copyright policy states that “NASA material is not protected by copyright unless noted”. http://dayton.hq.nasa.gov/IMAGES/LARGE/GPN-2000-001849.jpg</footnote>

<footnote n="3">³Image licensed under the Creative Commons Attribution-Share Alike 3.0 Unported license, taken by Wikimedia Commons user Olaf1541. http://commons.wikimedia.org/wiki/File:Kleiner_hecht.jpg</footnote>
<figure><img src="image_62.png" /><figcaption>FIGURE 10. Denoising results of real images with unknown noise model. Up: detail of the noisy image. Middle: detail of the denoised image. Bottom: difference image (removed noise). The color spots in the difference image and its random geometry at zones with the same intensity indicate that the denoising algorithm removed colored noise and, since details are kept at the denoised image, it can be concluded that the noise estimation was accurate. In order to see clearly the low-frequency noise and the denoising results, the reader is invited to look at the images on the screen with a 400% zoom at least.</figcaption></figure>

in the file format. Of course, for old photographs the metadata information might not exist at all (analogic photography) or may have been lost after image manipulation and re-encoding.

To determine if two patches contain the same geometric information, a new similarity function based on the sparse distance between the patches introduced. It exploits the sparsity between patches to give a likelihood measure which is robust to noise. The use of the MAD estimator to compute the STD confers more robustness to the method.

The method was validated by showing that the STD obtained at the temporal series (the ground-truth) coincides with the spatial STD given by the proposed algorithm, for both raw and JPEG images. The denoising results show that indeed the noise estimator is able to give an accurate estimation, as low frequency noise is removed and most of the fine details are kept. Our next endeavor would be to include an impulse noise estimator to the non-parametric noise estimation model. Old photographs can indeed present this sort of noise.

Nevertheless this estimation algorithm cannot be applied to *any* noisy image. For example, it does not apply if the noise is space dependent (and not only signal dependent), as can be observed in some synthetic images. Another limitation for estimating highly-correlated noise is the size of
<figure><img src="image_63.png" /><figcaption>Original noisy image *Apollo* (a)</figcaption></figure>

<figure><img src="image_64.png" /><figcaption>Low-freqs., noisy (b)</figcaption></figure>

<figure><img src="image_65.png" /><figcaption>Low-freqs., denoised (c)</figcaption></figure>

<figure><img src="image_66.png" /><figcaption>High-freqs., noisy (d)</figcaption></figure>

<figure><img src="image_67.png" /><figcaption>High-freqs., denoised (e)</figcaption></figure>

FIGURE 11. Noise curves corresponding to the low and high frequencies of the *Apollo* image (a) for which a detail is shown at the second column of Figure 10 using DCT blocks of 4 × 4 coefficients. (b) and (c): mean noise curve at the low-frequencies before (b) and after (c) denoising. (d) and (e): mean noise curves at the high-frequencies before (d) and after (e) denoising. Most of the noise is at the low-frequencies.

the noisy image, since two scales of the image are needed and the number of available samples (pixels) and the second scale is one quarter of the pixels in the noisy image. If the image is small and contains highly correlated noise, it may not be possible to characterize it properly.

Chapter 8 describes a new denoising algorithm able to denoised even JPEG-encoded images, using the noise estimation algorithm presented in this chapter.
<figure><img src="image_68.png" /><figcaption>Original noisy image *Kleiner* (a)</figcaption></figure>

<figure><img src="image_69.png" /><figcaption>FIGURE 12. Noise curves corresponding to the low and high frequencies of the *Kleiner* image (a) for which a detail is shown at the last column of Figure 10 using DCT blocks of 4 × 4 coefficients. (b) and (c): mean noise curve at the low-frequencies before (b) and after (c) denoising. (d) and (e): mean noise curves at the high-frequencies before (d) and after (e) denoising. Most of the noise is at the low-frequencies. The color of each of the curves corresponds to each color channel of the image (red, green, and blue).</figcaption></figure>
<figure><img src="image_70.png" /></figure>
Part 2

PATCH DENOISING
Once the noise estimation models and methods have been studied in the first part of the dis-
sertation, in this second part we focus on the problem of denoising. In Chapter 5 we discuss the
Bayesian patch-based methods. This formalism permits to obtain a closed formula that given a
noisy patch in the image gives a denoised version from it.

Regardless of the block-based denoising method used, it is possible to improve dramatically the
results by using the three generic tools explained in Chapter 6: aggregation of estimates, iteration
and "oracle" filters, and colorspace transform. This remark is valid for all block-based denoising
principles and will be used in our final Noise Clinic (Chapter 8). The effect of this tools is demon-
strated with an example of a simple DCT transform threshold, where the three generic tools are
applied successively.

In Chapter 7 we give a detailed description and analysis of two denoising methods for which reliable faithful implementations are available: the Non-local means and Non-local Bayes methods. The algorithmic descriptions of both methods is provided. We shall extend this non-local Bayes algorithm in our final Chapter 8 on the Noise Clinic.

Chapter 8 presents a new algorithm for both blind noise estimation and denoising for images where the noise model is imperfectly known or even unknown. For such images, it is possible to estimate from a single image a noise model which is simultaneously signal and frequency dependent, as we showed in Chapter 4. We propose a multiscale denoising algorithm adapted to this broad noise model, based on the NL-Bayes denoising algorithm (Chapter 7). This brings to a blind denoising algorithm which we demonstrate on real JPEG images and on scans of old photographs for which the formation model is unknown. This algorithm is finally compared to the unique state of the art previous blind denoising method, based on Gaussian Scale Mixtures In The Wavelet Domain.
## Bayesian patch-based methods

Bayesian patch-based methods give an optimal formulation under the assumption that the patches similar to a given image patch follow a stochastic model. Given a noiseless patch P of the ideal noiseless image u, and Gaussian model and assuming a Gaussian noise model with independence of pixel values, the probability $P(\tilde{P}|P)$ (the probability of observing a noisy patch $\tilde{P}$ with the prior information that the noiseless version of the noisy patch is P) can be obtained directly. Thanks to Bayes' formula, it is possible to obtain $P(P|\tilde{P})$, the probability of the noiseless patch P knowing that the noisy patch $\tilde{P}$ has been observed. The goal is to deduce P from $\tilde{P}$ by finding P which maximizes $P(P|\tilde{P})$. This chapter concludes with a closed formula that gives a restored patch $\hat{P}_1$ from an observed noisy patch $\tilde{P}$, using the covariance matrix $C_{\tilde{P}}$ of $\tilde{P}$.

### 1. Obtaining a restored patch $\hat{P}_1$ from an observed noisy patch $\tilde{P}$

In this chapter, we will review the Bayesian patch-based method applied to denoising, under the additive white Gaussian noise model. The Bayesian principle is coupled with a Gaussian (or a mixture of Gaussians) model for noiseless patches. Given u the noiseless ideal image and $\tilde{u}$ the noisy image corrupted with Gaussian noise of standard deviation $\sigma$ so that

$$ (5) \qquad \tilde{u} = u + n, $$

the conditional distribution $P(\tilde{u} | u)$ is

$$ (6) \qquad P(\tilde{u} | u) = \frac{1}{(2\pi\sigma^2)^{\frac{M}{2}}} \exp\left(-\frac{\|u - \tilde{u}\|^2}{2\sigma^2}\right), $$

where M is the total number of pixels in the image.

In order to compute the probability of the original image given the degraded one, $P(u | \tilde{u})$, we need to introduce a prior on u. In the first models [30], this prior was a parametric image model describing the stochastic behavior of a patch around each pixel by a Markov random field, specified by its Gibbs distribution. A Gibbs distribution for an image u takes the form

$$ P(u) = \frac{1}{Z} \exp(-E(u)/T), $$

where Z and T are constants and E is called the energy function and writes

$$ E(u) = \sum_{C \in \mathcal{C}} V_C(u), $$
where $\mathcal{C}$ denotes the set of cliques associated to the image and $V_C$ is a potential function. The maximization of the a posteriori distribution writes by Bayes formula

$$ \operatorname{Arg\,max}_{u} \mathbb{P}(u | \tilde{u}) = \operatorname{Arg\,max}_{u} \mathbb{P}(\tilde{u} | u)\mathbb{P}(u), $$

which is equivalent to the minimization of $-\log \mathbb{P}(u | \tilde{u})$,

$$ \operatorname{Arg\,min}_{u} \|u - \tilde{u}\|^2 + \frac{2\sigma^2}{T} E(u). $$

This energy writes as a sum of local derivatives of pixels in the image, thus being equivalent to a classic Tikhonoff regularization [30, 88].

Recent Bayesian methods have abandoned as too simplistic the global patch models formulated by an a priori Gibbs energy. Instead, the methods build local non parametric patch models learnt from the image itself, usually as a local Gaussian model around each given patch, or as a Gaussian mixture. The term “patch model” is now preferred to the terms “neighborhood” or “clique” previously used for the Markov field methods. In the nonparametric models, the patches are larger, usually 8 × 8, while the cliques were often confined to 3 × 3 neighborhoods. Given a noiseless patch $P$ of $u$ with dimension $\kappa \times \kappa$, and $\tilde{P}$ an observed noisy version of $P$, the same model gives by the independence of noise pixel values

$$ (7) \qquad \mathbb{P}(\tilde{P}|P) = c \cdot \exp \left( - \frac{\|\tilde{P} - P\|^2}{2\sigma^2} \right) $$

where $P$ and $\tilde{P}$ are considered as vectors with $\kappa^2$ components, $||P||$ denotes the Euclidean norm of $P$, $\sigma^2$ the variance of the Gaussian noise, and $c$ is the normalizing constant. Knowing $\tilde{P}$, our goal is to deduce $P$ by maximizing $\mathbb{P}(P|\tilde{P})$. Using Bayes’ rule, we can compute this last conditional probability as

$$ (8) \qquad \mathbb{P}(P|\tilde{P}) = \frac{\mathbb{P}(\tilde{P}|P)\mathbb{P}(P)}{\mathbb{P}(\tilde{P})}. $$

$\tilde{P}$ being observed, this formula can in principle be used to deduce the patch $P$ maximizing the right term, viewed as a function of $P$. This is only possible if we have a probability model for $P$, and these models will be generally learnt from the image itself, or from a set of images. For example [89] applies a clustering method to the set of patches of a given image, and [90] applies it to a huge set of patches extracted from many images. Each cluster of patches is thereafter treated as a set of Gaussian samples. This permits to associate to each observed patch its likeliest cluster, and then to denoise it by a Bayesian estimation in this cluster. Another still more direct way to build a model for a given patch $\tilde{P}$ is to group the patches similar to $\tilde{P}$ in the image. Assuming that these similar patches are samples of a Gaussian vector yields a standard Bayesian restoration [47, 91]. We shall now discuss this particular case, where all observed patches are noisy.

Why Gaussian? As usual when we dispose of several observations but of no particular guess on the form of the probability density, a Gaussian model is adopted. In the case of the patches $Q$ similar to a given patch $P$, the Gaussian model has some pertinence, as it is assumed that many contingent random factors explain the difference between $Q$ and $P$: other details, texture, slight
lighting changes, shadows, etc. The Gaussian model in presence of a combination of many such random and independent factors is heuristically justified by the central limit theorem. Thus, for good or bad, assume that the patches $Q$ similar to $P$ follow a Gaussian model with (observable, empirical) covariance matrix $C_P$ and (observable, empirical) mean $\overline{P}$. This means that

$$ (9) \qquad \mathbb{P}(Q) = c \cdot \exp \left( \frac{-(Q - \overline{P})^t C_P^{-1} (Q - \overline{P})}{2} \right) $$

From (6) and (8) we obtain for each observed $\tilde{P}$ the following equivalence of problems:

$$ \begin{align*} \max_P \mathbb{P}(P|\tilde{P}) &\Leftrightarrow \max_P \mathbb{P}(\tilde{P}|P)\mathbb{P}(P) \\ &\Leftrightarrow \max_P \left[ \exp\left(\frac{-\|P - \tilde{P}\|^2}{2\sigma^2}\right) \exp\left(\frac{-(P - \overline{P})^t C_P^{-1} (P - \overline{P})}{2}\right) \right] \\ &\Leftrightarrow \min_P \left[ \frac{\|P - \tilde{P}\|^2}{\sigma^2} + (P - \overline{P})^t C_P^{-1} (P - \overline{P}) \right]. \end{align*} $$

This expression does not yield an algorithm. Indeed, the noiseless patch $P$ and the patches similar to $P$ are not observable. Nevertheless, we can observe the noisy version $\tilde{P}$ and compute the patches $\tilde{Q}$ similar to $\tilde{P}$. An empirical covariance matrix can therefore be obtained for the patches $\tilde{Q}$ similar to $\tilde{P}$. Furthermore, using (5) and the fact that $P$ and the noise $n$ are independent,

$$ (10) \qquad C_{\tilde{P}} = C_P + \sigma^2 \mathbf{I}; \quad E\tilde{Q} = \overline{P}. $$

Notice that these relations assume that we searched for patches similar to $\tilde{P}$ at a large enough distance, to include all patches similar to $P$, but not too large either, because otherwise it can contain outliers. Thus the safe strategy is to search similar patches in a distance slightly larger than the expected distance caused by noise. If the above estimates are correct, our MAP (maximum a posteriori estimation) problem finally boils down by (10) to the following feasible minimization problem:

$$ \max_P \mathbb{P}(P|\tilde{P}) \Leftrightarrow \min_P \left[ \frac{\|P - \tilde{P}\|^2}{\sigma^2} + (P - \tilde{P})^t (C_{\tilde{P}} - \sigma^2 \mathbf{I})^{-1} (P - \tilde{P}) \right]. $$

Differentiating this quadratic function with respect to $P$ and equating to zero yields

$$ P - \tilde{P} + \sigma^2 (C_{\tilde{P}} - \sigma^2 \mathbf{I})^{-1} (P - \tilde{P}) = 0. $$

Taking into account that $\mathbf{I} + \sigma^2(C_{\tilde{P}} - \sigma^2 \mathbf{I})^{-1} = (C_{\tilde{P}} - \sigma^2 \mathbf{I})^{-1}C_{\tilde{P}}$, this yields

$$ (C_{\tilde{P}} - \sigma^2 \mathbf{I})^{-1} C_{\tilde{P}} P = \tilde{P} + \sigma^2 (C_{\tilde{P}} - \sigma^2 \mathbf{I})^{-1} \overline{P}. $$

and therefore

$$ \begin{align*} P &= C_{\tilde{P}}^{-1} (C_{\tilde{P}} - \sigma^2 \mathbf{I}) \tilde{P} + \sigma^2 C_{\tilde{P}}^{-1} \overline{P} \\ &= \tilde{P} + \sigma^2 C_{\tilde{P}}^{-1} (\overline{P} - \tilde{P}) \\ &= \overline{P} + [\mathbf{I} - \sigma^2 C_{\tilde{P}}^{-1}] (\tilde{P} - \overline{P}) \\ &= \overline{P} + [C_{\tilde{P}} - \sigma^2 \mathbf{I}] C_{\tilde{P}}^{-1} (\tilde{P} - \overline{P}) \end{align*} $$
Thus we have proved that a restored patch $\hat{P}_1$ can be obtained from the observed patch $\tilde{P}$ by the one step estimation

$$ (11) \qquad \hat{P}_1 = \tilde{P} + [C_{\tilde{P}} - \sigma^2 \mathbf{I}] C_{\tilde{P}}^{-1} (\tilde{P} - \tilde{P}), $$

which resembles a local Wiener filter.

REMARK 1. It is easily deduced that the expected estimation error is

$$ E||P - \hat{P}_1||^2 = \mathrm{Tr} \left[ \left( C_{\tilde{P}}^{-1} + \frac{\mathbf{I}}{\sigma^2} \right)^{-1} \right]. $$

Chapter 7 will examine two deriving patch-based denoising algorithms from variants of (11). The first question when looking at this formula is obviously how the matrix $C_{\tilde{P}}$ can be learnt from the image itself. Each method proposes a different notion to learn the patch model.

Of course, other, non Gaussian, Bayesian models are possible, depending on the patch density assumption. For example [11] assumes a local exponential density model for the noisy data, and gives a convergence proof to the optimal (Bayes) least squares estimator as the amount of data increases.
Generic tools for noise reduction

This chapter describes three generic tools (aggregation of estimates, iteration and “oracle” filters, and colorspace transform) that permit to increase the performance of any denoising princi- ple. We shall illustrate them on DCT denoising. Starting from the application of a simple DCT transform threshold, the three generic tools will applied successively. We shall observe a dramatic improvement of the denoising performance. This remark is valid for all denoising principles and will be used in our final Noise Clinic (Chapter 8).

1. Aggregation of estimates

Aggregation techniques combine for any pixel a set of *m* possible estimates. If these estimates were independent and had equal variance, then a uniform average would reduce this estimator variance by a factor *m*. Such an aggregation strategy was the main proposition of the *translation invariant wavelet thresholding algorithm* [92]. This method denoises several translations of the image by a wavelet thresholding algorithm and averages these different estimates once the inverse translation has been applied to the denoised images.

An interesting case is when one is able to estimate the variance of the *m* estimators. Statistical arguments lead to attribute to each estimator a weight inversely proportional to its variance [93]. For most denoising methods the variance of the estimators is high near image edges. When applied without aggregation, the denoising methods leave visible “halos” of residual noise near edges. For example in the sliding window DCT method, patches containing edges have many large DCT coefficients which are kept by thresholding. In flat zones instead, most DCT coefficients are canceled and the noise is completely removed. The proposition of [94] is to use the aggregation for DCT denoising, approximating the variance of each estimated patch by the number of non zero coefficients after thresholding. In the online paper [95] one can test an implementation of DCT denoising. It actually uses an aggregation with uniform weights: “translation invariant DCT denoising is implemented by decomposing the image to sliding overlapping patches, calculating the DCT denoising in each patch, and then aggregating the denoised patches to the image averaging the overlapped pixels. The translation invariant DCT denoising significantly improves the denoising performance, typically from about 2 to 5 dB, and removes the block artifact”.

The same risk of “halo” occurs with non-aggregated NL-means (Chapter 7), since patches
containing edges have many less similar instances in the image than flat patches. Thus the non-local
averaging is made over less samples, and the final result keeps more noise near image edges. The
same phenomenon occurs with BM3D if the aggregation step is not applied [36]. As a consequence,
an aggregation step is applied in all patch-based denoising algorithms. This weighted aggregation favors, at each pixel near an edge, the estimates given by patches which contain the pixel but do not meet the edge.

Aggregation techniques aim at a superior noise reduction by increasing the number of values being averaged for obtaining the final estimate or selecting those estimates with lower variance. Kervrann et al [96] considered the whole Bias+Variance decomposition in order to also adapt the search zone of neighborhood filters or of NL-means. Since the bias term depends on the original image, it cannot be computed in practice, and Kervrann et al. proposed to minimize both bias and variance by choosing the smallest spatial neighborhood attaining a stable noise reduction.

Another type of aggregation technique considers the risk estimate rather than the variance to locally attribute more weight to the estimators with small risks. In [12], Van De Ville and Kocher give a closed form expression of Stein's Unbiased Estimator of the Risk (SURE) for NL-Means. (See also generalizations of the SURE estimator to the non-Gaussian case in [97].) The aim is to select globally the best bandwidth for a given image. In [98], Duval et al. also use the SURE technique for minimizing the risk by selecting locally the bandwidth. Deledalle et al. [13] apply the same technique for combining the results of NL-means with different window sizes and shapes. A similar treatment can be found in [11], but with the assumption of a local exponential density for the noisy patches.

## 2. Iteration and “oracle” filters

Iterative strategies to remove residual noise would drift from the initial image. Instead, a first step denoised image can be used to improve the reapplication of the denoising method to the initial noisy image. In a second step application of a denoising principle, the denoised DCT coefficients, or the patch distances, can be computed in the first step denoised image. They are an approximation to the true measurements that would be obtained from the noise-free image. Thus, the first step denoised image is used as an “oracle” for the second step.

For averaging filters such as neighborhood filters or NL-means, the image $u$ can be denoised in a first step by the method under consideration. This first step denoised image denoted by $\hat{u}_1$ is used for computing more accurate color distances between pixels. Thus, the second step neighborhood filter is

$$ YNF_{h,\rho}\tilde{u}(\mathbf{i}) = \frac{1}{C(\mathbf{i})} \sum_{j \in B_{\rho}(\mathbf{j})} \tilde{u}(\mathbf{j})e^{-\frac{|\tilde{u}_1(\mathbf{j}) - \tilde{u}_1(\mathbf{i})|^2}{h^2}}, $$

where $\tilde{u}$ is the observed noisy image and $\hat{u}_1$ the image previously denoised by Equation (12).

$$ (12) \qquad YNF_{h,\rho}\tilde{u}(\mathbf{i}) = \frac{1}{C(\mathbf{i})} \sum_{j \in B_{\rho}(\mathbf{i})} \tilde{u}(\mathbf{j})e^{-\frac{|\tilde{u}(\mathbf{i}) - \tilde{u}(\mathbf{j})|^2}{h^2}}. $$

Similarly, for linear transform Wiener-type methods, the image is first denoised by its classic definition, which amounts to approximate the oracle coefficients of Theorem 1 using the noisy ones.
In a second iteration, the coefficients of the denoised image approximate the true coefficients of the noise-free image. Thus, the second step filter following the first step (Equation 13) is

$$ \mathbf{D}\tilde{U} = \sum_{i} a(i) \langle \tilde{U}, G_i \rangle G_i, \quad \text{with} \quad a(i) = \frac{|\langle \hat{U}_1, G_i \rangle|^2}{|\langle \hat{U}_1, G_i \rangle|^2 + \sigma^2}, $$

where $\hat{U}_1$ is the denoised image by applying a first time the thresholding algorithm to the observed image $\tilde{U}$.

$$ (13) \qquad \mathbf{D}\tilde{U} = \sum_{i=1}^{M} a(i) \langle \tilde{U}, G_i \rangle G_i. $$

**THEOREM 1.** Operator $\mathbf{D}_{inf}$ minimizing the mean square error (MSE),

$$ \mathbf{D}_{inf} = \arg\min_{\mathbf{D}} E\{\|U - \mathbf{D}\tilde{U}\|^2\} $$

is given by the family $\{a(i)\}_i$, where

$$ (14) \qquad a(i) = \frac{|\langle U, G_i \rangle|^2}{|\langle U, G_i \rangle|^2 + \sigma^2}, $$

and the corresponding expected mean square error (MSE) is

$$ (15) \qquad E\{\|U - \mathbf{D}_{inf}\tilde{U}\|^2\} = \sum_{i=1}^{M} \frac{|\langle U, G_i \rangle|^2 \sigma^2}{|\langle U, G_i \rangle|^2 + \sigma^2}. $$

2.0.0.1. Alternatives and extensions: “twicing” and Bregman iterations. In the recent review paper [99], many denoising operators are formalized in a general linear framework, noting that they can be associated with a doubly stochastic diffusion matrix $W$ with nonnegative coefficients. For example in NL-means, this matrix is obtained by the symmetrization of the matrix of the NL-means weights $w_{\hat{P},\hat{Q}}$ defined in Algorithm 1. Unless it is optimal, as is the case with an ideal Wiener filter, the matrix $W$ associated with the denoising filter can be iterated. A study of MSE evolution with these iterations is proposed in [99] for several denoising operators, considering several different patch types (texture, edge, flat). Iteration is, however, different from the oracle iteration described above. In the oracle iteration, the matrix $W$ is changed at each step, using its better estimate given by the previously denoised image. One does not generally observe much improvement by iterating the oracle method more than once. [99] points out another generic tool, used at least for total variation denoising, the so-called “twicing”, term due to Tukey [100]. Instead of repeated applications of a filter, the idea is to process the residual obtained as the difference between the estimated image and the initial image. If the residuals contain some of the underlying signal, filtering them should recover part of it. The author shows that the Bregman iterations [101] used for improving total variation based denoising are a twicing and so is the matching pursuit method used in the K-SVD filter.
### 3. Dealing with color images

The straightforward strategy to extend denoising algorithms to color or multivalued images is to apply the algorithm independently to each channel. The use of this simple strategy often introduces color artifacts, easily detected by the eye. Two different strategies are observable in state of the art denoising algorithms.

Depending on the algorithm formulation, a vector-valued version dealing at the same time with all color channels can be proposed. This solution is adopted by averaging filters like neighborhood filters or NL-means. These algorithms compute color differences directly in the vector valued image, thus yielding a unified weight configuration which is applied to each channel.

The alternative option is to convert the usual RGB image to a different color space where the independent denoising of each channel does not create noticeable color artifacts. Most algorithms use the YUV system which separates the geometric and chromatic parts of the image. This change writes as a linear transform by multiplication of the RGB vector by the matrix

$$YUV = \begin{pmatrix} 0.30 & 0.59 & 0.11 \\ -0.15 & -0.29 & 0.44 \\ 0.61 & -0.51 & -0.10 \end{pmatrix}, \quad Y_o U_o V_o = \begin{pmatrix} \frac{1}{3} & \frac{1}{3} & \frac{1}{3} \\ \frac{1}{2} & 0 & -\frac{1}{2} \\ \frac{1}{4} & -\frac{1}{2} & \frac{1}{4} \end{pmatrix}$$

The second color transform to the space $Y_o U_o V_o$ is an orthogonal transform. It has the advantage of maximizing the noise reduction of the geometric component, since this component is an average of the three color. The geometric component is perceptually more important than the chromatic ones, and the presence of less noise permits a better performance of the algorithm in this part. It also permits a higher noise reduction on the chromatic components $U_o$ and $V_o$, due to their observable regularity.

This latter strategy is adopted by transform thresholding filters for which the design of an orthonormal basis coupling the different color channels is not trivial.

### 4. Trying all generic tools on an example

This chapter applies incrementally the previous generic denoising tools to the DCT sliding window to illustrate how these additional tools permit to drastically improve the algorithm performance. We start with the basic DCT “neighborhood filter” as proposed by Yaroslavsky [28]. Its principle is to denoise a patch around each pixel, and to keep only the central denoised pixel.

Figure 1 displays the denoised images obtained by incrementally applying each of the following ingredients:

- Basic DCT thresholding algorithm by the neighborhood filter technique (keeping only the central pixel of the window). Each color channel is treated independently.

- Use of an orthogonal geometric and chromatic decomposition color system $Y_o U_o V_o$; grey parts are better reconstructed and color artifacts are reduced.

- Uniform aggregation; the noise reduction is superior and isolated noise points are removed.
- Adaptive aggregation using the estimator variance; the noise reduction near edges is increased, "halo" effects are removed.

- Additional iteration using “oracle” estimation; residual noise is totally removed and the sharpness of details is increased.

The *PSNR's* obtained by incrementally applying the previous strategies respectively are 26.85, 27.33, 30.65, 30.73, 31.25. This experiment illustrates how the use of these additional tools is crucial to achieve competitive results. This last version of the DCT denoising algorithm, incorporating all the proposed generic tools, will be the one used in the comparison chapter. A complete description of the algorithm can be found in Algorithm 17. The color version of the algorithm applies the denoising independently to each $Y_oU_oV_o$ component. This version is therefore slightly better than the version online in [95], which does not use the oracle step.
**Algorithm 17** DCT denoising algorithm. DCT coefficients lower than $3\sigma$ are canceled in the first step and a Wiener filter is applied in the “oracle” second step. The color DCT denoising algorithm applies the current strategy independently to each $Y_o U_o V_o$ component.

1: **Input:** noisy image $\tilde{u}$, $\sigma$ noise standard deviation.

2: **Optional:** prefiltered image $\hat{u}_1$ for “oracle” estimation.

3: **Output:** output denoised image.

4: Set parameter $\kappa = 8$: size of patches.

5: Set parameter $h = 3\sigma$: threshold parameter.

6: **for** each pixel i **do**

7: Select a square reference patch $\tilde{P}$ around i of size $\kappa \times \kappa$.

8: **if** $\hat{u}_1$ **then**

9: Select a square reference patch $P_1$ around i in $\hat{u}_1$.

10: **end if**

11: Compute the DCT transform of $\tilde{P}$.

12: **if** $\hat{u}_1$ **then**

13: Compute the DCT transform of $P_1$.

14: **end if**

15: **if** $\hat{u}_1$ **then**

16: Modify DCT coefficients of $\tilde{P}$ as

$$ \tilde{P}(i) = \tilde{P}(i) \frac{P_1(i)^2}{P_1(i)^2 + \sigma^2} $$

17: **else**

18: Cancel coefficients of $\tilde{P}$ with magnitude lower than $h$.

19: **end if**

20: Compute the inverse DCT transform obtaining $\tilde{P}$.

21: Compute the aggregation weight $w_{\tilde{P}} = 1/\#\{\text{number of non-zero DCT coefficients}\}$.

22: **end for**

23: **for** each pixel i **do**

24: Aggregation: recover the denoised value at each pixel i by averaging all values at i of all denoised patches $\tilde{Q}$ containing i, weighted by $w_{\tilde{Q}}$.

25: **end for**
<figure><img src="image_71.png" /><figcaption>FIGURE 1. Top: original and noisy images with an additive Gaussian white noise of standard deviation 25. Below, from top to bottom and left to right: crop of denoised images by sliding DCT thresholding filter and incrementally adding use of a $Y_oU_oV_o$ color system, uniform aggregation, variance based aggregation and iteration with the “oracle” given by the first step. The corresponding PSNR are 26.85, 27.33, 30.65, 30.73, 31.25.</figcaption></figure>
<figure><img src="image_72.png" /></figure>
# Detailed analysis of the Non-Local Means and the Non-local Bayes methods

In this chapter, we give a detailed description and analysis of two denoising methods for which reliable faithful implementations are available: the Non-local means and Non-local Bayes methods. The algorithmic descriptions of both methods is provided. We shall extend this non-local Bayes algorithm in our final Chapter 8 on the Noise Clinic.

## 1. Non-local means

The Non-local means (NL-means) algorithm tries to take advantage of the redundancy of most natural images. The redundancy, or self-similarity hypothesis is that for every small patch in a natural image one can find several similar patches in the same image, as illustrated in figures 1 and 2. This similarity is true for patches whose centers stand at a one pixel distance of the center of the reference patch. In that case the self-similarity boils down to a local image regularity assumption. Such a regularity is guaranteed by Shannon-Nyquist's sampling conditions, which require the image to be blurry. In a much more general sense inspired by neighborhood filters, one can define as "neighborhood of a pixel i" any set of pixels j in the image such that a patch around j looks like a patch around i. All pixels in that neighborhood can be used for predicting the value at i, as was first shown in [102] for the synthesis of texture images. This self-similarity hypothesis is a generalized periodicity assumption. The use of self-similarities is actually well-known in information theory from its foundation. In his 1948 Mathematical Theory of Communication, Shannon [25] analyzed the local self-similarity (or redundancy) of natural written language, and gave probably the first stochastic text synthesis algorithm. The Efros-Leung texture synthesis method adapted this algorithm to images, and NL-Means [103] seems to be first adaptation of the same idea to denoising<fnref n="1" />

NL-means denoises a square reference patch $\tilde{P}$ around i of dimension $κ × κ$ by replacing it by an average of all similar patches $\tilde{Q}$ in a square neighborhood of i of size $λ × λ$. To do this, a normalized Euclidean distance between $\tilde{P}$ and $\tilde{Q}$, $d(\tilde{P}, \tilde{Q}) = \frac{1}{κ^2}||\tilde{P} - \tilde{Q}||^2$ is computed for all patches $\tilde{Q}$ is the search neighborhood. Then the weighted average is

<footnote n="1">¹Nevertheless, some researchers have pointed out to us the report [104] as giving an early intuition that intuition could use signal redundancy. This very short paper describes an experiment in a few sentences. It suggests that region redundancy on both sides of an edge can be detected, and used for image denoising. Nevertheless, no algorithm is specified in this paper.</footnote>
<figure><img src="image_73.png" /><figcaption>FIGURE 1. q1 (a) and q2 (c) have a large weight because their similarity windows are similar to that of p (b). On the other side, for q3 (d) weight w(p, q3) is much smaller because the intensity grey values in the similarity windows are very different.</figcaption></figure>

$$ \hat{P} = \frac{\sum_{\tilde{Q}} \tilde{Q} e^{-\frac{d(\tilde{P}, \tilde{Q})^2}{h^2}}}{\sum_{\tilde{Q}} e^{-\frac{d(\tilde{P}, \tilde{Q})^2}{h^2}}}. $$

The thing that helps NL-means over the neighborhood filters is the concentration of the noise law, as the number of pixels increases. Because the distances are computed on many patch samples instead of only one pixel, the fluctuations of the quadratic distance due to the noise are reduced.

**Related attempts:** [105] proposed a “universal denoiser” for digital images. The authors prove that this denoiser is universal in the sense “of asymptotically achieving, without access to any information on the statistics of the clean signal, the same performance as the best denoiser that does have access to this information”. In [106] the authors presented an implementation valid for binary images with an impulse noise, with excellent results. Awate and Whitaker [107] also proposed a method whose principles stand close to NL-means, since the method involves comparison between patches to estimate a restored value. The objective of the algorithm is to denoise the image by decreasing the randomness of the image.

**A consistency theorem for NL-means.** NL-means is intuitively consistent under stationarity conditions, namely if one can find many samples of every image detail. It can be proved [32] that if the image is a fairly general stationary and mixing random process, for every pixel i, NL-means
<figure><img src="image_74.png" /><figcaption>FIGURE 2. On the right-hand side of each pair, one can see the weight distribution used to estimate a centered patch of the left image by NL-means. (a) In flat zones, the weights are uniformly distributed, NL-means acts like a low pass isotropic filter. (b) On straight edges, the weights are distributed in the direction of the edge (like for anisotropic filters). (c) On curved edges, the weights favor pixels belonging to the same contour. (d) In a flat neighborhood, the weights are distributed in a grey-level neighborhood (exactly like for neighborhood filters). In the cases of (e) and (f), the weights are distributed across the more similar configurations, even though they are far away from the observed pixel. This behavior justifies the "non local" appellation.</figcaption></figure>

converges to the conditional expectation of *i* knowing its neighborhood, which is the best Bayesian estimate.

**NL-means as an extension of previous methods.** A Gaussian convolution preserves only flat zones, while contours and fine structure are removed or blurred. Anisotropic filters instead preserve straight edges, but flat zones present many artifacts. One could think of combining these two methods to improve both results. A Gaussian convolution could be applied in flat zones, while an anisotropic filter could be applied on straight edges. Still, other types of filters should be designed to specifically restore corners, or curved edges, or periodic texture. Figure 2 illustrates how NL-means chooses the right weight configuration for each sort of image self-similarity.
**Algorithm 18** NL-means algorithm (parameter values for $κ$, $λ$ are indicative).

1: **Input:** noisy image $\tilde{u}$, $\sigma$ noise standard deviation.

2: **Output:** output denoised image.

3: Set parameter $κ = 3$: size of patches.

4: Set parameter $λ = 31$: size of research zone for which similar patches are searched.

5: Set parameter $h = 0.6 \sigma$: bandwidth filtering parameter.

6: **for** each pixel **i** **do**

7: Select a square reference patch $\tilde{P}$ around **i** of dimension $κ \times κ$.

8: Set $\hat{P} = 0$ and $\hat{C} = 0$.

9: **for** each patch $\tilde{Q}$ in a square neighborhood of **i** of size $λ \times λ$ **do**

10: Compute the normalized Euclidean distance between $\tilde{P}$ and $\tilde{Q}$, $d(\tilde{P}, \tilde{Q}) = \frac{1}{κ^2} \|\tilde{P} - \tilde{Q}\|^2$.

11: Accumulate $\tilde{Q}e^{-\frac{d(\tilde{P}, \tilde{Q})^2}{h^2}}$ to $\hat{P}$ and $e^{-\frac{d(\tilde{P}, \tilde{Q})^2}{h^2}}$ to $\hat{C}$.

12: **end for**

13: Normalize the average patch $\hat{P}$ by dividing it by the sum of weights $\hat{C}$

14: **end for**

15: **for** each pixel **x** **do**

16: Aggregation: recover the denoised value at each pixel **i** by averaging all values at **i** of all denoised patches $\hat{Q}$ containing **i**

17: **end for**
## 2. Non-local Bayesian denoising

It is apparent that Equation (11) given in Chapter 5,

$$ \hat{P}_1 = \bar{P} + [C_{\tilde{P}} - \sigma^2 \mathbf{I}] C_{\tilde{P}}^{-1} (\tilde{P} - \bar{P}), $$

gives by itself a denoising algorithm, provided we can compute the patch expectations and patch covariance matrices. We shall now explain how the Non-local Bayes algorithm proposed in [47, 91] does it. Let $\mathcal{P}(\tilde{P})$ be the set of patches $\tilde{Q}$ similar to the patch $\tilde{P}$, which have obtained with a suitably chosen tolerance threshold, so that we can assume that they represent noisy versions of the patches similar to $P$. Then, by the law of large numbers, we have

$$ (16) \qquad C_{\tilde{P}} \approx \frac{1}{\#\mathcal{P}(\tilde{P}) - 1} \sum_{\tilde{Q} \in \mathcal{P}(\tilde{P})} (\tilde{Q} - \bar{P})(\tilde{Q} - \bar{P})^t, \quad \bar{P} \approx \frac{1}{\#\mathcal{P}(\tilde{P})} \sum_{\tilde{Q} \in \mathcal{P}(\tilde{P})} \tilde{Q}. $$

Nevertheless, the selection of similar patches at the first step is not optimal and can be improved in a second estimation step where the first step estimation is used as oracle. Thus, in a second step, where all patches have been denoised at the first step, all the denoised patches can be used again to obtain an estimation $C_{\hat{P}_1}$ for $C_P$, the covariance of the cluster containing $P$, and a new estimation of $\bar{P}$, the average of patches similar to $\tilde{P}$. Indeed, the patch similarity is better estimated with the denoised patches. Then it follows from (10) and (11) that we can obtain a second better denoised patch,

$$ (17) \qquad \hat{P}_2 = \bar{P}^1 + C_{\hat{P}_1} [C_{\hat{P}_1} + \sigma^2 \mathbf{I}]^{-1} (\tilde{P} - \bar{P}^1) $$

where

$$ (18) \qquad C_{\hat{P}_1} \approx \frac{1}{\#\mathcal{P}(\hat{P}_1) - 1} \sum_{\hat{Q}_1 \in \mathcal{P}(\hat{P}_1)} (\hat{Q}_1 - \bar{P}^1)(\hat{Q}_1 - \bar{P}^1)^t, \quad \bar{P}^1 \approx \frac{1}{\#\mathcal{P}(\bar{P}_1)} \sum_{\bar{Q}_1 \in \mathcal{P}(\bar{P}_1)} \bar{Q}. $$

We write the denoised patches $\bar{P}$ in (16) and $\bar{P}^1$ in (17). Indeed, in (17), the denoised version of $\tilde{P}$ computed as the average of noisy patches $\tilde{Q}$ whose denoised patch is similar to $\hat{P}_1$.

In short, the estimates (11) and (17) appear equivalent, but they are not in practice. $C_{\hat{P}_1}$, obtained after a first denoising step, is a better estimation than $C_{\tilde{P}}$. Furthermore, $\bar{P}^1$ is a more accurate mean than $\bar{P}$. It uses a better evaluation of patch similarities. All above quantities being computable from the noisy image, we obtain the two step algorithm 19.

As pointed out in [47, 91], the Nonlocal Bayes algorithm only is an interpretation (with some generic improvements like the aggregation) of the PCA based algorithm. This paper has a self-explanatory title: “Two-stage image denoising by principal component analysis with local pixel grouping.” It is equivalent to apply a PCA on the patches similar to $\tilde{P}$, followed by a Wiener filter on the coefficients of $\tilde{P}$ on this PCA, or to apply formula (11) with the covariance matrix of the similar patches. Indeed the PCA computes nothing but the eigenvalues of the empirical covariance
**Algorithm 19** Non local Bayes image denoising

1: **Input:** noisy image

2: **Output:** denoised image

3: **for** all patches $\tilde{P}$ of the noisy image **do**

4: Find a set $\mathcal{P}(\tilde{P})$ of patches $\tilde{Q}$ similar to $\tilde{P}$.

5: Compute the expectation $\tilde{\mathcal{P}}$ and covariance matrix $\mathcal{C}_{\tilde{\mathcal{P}}}$ of these patches by

$$ \mathcal{C}_{\tilde{\mathcal{P}}} \approx \frac{1}{\#\mathcal{P}(\tilde{\mathcal{P}}) - 1} \sum_{\tilde{Q} \in \mathcal{P}(\tilde{\mathcal{P}})} (\tilde{Q} - \tilde{\mathcal{P}})(\tilde{Q} - \tilde{\mathcal{P}})^t, \quad \tilde{\mathcal{P}} \approx \frac{1}{\#\mathcal{P}(\tilde{\mathcal{P}})} \sum_{\tilde{Q} \in \mathcal{P}(\tilde{\mathcal{P}})} \tilde{Q}. $$

6: Obtain the first step estimation:

$$ \hat{P}_1 = \bar{\mathcal{P}} + [\mathcal{C}_{\bar{\mathcal{P}}} - \sigma^2 \mathbf{I}] \mathcal{C}_{\bar{\mathcal{P}}}^{-1} (\tilde{P} - \bar{\mathcal{P}}). $$

7: **end for**

8: Obtain the pixel value of the basic estimate image $\hat{u}_1$ as an average of all values of all denoised patches $\hat{Q}_1$ which contain **i**.

9: **for** all patches $\tilde{P}$ of the noisy image **do**

10: Find a new set $\mathcal{P}_1(\tilde{P})$ of noisy patches $\tilde{Q}$ similar to $\tilde{P}$ by comparing their denoised “oracular” versions $Q_1$ to $P_1$.

11: Compute the new expectation $\bar{\mathcal{P}}^{-1}$ and covariance matrix $\mathcal{C}_{\bar{\mathcal{P}}_1}$ of these patches:

$$ \mathcal{C}_{\bar{\mathcal{P}}_1} \approx \frac{1}{\#\mathcal{P}(\bar{\mathcal{P}}_1) - 1} \sum_{\tilde{Q}_1 \in \mathcal{P}(\bar{\mathcal{P}}_1)} (\hat{Q}_1 - \bar{\mathcal{P}}^{-1})(\hat{Q}_1 - \bar{\mathcal{P}}^{-1})^t, \quad \bar{\mathcal{P}}^{-1} \approx \frac{1}{\#\mathcal{P}(\bar{\mathcal{P}}_1)} \sum_{\tilde{Q}_1 \in \mathcal{P}(\bar{\mathcal{P}}_1)} \tilde{Q}_1. $$

12: Obtain the second step patch estimate

$$ \hat{P}_2 = \bar{\mathcal{P}}^{-1} + \mathcal{C}_{\hat{P}_1} [\mathcal{C}_{\hat{P}_1} + \sigma^2 \mathbf{I}]^{-1} (\tilde{P} - \bar{\mathcal{P}}^{-1}). $$

13: **end for**

14: Obtain the pixel value of the denoised image $\hat{u}(\mathbf{i})$ as an average of all values of all denoised patches $\hat{Q}_2$ which contain **i**.

matrix. Thus, the method in [91] gets its Bayesian interpretation. A study on the compared performance of local PCA versus global PCA for TSID is actually proposed in [70].
# The “Noise Clinic”: a universal denoiser

Arguably several thousands papers are dedicated to image denoising. Most papers assume a fixed noise model, mainly white Gaussian or Poissonian (Chapters 1 and 2). This assumption is only valid for raw images. Yet in most images handled by the public and even by scientists, the noise model is imperfectly known or unknown. End users only dispose of the result of a complex image processing chain (Chapter 3) effectuated by uncontrolled hardware and software (and sometimes by chemical means). For such images, it is possible to estimate from a single image a noise model which is simultaneously signal and frequency dependent, as we showed in Chapter 4. In this chapter we propose a multiscale denoising algorithm adapted to this broad noise model, based on the NL-Bayes denoising algorithm (Chapter 7). This leads to a blind denoising algorithm which we demonstrate on real JPEG images and on scans of old photographs for which the formation model is unknown. This algorithm is finally compared to the unique state of the art previous blind denoising method, based on Gaussian Scale Mixtures In The Wavelet Domain [14].

## 1. Introduction

1.1. Motivations. Blind denoising is the conjunction of a thorough noise estimation method followed by the application of an adapted denoising method. To cope with the broad variety of observed imaging noises, the noise model must be far more comprehensive than the usual white Gaussian noise. Our lead example will be JPEG images from digital CCD or CMOS cameras, where the initial signal dependent white Poisson noise has undergone nonlinear transforms, linear filters and a quantization of its DCT coefficients. After such alterations, a signal, frequency and scale dependency is a minimal assumption for the remaining noise. This requires dealing with a noise model depending on hundreds of parameters, in contrast with the usual one-parameter Gaussian white noise and the two-parameter Poisson noise. A flexible denoising method must also be conceived to cope with this signal, scale, and frequency dependent noise model.

To be useful to all image users, who generally have only access to the end result of a complex processing chain, blind denoising must be able to cope with both raw and preprocessed images of all sorts. The archives of the online executions at the IPOL journal of six classic denoising methods, namely DCT denoising [108], TV denoising [109], K-SVD [110], NL-means [111], BM3D [112] and NL-Bayes [84] are replete with such puzzling noisy images. IPOL users are in principle requested to upload noiseless images, to which the noise is added on line to test the performance of each algorithm. Yet, as one can observe in this public archive, the demand for a blind denoiser is so
strong that more than 10000 noisy images have been unduly uploaded. This shows how necessary
“blind” methods are, for diffusing image processing techniques in science and technology.

**1.2. Antecedents.** We found only a few references on blind denoising approaches: Portilla [37], [14], Rabie [113] and Liu, Freeman, Szeliski and Kang [52]. Portilla's method is an adaptation of the famous BLS-GSM algorithm, which models wavelet patches at each scale by a Gaussian scale mixture (GSM), followed by a Bayesian least square (BLS) estimation for wavelet patches. This method is in principle adapted to homogeneous, Gaussian or mesokurtic noise. Yet, according to the author, the GSM model provides an automatic way to separate noise from signal. Indeed, for natural images, a GSM captures for the wavelet coefficients both high kurtosis marginals and a positive covariance between neighbor coefficient amplitudes. These coefficients are not shared by Gaussian or lower kurtosis noise sources. Then, for each wavelet subband a correlated Gaussian model can be used to estimate the noise and a correlated GSM is used for the signal. This algorithm is fully automatic, and will be compared to our results in Chapter 6.3. Our proposed solution shares many features with Portilla's method. Our noise model is nonetheless more general, being signal dependent, and our patch model is local, while the GSM wavelet patch model is global. (A recent local version of BLS-GSM [38] obtains a better performance than BLS-GSM.)

Liu, Freeman, Szeliski and Kang [52] proposed a unified denoising framework for JPEG images with two tasks in view: 1) automatic estimation and 2) removal of colored noise from a single image. These steps are performed by involving a piecewise smooth image model and a segmentation. The authors introduce the so called “noise level functions” (NLF) to estimate the noise level as a function of the image grey level. The obtained noise curve by their algorithm is an estimate of an upper bound of the real NLF, done by fitting a lower envelope to the standard deviations of per-segment image variances. In their denoising procedure, the chrominance of the colored noise is significantly removed by projecting pixel values onto a line fitted to the RGB values in each segment. Then, a Gaussian conditional random field is constructed to obtain the underlying clean image from the noisy input. Unfortunately no code is available for this complex procedure.

The method proposed by Rabie [113] seems less effective and works only for Gaussian noise. Here the blind denoising filter is based on the theory of robust statistics. The denoising part is done by minimizing a stationary cost function. For an adaptive window around the pixel of interest, noise pixels are seen as outlier pixels and rejected according to the Lorentzian robust estimator. The noise is basically estimated over a flat area of the noisy image. “Optimal-size” adaptive window are used to obtain the largest area containing relatively uniform structures around each pixel of interest. The uniformity is based on local signal variance estimate. This method seems less general than Portilla’s method, since it can only deal with a signal-independent Gaussian noise. Observing the results shown in [113], indicates that this method mainly works on images with large homogeneous areas. An entropy-based noise level estimator has been proposed in [114], which may work for any sort of noise. Unfortunately it delivers a noise level but not a noise model. So we could not use it for noise estimation. Our denoising method will be based on a noise signal and frequency noise estimator proposed by Colom et al. [10], relying on a Ponomarenko et al. general
principle [83] to build a noise patch model. This method is proved in the aforementioned reference to estimate accurately the variances of DCT coefficients of noise patches in a JPEG image. We shall see that it can be easily extended to cope with a scale dependency.

## 2. A generalized nonlocal Bayesian algorithm

Most denoising methods in the literature focus on Gaussian white noise, which is a reasonable simplification of the problem, since for example Poisson noise can be transformed into approximately white Gaussian noise by the Anscombe transform [85]. In this chapter we show that one of them, the NL-Bayes method, designed for Gaussian white noise, can be extended to deal with a signal, scale and frequency dependent noise. NL-Bayes only requires the knowledge of a local Gaussian patch model and of a Gaussian noise model. It is therefore possible to extend the noise model to make obtain a denoising method compatible with a scale and signal dependent.

Like other patch based denoising methods, NL-Bayes denoises all noisy square patches extracted from the noisy image $\tilde{u}$ and then obtains the final denoised image $\hat{u}$ by replacing every image pixel value by an average of the denoised values obtained for this pixel in all denoised patches containing it. We shall denote by $\tilde{P}$ a reference patch extracted from the image, and by $\mathcal{P}(\tilde{P})$ a set of patches $\tilde{Q}$ similar to the reference patch $\tilde{P}$. Assuming that $\tilde{Q}$ follows a Gaussian model, a first basic estimation of the denoised patch $P$ can be obtained [78] by

$$ (19) \qquad P^{\text{basic}} = \bar{P} + [C_{\tilde{P}} - C_n] C_{\tilde{P}}^{-1} (\tilde{P} - \bar{P}) $$

where

* $\bar{P}$ is the empirical average of the patches similar to $\tilde{P}$:

$$ (20) \qquad \bar{P} \approx \frac{1}{\#\mathcal{P}(\tilde{P})} \sum_{\tilde{Q} \in \mathcal{P}(\tilde{P})} \tilde{Q}; $$

* $C_n$ is the covariance matrix of the noise;

* $C_{\tilde{P}}$ is the empirical covariance matrix of the patches similar to $\tilde{P}$, which may be obtained by

$$ (21) \qquad C_{\tilde{P}} \approx \frac{1}{\#\mathcal{P}(\tilde{P}) - 1} \sum_{\tilde{Q} \in \mathcal{P}(\tilde{P})} (\tilde{Q} - \bar{P})(\tilde{Q} - \bar{P})^t. $$

For pure Gaussian signal-independent noise with variance $\sigma^2$, we simply have $C_n = \sigma^2 I$. The above estimate would the optimal Bayesian estimate, if $C_{\tilde{P}}$ and $\bar{P}$ were the true covariance matrix and expectation of the patches similar to $\tilde{P}$. In a second step, all the denoised patches obtained after the previous first step estimation can be reused by a classic Wiener argument to obtain a better unbiased estimation $C_{\tilde{P}}^{\text{basic}}$ for the covariance of the 3D group containing $P$. Similarly, a new estimation $\bar{P}^{\text{basic}}$ of the average of patches similar to $P$ can be obtained. This leads to a second Wiener-Bayes estimate

$$ (22) \qquad P^{\text{final}} = \bar{P}^{\text{basic}} + C_{\tilde{P}}^{\text{basic}} [C_{\tilde{P}}^{\text{basic}} - C_n]^{-1} (\tilde{P} - \bar{P}^{\text{basic}}). $$
2.0.0.1. Adaptation to Signal-Dependent Noise. As formulas (19) and (22) show, the above Bayesian principle is compatible with a patch noise model $C_n$ depending on each patch $\tilde{P}$. The above formulas only require a good estimate of the covariance matrix of the noise associated with each group of similar patches. The algorithm computing this matrix is given in Chapter 3. The noise model being signal dependent, for each intensity $i$ in the range intensity [0, 255] of the image a noise covariance matrix $C_{ni}$ will be available. The noise model for each group of patches similar to $\tilde{P}$ will depend on $\tilde{P}$ through their mean $i$. The reference intensity for the current 3D group $\mathcal{P}(\tilde{P})$ must therefore be estimated to apply formulas (19) and (22) with the appropriate noise covariance matrix. This intensity is simply estimated as the average of all pixels contained in $\mathcal{P}(\tilde{P})$.

2.0.0.2. Local Correction of the Covariance Matrix. The denoising performance strongly depends on the noise covariance matrices estimation. If the matrices $\{C_{ni}\}_{i \in [0,255]}$ are not accurate enough, denoising can cause ugly artifacts, particularly in the first step. The noise estimation procedure from the image is always at risk of an overestimation, particularly when the image is small or when it contains a uniform texture which becomes indistinguishable from colored noise. If $C_n$ is overestimated, then (19) risks adding “negative noise” to the image, because of the $-C_n$ term in this equation. Thus, a conservative estimation strategy must be applied on the first Bayesian step to avoid noise overestimation artifacts. This strategy ensures that the noise variances are always smaller than the noisy patch variances. This sanity check based on the diagonal values of both $C_{\tilde{P}}$ and $C_n$ covariance matrices leads to the following more conservative estimate of the diagonal elements of the patch covariance matrix used in (19):

$$ (23) \qquad \forall p \in [0, \kappa^2 - 1], C_{\tilde{P}}(p, p) = \max(C_{\tilde{P}}(p, p), C_n(p, p)). $$

2.0.0.3. Homogeneous Area Detection. The original NL-Bayes algorithm [78] has a statistical test to determine if a 3D group belongs to a homogeneous area, and in this case the estimation of all patches is replaced by the global mean over all pixels contained in the 3D group. This criterion is merely based on the comparison of the empirical standard deviation of all pixels of $\mathcal{P}(\tilde{P})$ with $\sigma^2$.

In our generalization of this algorithm, $\sigma$ doesn't exist since $C_n \neq \sigma^2 I$. So this criterion must be adapted to better take into account $C_n$ in the following way:

* First, compute the difference of the traces of both covariance matrices for each channel $c$,

$$ (24) \qquad \delta_c = \mathrm{Tr}(C_{\tilde{P}}) - \mathrm{Tr}(C_n). $$

* Denote by $\hat{Q}$ a first estimation of $\tilde{Q}$ obtained by (19). Then the basic estimate is $\forall \tilde{Q} \in \mathcal{P}(\tilde{P})$,

$$ (25) \qquad Q^{\text{basic}} = \begin{cases} \bar{\bar{\tilde{P}}} & \text{if } \delta_c < \alpha \mathrm{Tr}(C_n) \\ \hat{Q} & \text{if } \delta_c > \beta \mathrm{Tr}(C_n) \\ t\hat{Q} + (1-t)\bar{\bar{\tilde{P}}} & \text{otherwise.} \end{cases} $$
where

$$t = \frac{\delta_c - \alpha \mathrm{Tr}(\mathbf{C}_n)}{\beta \mathrm{Tr}(\mathbf{C}_n) - \alpha \mathrm{Tr}(\mathbf{C}_n)}$$

and

$$\bar{\bar{P}} = \frac{1}{\#\mathcal{P}(\tilde{P})\kappa^2} \sum_{\tilde{Q} \in \mathcal{P}(\tilde{P})} \sum_{p=1}^{\kappa} \sum_{q=1}^{\kappa} \tilde{Q}(p, q)$$

The thresholds $(\alpha, \beta)$ are chosen equal to $(-\frac{1}{3}, \frac{1}{3})$. This (optional) correction which generally increases the PSNR is only used for the first step of the finest scale of the multiscale algorithm.

### 3. Obtaining the covariance matrix of noise patches

Colom et al., [10], proposed an adaptation of the Ponomarenko et al. [83] method estimating a frequency dependent noise to estimate noise in JPEG images. Given a patch size $\kappa \times \kappa$, the method extracts from the image a set with fixed cardinality of sample blocks with very similar patches in DCT space, which are therefore likely to contain only noise. These noise blocks are transformed by a DCT, and an empirical standard deviation of their DCT coefficients is computed. This gives a noise model that is proved in [10] to be accurately consistent with noise observed in JPEG images. This algorithm computes for every intensity $i$ with a multi-frequency noise estimate given by a $\kappa^2 \times \kappa^2$ matrix

$$ (26) \qquad \mathbf{M}_i := \mathbb{E} (\mathcal{D}N_i (\mathcal{D}N_i)^t) $$

where:

* $\mathcal{D}$ is the $\kappa^2 \times \kappa^2$ matrix of the discrete cosine transform (DCT);

* $N_i$ denotes the $\kappa \times \kappa$ stochastic noise patch model at intensity $i$.

#### 3.1. Are noise covariances negligible in the block DCT space?
The method of the preceding chapter only estimates the variances of the DCT coefficients of noise blocks and not their covariances. The covariance matrices are therefore assumed to be diagonal, which amounts to assume that the DCT decorrelates the noise. A formal argument can be given in favor of this assumption. Assume that the initial image noise was white Gaussian, and that the image has undergone a symmetric, real, periodic linear filter $H$. Then this filter corresponds to applying a diagonal operator to the image in the DCT frequency domain. Thus the noise covariance of the filtered noise remains diagonal in the DCT domain. Yet, this argument is only valid for a global image DCT. Here, because we need a signal dependent noise model, we are estimating it on local DCTs applied to each block. It is therefore no more true that the blocks have undergone a periodic convolution filter. Thus, it cannot be exactly true that after the application of a global linear filter, the noise block DCTs have a diagonal covariance. To check nonetheless the quantitative validity of this assumption, we tested three different filters applied to a white noise:

* $\mathbf{H}_1$ with coefficients $\frac{1}{4} \begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix}$ supported by the pixels $(-\frac{1}{2}, -\frac{1}{2}), (\frac{1}{2}, \frac{1}{2}), (\frac{1}{2}, -\frac{1}{2}), (-\frac{1}{2}, \frac{1}{2})$;
<table><caption>TABLE 1. Statistics of the estimated DCT covariance matrix of noise filtered by $\mathbf{H}_1$.</caption><thead><tr><th>κ</th><th>4</th><th>6</th><th>8</th><th>16</th></tr></thead><tbody><tr><td>mean {|C<sub>i,j</sub>|}<sub>i≠j</sub></td><td>0.83</td><td>0.48</td><td>0.31</td><td>0.10</td></tr><tr><td>mean {|C<sub>i,j</sub>|}<sub>i=j</sub></td><td>24.89</td><td>25.31</td><td>24.95</td><td>24.73</td></tr><tr><td>median {|C<sub>i,j</sub>|}<sub>i≠j</sub></td><td>0.04</td><td>0.03</td><td>0.02</td><td>0.01</td></tr><tr><td>median {|C<sub>i,j</sub>|}<sub>i=j</sub></td><td>19.42</td><td>17.14</td><td>16.28</td><td>14.41</td></tr></tbody></table>

<table><caption>TABLE 2. Statistics of the estimated DCT covariance matrix of noise filtered by $\mathbf{H}_2$.</caption><thead><tr><th>κ</th><th>4</th><th>6</th><th>8</th><th>16</th></tr></thead><tbody><tr><td>mean {|C<sub>i,j</sub>|}<sub>i≠j</sub></td><td>0.48</td><td>0.28</td><td>0.19</td><td>0.06</td></tr><tr><td>mean {|C<sub>i,j</sub>|}<sub>i=j</sub></td><td>14.59</td><td>13.95</td><td>14.45</td><td>14.23</td></tr><tr><td>median {|C<sub>i,j</sub>|}<sub>i≠j</sub></td><td>0.010</td><td>0.008</td><td>0.005</td><td>0.002</td></tr><tr><td>median {|C<sub>i,j</sub>|}<sub>i=j</sub></td><td>6.75</td><td>4.50</td><td>3.77</td><td>2.35</td></tr></tbody></table>

• $\mathbf{H}_2$ the centered filter with coefficients $\frac{1}{16}$
$$ \begin{pmatrix} 1 & 2 & 1 \\ 2 & 4 & 2 \\ 1 & 2 & 1 \end{pmatrix}; $$

• $\mathbf{H}_3$ the centered filter with coefficients $\frac{1}{88}$
$$ \begin{pmatrix} 1 & 2 & 4 & 8 & 4 & 2 & 1 \\ 2 & 4 & 8 & 16 & 8 & 4 & 2 \\ 1 & 2 & 4 & 8 & 4 & 2 & 1 \end{pmatrix}. $$

The noise image $\tilde{u}$ was a 256 × 256 Gaussian white noise with mean 128, and standard deviation $\sigma = 20$. After convolution, we extracted $N$ distinct $\kappa \times \kappa$ patches $\{P_n\}_{n \in N}$ from the image and a 2D normalized DCT was applied on them. Finally, their empirical $\kappa^2 \times \kappa^2$ covariance matrix $\mathbf{C}$ was computed as

$$ \begin{array}{l} \forall (p, q), (i, j) \in [0, \kappa - 1]^2, \\ (27) \qquad \mathbf{C}(p, q, i, j) = \frac{1}{N} \sum_{n=1}^{N} \hat{P}(p, q) \hat{P}(i, j) \\ \qquad \qquad - \frac{1}{N^2} \left( \sum_{n=1}^{N} \hat{P}(p, q) \right) \left( \sum_{n=1}^{N} \hat{P}(i, j) \right) \end{array} $$

These covariances matrices can be visualized by the absolute value of their coefficients $|C_{i,j}|$, normalized in $[0, 1]$ so that the largest coefficient is set equal to 1, and the smallest equal to 0. The following color code is used in the visualization: a coefficient appears in blue if it is near 0; in green if it is near 0.5 and in red if it is near 1. The results for various patch sizes are shown in Figure 1. This illustration and the quantitative tables 1, 2 and 3 confirm that the block DCT noise covariance matrices are nearly diagonal.

So from now on, only variance coefficients will be considered in DCT space.

**3.2. Covariance matrix filtering.** Since the noise covariance matrices can only be estimated for sparse bins in the intensity range, an interpolation must be applied to obtain a noise
<table><caption>TABLE 3. Statistics of the estimated DCT covariance matrix of noise filtered by $H_3$.</caption><thead><tr><th>κ</th><th>4</th><th>6</th><th>8</th><th>16</th></tr></thead><tbody><tr><td>mean {|C<sub>i,j</sub>|<sub>i≠j</sub></td><td>0.22</td><td>0.15</td><td>0.10</td><td>0.04</td></tr><tr><td>mean {|C<sub>i,j</sub>|<sub>i=j</sub></td><td>9.28</td><td>8.94</td><td>9.04</td><td>8.51</td></tr><tr><td>median {|C<sub>i,j</sub>|<sub>i≠j</sub></td><td>0.020</td><td>0.016</td><td>0.010</td><td>0.003</td></tr><tr><td>median {|C<sub>i,j</sub>|<sub>i=j</sub></td><td>3.32</td><td>2.86</td><td>2.44</td><td>1.73</td></tr></tbody></table>

<figure><img src="image_75.png" /><figcaption>FIGURE 1. Visualization of the noise covariance matrices in DCT space after applying filter $H_1$ to illustrate that it is almost diagonal. From left to right, top to bottom patch size $κ = 4, 6, 8$.</figcaption></figure>

covariance matrix of the noise for each given intensity. The covariance matrices must be smoothed before such an interpolation. This can be obtained by a regularization of the covariance matrices in DCT space before applying the inverse DCT to get back a covariance matrix in the image domain. We found that a robust regularization could be performed in the following two steps:

(1) For each frequency independently, perform a linear interpolation between the bin values to obtain a noise curve for this frequency, giving the variance as a function of the signal $i$. Smooth this curve by applying a sliding average;

(2) For every bin, replace each matrix coefficient by the median of its four neighbors and itself.

Since the filtering is channel independent, the pseudo-code only describes the filtering for one channel.

3.2.0.1. Getting Back to the Space Domain. For a given intensity $i$, the covariance matrix of the noise is by definition

$$ \mathrm{Cov}(N_i) = \mathbb{E}(N_i N_i^t) $$

which leads to

$$ \begin{align} \mathcal{D}\mathrm{Cov}(N_i)\mathcal{D}^t &= \mathcal{D}\mathbb{E}(N_i N_i^t)\mathcal{D}^t \nonumber \\ &= \mathbb{E}(\mathcal{D}N_i N_i^t\mathcal{D}^t) \nonumber \\ &= \mathbb{E}(\mathcal{D}N_i(\mathcal{D}N_i)^t) \nonumber \\ &= \mathbf{M}_i \tag{28} \end{align} $$

thanks to equation (26). Since $\mathcal{D}^{-1} = \mathcal{D}^t$, then from equation (28) we get

$$ \mathrm{Cov}(N_i) = \mathcal{D}^t \mathbf{M}_i \mathcal{D}. \tag{29} $$
4. The multiscale algorithm

**4.1. Why a multiscale algorithm?** Classic denoising algorithms such as BM3D (Dabov et al. [36]), NL-means (Buades et al. [86]), K-SVD (Mairal et al. [115], [116]), Wiener filters applied on DCT (Yaroslavsky et al. [117], [118]) or on wavelet transform (Donoho et al. [119]) and the total variation minimization (Rudin et al. [120]) achieve good results for moderate noise ($\sigma \le 20$). Yet for larger noise artifacts inherent to each method (and different for each method) start appearing. In particular all keep an often disturbing low frequency noise. A natural idea to deal with low frequency noise is to involve a coarse to fine multiscale procedure, which promises three improvements:

(1) in the patch-based methods, it favors a better patch comparison, because the patch low frequencies are denoised *before* grouping them by similarity for denoising their higher frequencies;

(2) at coarse scales the noise decreases by zoom out, and state-of-the-art algorithms work better;

(3) subsampling the image before denoising amounts to enlarge the size of the neighborhood on which the denoising is performed, thus permitting to grab and remove low frequency noise on larger regions.

A still stronger argument in favor of a multiscale procedure is that in most images submitted by users, the main bulk of the noise is contained in the low frequencies. This is explainable by several factors. In accurately scanned old photographs, the chemical noise is over-sampled and its grain has low frequency components. In JPEG images, compression has strongly attenuated high frequency noise components, but the low frequency components after the third octave are intact.

<figure><img src="image_76.png" /><figcaption>FIGURE 2. A multiscale process is required to remove the low frequency noise. This is particularly apparent in the flat image regions. From left to right: Noisy image ($\sigma = 30$), result of the “Classic NL-Bayes”, result of the multiscale (three scales) NL-Bayes.</figcaption></figure>

To define a coarse to fine multiscale structure, we proceed by a classic oversampled wavelet denoising strategy [121]. The image is convolved by a Haar “mother wavelet”, which is nothing
but a box-filter **F** where each lower scale pixel is the mean of four samples in the higher scale. This cumulates the advantage of dividing the noise standard deviation by two and of maintaining the independence of the samples after down-sampling. By this process a white noise remains white after subsampling. A classic objection to this wavelet method is that the sub-sampled image is aliased and cannot be up-sampled after denoising. The classic wavelet method avoids this obstacle by denoising simultaneously the three wavelet components obtained by convolving the image with the three Haar wavelets, before reconstructing the finer scale. Yet when dealing with patch based methods, it is better to keep all frequency components together to perform a better nonlocal patch comparison. For this reason the proposed multiscale algorithm keeps and processes four channels that are partly redundant. The four channels are obtained by moving the sub-sampling grid by respectively (0, 0), (1, 0), (0, 1), (1, 1). In that way there is enough information for up-sampling after denoising the denoised images at the lower scale.

The above method is multiscale but does not take advantage of the sub-sampling in the lower scales to increase the algorithm speed. A normal multiscale algorithm is only $1 + \frac{1}{4} + \frac{1}{16} + \cdots = \frac{4}{3}$ more complex than the single scale algorithm. Instead a multiscale algorithm keeping all sub-images when sub-sampling will be twice to five times slower, depending on the number of scales involved, (by default two). Yet, the redundancy of this denoising at lower scales notably increases the restoration quality. This is particularly important, as any denoising error on a down-sampled image is amplified by a four-factor after upsampling.

**4.2. The mean sub-sampling method.** We shall denote by $s$ the current dyadic scale of the multiscale algorithm. For the particular case of white noise, the aim of the sub-sampling is to obtain from $\tilde{u}_s$ an image $\tilde{u}_{s+1}$ where the standard deviation of the noise has been divided by two compared to the noise contained in $\tilde{u}_s$. To get this result, one can use a filter $f(i, j)$ satisfying

$$ \sum_{i,j} f(i,j) = 1 \quad \text{and} \quad \sum_{i,j} f(i,j)^2 = \frac{1}{4}. $$

The simplest filter coping with these conditions is the average filter **F**, defined by

$$ \mathbf{F}(i, j) = \begin{cases} \frac{1}{4} & \text{if } (i, j) \in [(0, 0), (0, 1), (1, 0), (1, 1)], \\ 0 & \text{otherwise.} \end{cases} $$

which averages each group of four adjacent neighboring pixels. There are four different filter+sub-sample results, as shown in figure 3. Moreover if the image $\tilde{u}_s$ is well-sampled, so is $\tilde{u}_s * \mathbf{F}$. Thus, the difference image is not aliased. Since all sub-sampled images are available, the noise estimation can work with the same amount of samples at every scale, which favors a good precision on the noise estimation at lower scales. All sub-sampled images must also be denoised. To avoid handling them separately, we introduce here a new procedure to process them jointly in a single image, while avoiding creating artificial borders. The four sub-sampled images are regrouped in one mosaic image, as shown in figure 4. The boundaries of the sub-images are in that way better denoised, because they are included in a smooth larger image.
<figure><img src="image_77.png" /><figcaption>FIGURE 3. Four different ways to average red neighbors of the yellow reference pixel.</figcaption></figure>

<figure><img src="image_78.png" /><figcaption>FIGURE 4. Left: mosaic of the scale 1 sub-images. Right: mosaic of the scale 2 sub-images, The input image has scale 0.</figcaption></figure>

**4.3. The mean up-sampling method.** The aim of the up-sampling is to go back to the upper scale, after denoising the four sub-images obtained by sub-sampling as seen in Chapter 4.2. The four sub-images $\tilde{u}_1, \tilde{u}_2, \tilde{u}_3$ and $\tilde{u}_4$ have their pixel center (resp. in red, purple, green and blue in figure 5) located at the center of four pixels of $\tilde{u}$ (in black in figure 5). Thus they are shifted by $\pm \frac{1}{2}$ in both coordinate directions. The reconstruction of the pixels of $\tilde{u}$ (see the example of the pixel in yellow in figure 5) will be done by averaging their four neighbors, each one belonging to each sub-image.

**4.4. Noise estimation.** If the input noisy image had pure Gaussian noise, then after each sub-sampling the noise should be divided by two and remain white. For raw images it is the case, since (almost) no alteration nor transformations are applied to the original noisy pixels. Then the noise is a Poisson random process, which can be approximated by a signal-dependent Gaussian noise.

However, the proposed algorithm must deal with all kinds of noisy images. A large majority of them are JPEG images where JPEG has quantized DCT coefficients, making the energy decrease as the frequency increases. In such images the noise increases at lower scales, as illustrated in Figure 6, which are the noise curves of the image shown in Figure 16. This figure displays average
<figure><img src="image_79.png" /><figcaption>FIGURE 5. Position of the center of pixels in the original image $\tilde{u}$ in black, in the four sub-images $\tilde{u}_1$ in red, $\tilde{u}_2$ in purple, $\tilde{u}_3$ in green and $\tilde{u}_4$ in blue. The yellow pixel will be reconstructed by averaging the top left red pixel, the top right purple pixel, the bottom left green pixel and the bottom right blue pixel of its four pixel neighborhood.</figcaption></figure>

noise curves for high and low frequencies respectively, in the three scales noise estimation from a JPEG image. The low-frequency noise is not altered by JPEG and becomes a high-frequency noise after three<fnref n="1" /> subsampling operations.

In our redundant noise estimation, the noise covariance matrices are estimated at each dyadic scale. Chapter 4.2 explains how the noise estimation is applied on the mosaic image composed of all sub-images. Then for every scale the same number of samples is available, which allows the noise estimation to retain a decent accuracy even at coarse scales. At each given scale, all sub-images of the mosaic are denoised with the same set of noise covariance matrices.

The whole coarse to fine multiscale procedure is summarized in Algorithm 20. During the sub-sampling the four sub-images are kept and assembled in a mosaic to be denoised together. It follows that for each scale, the mosaic keeps the original image size. Thus the complexity for the whole algorithm is approximately equal to $N$ times the complexity of the one scale algorithm. In the sequel we shall call our proposed algorithm the “Noise Clinic” as it combines a diagnose of the image illness with an immediate cure.

## 5. Validation

Blind denoising is designed mainly for images where the image history is unknown and no ground truth available. But we can test the denoising performance of the Noise Clinic after simulating a whole image processing chain on a Poisson noisy image for which the ground truth is available. One of the worst possible noise distortion is provided by the image processing chain applied in the camera hardware and generally ending with JPEG compression. This chain includes nonlinear corrections on the raw image, followed by some denoising, demosaicing, gamma-correction, white balance and JPEG compression, namely the quantization of local block DCT coefficients. To see to which extent the method works, we started with perturbations consistent with our noise model

<footnote n="1">¹Since JPEG transform is based on the 8 × 8 DCT transform, after three subsamplings the 8 × 8 pixels patches become a single pixel. Thus, at the third scale the noise is only high-frequency and uncorrelated.</footnote>
<figure><img src="image_80.png" /><figcaption>FIGURE 6. Average noise curves for a typical JPEG-encoded image (shown in Figure 16). From left to right: low frequencies, high frequencies. From top to bottom: scale 2, scale 1, scale 0. Instead of being divided by two at each scale (as it should happen with white noise), the noise grows in lower scales, where JPEG has not removed it.</figcaption></figure>

and then simulated a typical camera image processing chain ending with JPEG compression. We first obtained a noise-free raw image $u_{\text{raw}}$ by subsampling a high quality outdoor image. Then a Poisson noisy $\tilde{u}_{\text{raw}}$ was simulated from it. Four validation experiments were performed.

First, we computed a reference denoised version of the image:

* the Noise Clinic was directly applied on $\tilde{u}_{\text{raw}}$ to get $\hat{u}_{\text{raw}}$;
**Algorithm 20 Noise Clinic**

1: **Input**: Noisy image $\tilde{u}_0$

2: **Input**: Number of scales $N$

3: **Output**: Denoised image $\hat{u}_0$

4: **Part 1**: Builds the image scale pyramid and records the difference images

5: **for** each scale $s = 1$ to $N - 1$ **do**

6: Let $\{\tilde{u}_{s-1}^k\}_{k \in [1,4^{s-1}]}$ be the set of noisy subsampled images obtained at the previous scale.
(For scale $s = 1$, it is $\tilde{u}_0$);

7: **for** $k = 1$ to $4^{s-1}$ **do**

8: Downsample $\tilde{u}_{s-1}^k$ into 4 sub-images: $\{\tilde{u}_s^{4(k-1)+i}\}_{i \in [1,4]}$ as described in Chapter 4.2;

9: Save difference images for this scale:

$$ \tilde{d}_{s-1}^k = \tilde{u}_{s-1}^k - \mathbf{U} \left( \{\tilde{u}_s^{4(k-1)+i}\}_{i \in [1,4]} \right)^k . $$

10: **end for**

11: **if** $s = N - 1$ **then**

12: Set $\{\tilde{v}_{N-1}^k\}_k = \{\tilde{u}_{N-1}^k\}_k$

13: Build the noisy mosaic image $\tilde{m}_{N-1}$ from the set of sub-images $\{\tilde{v}_{N-1}^k\}_{k \in [1,4^{N-1}]}$.

14: **end if**

15: **end for**

16: **Part 2**: Estimates noise and denoises bottom-up in the pyramid

17: **for** $s = N - 1$ to 0 **do**

18: Estimate the noise covariance matrices on the mosaic $\tilde{m}_s$ as explained in Chapter 3;

19: Denoise the noisy mosaic image $\tilde{m}_s$ with the NL-Bayes extension of Chapter 2, using the noise covariance matrices $\{\mathcal{D}^t\mathbf{M}_i\mathcal{D}\}_i$ to obtain $\hat{m}_s$;

20: **if** $s > 0$ **then**

21: **for** $k = 1$ to $4^{s-1}$ **do**

22: Up-sample $\{\hat{u}_{s,4(k-1)+i}\}_{i \in [1,4]}$ and add the saved details $\tilde{d}_{s-1}^k$ to get $\tilde{v}_{s-1}^k$ as described in Chapter 4.3

23: **end for**

24: Construct the mosaic image of the next scale $\tilde{m}_{s-1}$ from the set of sub-images $\{\tilde{v}_{s-1}^k\}$.

25: **else**

26: $\hat{u}_0 = \hat{u}_0^1$

27: **end if**

28: **end for**

* a white balance and a gamma correction were applied on $u_{\text{raw}}$, $\tilde{u}_{\text{raw}}$ and $\hat{u}_{\text{raw}}$ to get $u_{\text{rgb}}$, $\tilde{u}_{\text{rgb}}$ and $\hat{u}_{\text{rgb}}$.

Those images will be used as reference, to see how other parts of the image processing chain (such as the demosaicing and the JPEG compression) impact the result of the denoising. Table 4 shows RMSEs between the noisy and denoised images and the reference one. One can also remark that
<table><caption>TABLE 4. RMSE between noisy/denoised images and corresponding reference image ($u_{raw}$) when two and three scales are used. $v_{rgb}$ denotes the result of Blind BLS-GSM for this experiment.</caption><thead><tr><th>ũ<sub>rgb</sub></th><th>ũ<sub>rgb</sub>|<sup>s2</sup></th><th>ũ<sub>rgb</sub>|<sup>s3</sup></th><th>ṽ<sub>rgb</sub></th></tr></thead><tbody><tr><td>8.62</td><td>3.63</td><td>3.65</td><td>6.46</td></tr></tbody></table>

<table><caption>TABLE 5. RMSE between noisy/denoised images and corresponding reference image ($u_d$) when two and three scales are used for the demosaicking experiment. $v_d$ denotes the result of Blind BLS-GSM for this experiment.</caption><thead><tr><th>ũ<sub>d</sub></th><th>ũ<sub>d</sub>|<sup>s2</sup></th><th>ũ<sub>d</sub>|<sup>s3</sup></th><th>ṽ<sub>d</sub></th></tr></thead><tbody><tr><td>8.64</td><td>4.84</td><td>4.84</td><td>6.43</td></tr></tbody></table>

the best result of the denoising (both in terms of RMSEs and visual aspects) is obtained when the Noise Clinic is applied directly before any transformation.

Second, a demosaicing algorithm was added to the image processing chain before calling the denoising part:

* extract the mosaic<fnref n="2" /> of the noise-free image: $u_m = \text{Mosaic}(u_{raw})$;

* do the same for the noisy image: $\tilde{u}_m = \text{Mosaic}(\tilde{u}_{raw})$;

* apply a classic demosaicing method<fnref n="3" /> on both images, followed by a white balance and a gamma correction to get $u_d$ and $\tilde{u}_d$;

* finally apply the Noise Clinic on $\tilde{u}_d$ to get $\hat{u}_d$.

Table 5 shows RMSEs for this experiment. One may notice that after a demosaicing the noise is no more white, and some structures appear in the noise. These structures are preserved and sometimes enhanced by the denoising algorithm, since it is seen as structure and not as noise. This explains why RMSEs are less favorable than when the denoising is directly applied on the raw images.

Third, a complete image processing chain was simulated to obtain a final JPEG compressed image:

* apply a JPEG compression of quality 92 over both $u_d$ and $\tilde{u}_d$ to get $u_{jpeg}$ and $\tilde{u}_{jpeg}$;

* apply the Noise Clinic to get $\hat{u}_{jpeg}$.

Table 6 shows RMSEs for this experiment. Of course, as JPEG compression creates more artifacts and structured noise, results are worse than with the first two experiments. This only means that the denoising should be applied as soon as possible in the whole image processing chain. However,

<footnote n="2">²The mosaic image is obtained by keeping only the Bayer (R Gr Gb B) over a group of four pixels instead of all RGB values.</footnote>

<footnote n="3">³The demosaicing algorithm used in this experiment was Self-similarity Driven Demosaicking algorithm [122], available on IPOL.</footnote>
<table><caption>TABLE 6. RMSE between noisy/denoised images and corresponding reference image ($u_{jpeg}$) for the JPEG experiment, with compression quality of 92. $\hat{v}_{jpeg}$ denotes the result of Blind BLS-GSM for this experiment.</caption><thead><tr><th>tilde u_jpeg</th><th>hat u_jpeg</th><th>s2</th><th>s3</th><th>hat v_jpeg</th></tr></thead><tbody><tr><td>8.70</td><td>5.34</td><td>5.53</td><td>6.30</td><td></td></tr></tbody></table>

<table><caption>TABLE 7. RMSE between noisy/denoised images and corresponding reference image ($u_f$) for the filtered experiment. $\hat{v}_{f1}$ and $\hat{v}_{f2}$ denote results of Blind BLS-GSM for this experiment.</caption><thead><tr><th>tilde u_f</th><th>hat u_f1</th><th>s2</th><th>hat u_f1</th><th>s3</th><th>hat u_f2</th><th>s2</th><th>hat u_f2</th><th>s3</th><th>hat v_f1</th><th>hat v_f2</th></tr></thead><tbody><tr><td>1.58</td><td>0.75</td><td>0.82</td><td>0.75</td><td>0.75</td><td>1.24</td><td>1.28</td><td>1.24</td><td>1.28</td><td>1.24</td><td>1.28</td></tr></tbody></table>

results are not very far from the ideal case, which confirms the interest and the strength of the Noise Clinic.

Fourth, the filter $H_2$ seen in Chapter 3.1 was used to get:

* a reference filtered image: $u_f = H_2 * u_{raw}$;

* a noisy filtered image: $\tilde{u}_f = H_2 * \tilde{u}_{raw}$;

* the result of the Noise Clinic of the noisy filtered image: $\hat{u}_{f1} = NC(H_2 * \tilde{u}_{raw})$;

* the filtered result of the Noise Clinic of the noisy image: $\hat{u}_{f2} = H_2 * NC(\tilde{u}_{raw})$.

Table 7 shows RMSEs associated to this experiment. Of course after this filtering, there only remains low frequency noise, which explains why RMSEs values are better than in the ideal case. However, the Noise Clinic is still able to give good results.

Figure 7 (resp. 8 and 9) shows results associated of the raw experiment (resp. demosaicking and JPEG).

Figure 10 (resp. 11 and 12) shows a comparison between the Noise Clinic and Blind BLS-GSM for the raw experiment (resp. demosaicking and JPEG).

## 6. Results

**6.1. Detailed results.** In this chapter we applied the blind denoising to real noisy images for which no noise model was available. To illustrate the algorithm structure and its action at each scale, we present for each experiment the noisy input image and for each scale:

* the noisy image where noise has already been removed at coarser scales;

* the denoised image at this scale;

* the difference image = noisy - denoised at this scale;

* the average noise curve over high frequencies;

* the average noise curve over low frequencies.
<figure><img src="image_81.png" /><figcaption>FIGURE 7. Visual results of the reference (first) experiment. From top to bottom, and left to right: full noise-free image, crop of the noise-free image $u_{\text{rgb}}$, crop of the noisy image $\tilde{u}_{\text{rgb}}$, crop of the result of the Noise Clinic using two scales $\hat{u}_{\text{rgb}}|^{s^2}$ and crop of the result of the Noise Clinic using three scales $\hat{u}_{\text{rgb}}|^{s^3}$.</figcaption></figure>
<figure><img src="image_82.png" /><figcaption>FIGURE 8. Visual results of the demosaicking (second) experiment. From top to bottom, and left to right: crop of the noise-free image $u_d$, crop of the noisy image $\tilde{u}_d$, crop of the result of the Noise Clinic using two scales $\hat{u}_d$$^{|s^2}$ and crop of the result of the Noise Clinic using three scales $\hat{u}_d$$^{|s^3}$.</figcaption></figure>

For each scale larger than 1, the subsampled images are up-sampled to keep the original image size. Similarly, the noisy image shown at each scale is the sum of the upsampled version of the denoised sub-images of the previous scale and of the still noisy difference image kept in reserve. In other terms this image contains the remaining noise at the current scale; the noise at coarser scales has in principle already been removed. Visual results are shown in Figure 13. The corresponding noise curves are presented in Figure 14. The experiments made on JPEG photographs from unknown sources are obviously noisy but, as the noise curves illustrate, the noise is not white and is signal dependent. This is easily detected by the fact that the noise curves are not flat and that they are not divided by two from a scale to the next, as they should if the noise were white. A typical fact of JPEG images is that the noise increases at the lower scales. This confirms the necessity of a multiscale algorithm.
<figure><img src="image_83.png" /><figcaption>FIGURE 9. Visual results of the JPEG (third) experiment. From top to bottom, and left to right: crop of the noise-free image $u_{jpeg}$, crop of the noisy image $\hat{u}_{jpeg}$, crop of the result of the Noise Clinic using two scales $\hat{u}_{jpeg}|^{s^2}$ and crop of the result of the Noise Clinic using three scales $\hat{u}_{jpeg}|^{s^3}$.</figcaption></figure>

**6.2. Influence of the number of scales.** Theoretically any number of scales could be used. Indeed at a very coarse scale the noise should be almost null and estimated as such, so that no denoising eventually would occur at very coarse scales. In practice however, some structure of the image may be confused with noise in the noise estimation step. Indeed the noise estimation method is tight on very large images on which pure noise samples in large numbers can be found [83]. After several subsamplings, the image becomes too small, and the risk of confusing texture with noise increases. In consequence applying a blind denoising on a small image is increasingly at risk of removing detail when the scale increases. Thus, it is almost always better to use a minimal number of scales, in most cases not more than two. However, we found that for some images with large low frequency noise it is sometimes better to use up to five scales. From that point of view our “blind denoising” is not fully blind and requires an user evaluation of the number of scales involved.
<figure><img src="image_84.png" /><figcaption>FIGURE 10. Visual comparison of the reference (first) experiment. From left to right: crop of the result of the Noise Clinic using three scales $\hat{u}_{rgb}|^{s^3}$ and crop of the result of the Blind BLS-GSM algorithm $\hat{v}_{rgb}$.</figcaption></figure>

<figure><img src="image_85.png" /><figcaption>FIGURE 11. Visual comparison of the demosaicking (second) experiment. From left to right: crop of the result of the Noise Clinic using three scales $\hat{u}_d|^{s^3}$ and crop of the result of the Blind BLS-GSM algorithm $\hat{v}_d$.</figcaption></figure>

Nevertheless our default value is two, and works on a large majority of the images. Illustrations of the use of the “right” number of scales are presented in Figure 15.

For the “Palace” image in Figure 15, five scales are needed to obtain a noise-free result because of the huge low-frequency noise. In the difference image using five scales one can see that some image structure has been included in the noise. Yet, this low frequency loss is harmless, being undetectable in the resulting denoised image.

6.2.0.1. Result on typical low-light JPEG image. The amount of noise is directly related to the amount of light during the acquisition. Images as shown in Figure 16, taken in a bar with low light conditions are typically very difficult to denoise, even if we had directly access to the RAW
<figure><img src="image_86.png" /><figcaption>FIGURE 12. Visual comparison of the JPEG (third) experiment. From left to right: crop of the result of the Noise Clinic using three scales $\hat{u}_{jpeg}^{s3}$ and crop of the result of the Blind BLS-GSM algorithm $\hat{v}_{jpeg}$.</figcaption></figure>

image, due to the huge amount of noise. One can observe big colored spots caused by demosaicing. JPEG compression ends up creating structured noise. The big colored spots are well attenuated by blind denoising, but the structure created by JPEG is partly left. This is easily explained. These artifacts present sharp recurrent structures which are necessarily confused with signal in an algorithm based on image self-similarity.

6.2.0.2. Results on Old Photographs. Scanned old photographs form a vast image corpus for which the noise model can't be anticipated. The noise is chemical, generally with big grain and further altered by the scanning and JPEG encoding. Figures 17 and 18 show results obtained by the Noise Clinic over this kind of noisy images.

**6.3. Comparison to one of the very few available blind denoising algorithms.** We end this experimental chapter with a comparison of the Noise Clinic with blind BLS-GSM introduced in [14] and [37], a state-of-the-art blind denoising algorithm. The comparison was performed on several images with various noise models. BLS-GSM also is a multiscale algorithm modeling wavelet coefficient patches at each scale and making a global sophisticated Bayesian estimation of them as a Gaussian mixture. NL-Bayes instead has a simpler, but local patch Gaussian model. The global patch model in BLS-GSM has to be more complex to cope with the global patch variability.

In Figure 19 noisy images present strongly structured periodic noise, which is remarkably removed by the blind BLS-GSM algorithm, whereas our blind denoising keeps it and even reinforces it. However one can argue that this structured noise may be seen as a repetitive texture belonging to the image and therefore must be treated as detail and not as noise.

In Figure 20 the noise is more “normal” and closer to what can be expected from a natural image, and our blind denoising performs better. Blind BLS-GSM manages to remove some noise, but a slightly structured noise still remains, appearing in horizontal strips.
<figure><img src="image_87.png" /><figcaption>FIGURE 13. Illustration of blind denoising of a JPEG image, the “Frog” image. It is advised to zoom in the high quality .pdf to see detail.</figcaption></figure>

Figures 21, 22 and 23 show comparisons for low-light JPEG image and old Photographs presented in Chapters 6.2.0.1 and 6.2.0.2

7. Discussion

Blind denoising can be performed with minimal assumptions on the nature of the noise. We observed good results on almost any natural image, even if it had been modified by destructive applications such as JPEG compression or chemical processes. Particularly in old photographs,
<figure><img src="image_88.png" /><figcaption>FIGURE 14. Noise estimation of the “Frog” image: The noise in this image is clearly colored: it increases with descending octaves instead of being divided by two, as it should if it were white.</figcaption></figure>

noise can acquire a thick grain which is only efficiently denoised at low scales. This method does not apply to impulse or multiplicative noise and should be extended to such alterations. Also our local noise estimation procedure did not detect the strength of the fully structured noise present
<figure><img src="image_89.png" /><figcaption>FIGURE 15. Blind denoising when varying the number of scales on “Palace”.</figcaption></figure>
<figure><img src="image_90.png" /><figcaption>FIGURE 16. Blind denoising on “Bar”, using three scales. From left to right, top to bottom : input noisy image, crop of the noisy image, crop of the output denoised image, crop of the difference image.</figcaption></figure>

in the third infrared image of Figure 19. The case of a globally frequency dependent noise is of course better treated by Portilla’s method which assumes a global noise model.

We wrote that the proposed method was “signal, scale and frequency” dependent. In fact as indicated by the preceding caveat, the method estimates and processes noise frequencies in the DCT of small blocks. So these frequency coefficient are far less precise than global image frequencies. Furthermore they are scale dependent, since we applied a dyadic subsampling procedure. Since at each dyadic scale, frequencies are estimated for blocks with at least 4 × 4 size, it follows that these scale dependent frequencies overlap. This leads to a redundant denoising since left-over noise at a coarse scale can be estimated again, and removed again at the overlapping finer dyadic scale. This redundancy of estimators is particularly necessary for such a complex noise model. The fact that JPEG images can be denoised in that way was far from granted. Indeed, it is impossible to really model noise in JPEG images, which are the result of a chain of nonlinear operators. It can be argued that our noise signal, frequency and scale dependent noise estimation is not yet general enough to cope with such alterations. This objection is definitely valid for block artifacts apparent in strong JPEG compression. Thus, strongly compressed images where blocking effects dominate remain beyond our scope.
<figure><img src="image_91.png" /><figcaption>FIGURE 17. Blind denoising on “Marilyn”, using two scales. From left to right, top to bottom : input noisy image, crop of the noisy image, crop of the output denoised image, crop of the difference image.</figcaption></figure>
<figure><img src="image_92.png" /><figcaption>FIGURE 18. Blind denoising on “Solvay conference, 1927”, using three scales. From left to right, top to bottom : input noisy image, crop of the noisy image, crop of the output denoised image, crop of the difference image.</figcaption></figure>
<figure><img src="image_93.png" /><figcaption>FIGURE 19. Results of our blind denoising and of Blind BLS-GSM on several images from [14]. From left to right: Noisy image, result of the Noise Clinic, result of the Blind BLS-GSM algorithm. It is advised to zoom in by a 300% factor the digital document to examine details.</figcaption></figure>
<figure><img src="image_94.png" /><figcaption>FIGURE 20. Comparing our blind denoising with Blind BLS-GSM on several images. It is advised to zoom in by a 400% factor the digital document to examine details. From left to right: Noisy image, result of the Noise Clinic, result of the Blind BLS-GSM algorithm.</figcaption></figure>
<figure><img src="image_95.png" /><figcaption>FIGURE 21. Blind denoising on “Bar”. From left to right: crop of the result of the Noise Clinic by using three scales and crop of the result of the Blind BLS-GSM algorithm.</figcaption></figure>

<figure><img src="image_96.png" /><figcaption>FIGURE 22. Blind denoising on Marilyn”. From left to right: crop of the result of the Noise Clinic by using two scales and crop of the result of the Blind BLS-GSM algorithm.</figcaption></figure>
<figure><img src="image_97.png" /><figcaption>FIGURE 23. Blind denoising on “Solvay conference, 1927”. From left to right: crop of the result of the Noise Clinic by using three scales and crop of the result of the Blind BLS-GSM algorithm.</figcaption></figure>
Part 3

REPRODUCIBLE RESEARCH
CONTRIBUTIONS
This part contains a fully detailed and reproducible account of the three main basic noise estimation methods that we extended in our work to reach more and more general noise estimation methods, culminating in the SFD noise model estimation. We present in order a general method to make any homoscedastic noise estimation method into a signal-dependent noise estimation (Chapter 9).

The three methods detailed here are:

* The Ponomarenko et al. method (Chapter 10), which estimates the standard deviation of the noise using the high-frequency coefficients on a low quantile of the blocks in the image whose energy measured in the low-frequency coefficients is minimal.

* The Percentile method (Chapter 11), which estimates the standard deviation of the noise from the blocks of a high-passed version of the noisy image, using only the blocks whose standard deviation is under a small quantile.

* The PCA method (Chapter 12), which estimates the variance of the noise from the eigenvalues of the covariance matrix of the blocks of the noisy image.

In Chapter 13 these three methods are evaluated in different scenarios.
How to adapt homoscedastic noise estimators to
signal-dependent noise

Most block-based homoscedastic noise estimators can be easily adapted to deal with signal-
dependent noise. Even if the signal-dependent noise model by itself is not sufficient to characterize
the correlated noise (see Chapters 3 and 4), yet it is useful to obtain the noise curves of raw images,
where the noise is only-signal dependent. Section 1.1 explains in detail the procedure to adapt a
block-based noise estimator to measure signal-dependent noise.

Once the noise curve has been obtained, it might happen that some of its control points exhibit an overestimation for some particular intensities. This happens when all blocks in a small intensity range belong to a texture. Section 1.2 proposes a filtering algorithm to minimize the overestimation in the noise curve. Another problem that block-based noise estimators are expected to find is the presence of saturated points in the image. This special points appear at extreme intensities as isolated control points that distort the noise curve when interpolation is performed in between. Section 1.3 presents a procedure to avoid taking into account completely saturated pixels in the noise estimation.

The techniques presented in this chapter are general, in the sense that they can be applied
to almost any block-based noise estimator. In Chapters 10, 11, and 12, we apply them to the
Ponomarenko et al., Percentile and PCA methods, namely. The articles presented in this third
part were published in the Image Processing On Line (IPOL) journal. It publishes image processing
and image analysis algorithms, described in accurate literary form, coupled with code. This way,
scientists are allowed to check directly the published algorithms online with any uploaded image. It
also promotes reproducible research, and the establishment of a state of the art verifiable by all,
and on any image.

1. General techniques for adapt to signal-dependent noise estimation

1.1. Extension to signal-dependent noise. Most noise estimation methods found in the literature assume that the noise in the image is additive, signal-independent, and Gaussian. Note that in this context *uniform* means that the variance of the Gaussian noise is fixed and it does not depend on the intensity of the pixels of the ideal image. This assumption is not realistic because of the quantum nature of light itself and the way a CCD or CMOS detector responds to light. It is well-known that the emission of photons by a body follows a Poisson distribution. This distribution can be approximated by a Gaussian distribution when the number of photons is large enough. For very dark regions of the image this assumption does not hold. We consider
an image pixel $\tilde{U}(x, y)$ as a Poisson variable with variance and mean $U(x, y)$. The Poisson noise has therefore a standard deviation of $\sqrt{U(x, y)}$. (An image is nothing but a noise whose mean would be the ideal image.) This noise adds up to a thermal noise and to an electronic noise which are approximately additive and white. On a motionless scene with constant lighting, the expected value $U$ can be approximated by simply accumulating photons for a long exposure time, and then by taking the temporal average of this photon count. Any noise estimation algorithm assuming that the noise is uniform is unrealistic. Fortunately, most block-based methods are easily adapted to signal-dependent noise.

For a signal-dependent noise, a “noise curve” must be established. This noise curve associates with each image value $U(x, y)$ an estimation of the standard deviation of the noise associated with this value. Thus, for each block in the image, its mean must be computed and will give an estimation of a value in $U$. The measurement of the variation of the block (for example, its variance) will also be stored. The means are classified into a disjoint union of variable intervals or bins, in such a way that each interval contains a large enough number of elements. These measurements allow for the construction of a list of block variances whose corresponding means belong to the given bin. To find the number of samples/bin that minimizes the committed RMSE of the obtained noise curve compared to the ground-truth noise curve, we simulated signal-dependent noise with variance $\sigma^2 = 1 + 2U$ on a set of noise-free images (Figure 1) of 1080 × 808 pixels each. Then, the mean RMSE along all the images in the set was computed depending on the number of bins used. The number of bins that minimizes the RMSE depends on the method. For the Ponomarenko et al. and Percentile methods, the required number of samples/bin is 42000. Figure 2 shows the RMSE depending on the number of bins for the PCA method. The minimum for the PCA method is attained when using 5 bins. However, since the error committed when using 8 bins is not much worse than the error using 5 bins and a noise curve with 8 control points is more informative that the noise curve with 5 bins, we decided to use 8 bins for images of 1080 × 808 pixels. Therefore, the number of samples/bin is $\frac{1080 \times 808}{8} = 109080$. Experimental results with other noise-free natural images refined the minimum to 112000 samples/bin for the PCA method.

Therefore, it is possible to apply the generic noise estimator to each set of blocks associated with a given bin. In this way, an estimation of the noise for the intensities inside the limits of the bin is obtained. Because the set of bins is disjoint and there is no gap between bins, it is possible to deduce by interpolation a curve that relates the means of the blocks with their standard deviation, hence obtaining a signal-dependent *noise curve*. The intensity associated to each bin is given by the mean of the block at the percentile. The algorithmic description of the function building this histogram of block means can be found in Algorithm 21. This algorithm works as follows:

(1) It takes as input the number of bins that will be used (“bins” variable), the input data (the variances of the blocks, “data” variable), the associated intensities of the input data (the means of the blocks, “data” variable) and the total number of samples (“N” variable). The algorithm stores at the variable “samples_per_bin” the integer value of $N$/bins. In general, $samples_per_bin = 112000$ samples/bin. Since the last bin contain the remaining samples, it may contain less than $samples_per_bin$ samples.
<figure><img src="image_98.png" /><figcaption>FIGURE 1. Set of noise-free images used to determine empirically the optimal number of bins that the algorithm should use, as a function of the size of the image. Each image is 1080 × 808 pixels.</figcaption></figure>

<figure><img src="image_99.png" /><figcaption>FIGURE 2. Mean RMSE of the estimation of the signal-dependent noise with variance $\sigma^2 = 1 + 2U$ along all the noise-free test images in Figure 1, for the PCA method.</figcaption></figure>

(2) It returns for each bin $b$ its intensity bounds ("limits_begin[$b$]" and "limits_end[$b$]" variables), the list of variances that belong to bin $b$ ("data bins[$b$]" variable) and the list of intensities (block means) that belong to bin $b$ ("datal_bin[$b$]" variable).

(3) For each bin $b$, the algorithm fills the `data_bin[$b$]` and `datal_bin[$b$]` buffers with the variances and intensities of the blocks, sorted by their mean.

(4) The lower and upper intensity bounds of the current bin $b$ are stored into the variables `limits_begin[$b$]` and `limits_end[$b$]`. Then, the next bin is processed.
**1.2. Filtering the noise curve.** Optionally, the noise curve obtained on real images can be filtered. Indeed, it may present peaks when some given gray level interval contains mostly means of blocks belonging to a highly textured region. In this case, the measured block variance would be caused by the signal itself and not by the noise and the noise variance would be overestimated.

Given the *i*-th control point of the noise curve ($\hat{\mu}_i, \hat{\sigma}_i$) a closed intensity interval centered at this bin is considered, that is, [$\hat{\mu}_i - D, \hat{\mu}_i + D$]. For each intensity $\mu$ inside the interval (assuming that $\mu$ starts at $\hat{\mu}_i - D$ and it is incremented with a step of 0.05 while it is less or equal to $\hat{\mu}_i + D$), it is obtained the interpolated standard deviation that corresponds to each intensity $\mu$. In order to avoid an excessive interpolation, if $\hat{\mu}_i - D < \hat{\mu}_0$ for the *i*-th bin, then the diameter $D$ is changed to the value $\hat{\mu}_b - \hat{\mu}_0$. In the same way, if $\hat{\mu}_i + D > \hat{\mu}_{b-1}$ (being $B$ the number of bins), then the diameter $D$ is changed to the value $\hat{\mu}_{b-1} - \hat{\mu}_b$. Since each $\hat{\mu}_i$ can be seen as an oscillation (given by the RMSE) around the ideal value, averaging the noise curve inside the interval [$\hat{\mu}_i - D, \hat{\mu}_i + D$] for each control point attenuates the oscillations and puts them closer to the ground-truth. Once the oscillations have been attenuated, it might happen that a control point corresponds to a peak caused by a texture. In that case, the action taken is to compute the average inside the interval [$\hat{\mu}_i - D, \hat{\mu}_i + D$] and to substitute the standard deviation $\hat{\mu}_i$ of the *i*-th control point by the average only if it is lower than the average of the intensities in the interval. The filtering procedure is iterated five times. In the first three iterations the control points are allowed to go up and down, thus canceling the oscillations around the ideal value. In the next two iterations the points are only allowed to go down, to attenuate the overestimation of the noise because of textures. The simple strategy presented here performs properly for most natural images and in general not more than five filtering iterations are needed to get a reliable estimation of the noise. Applying more than five iterations does not improve the results significantly and for certain images it could produce noise curves that are excessively smooth. A diameter $D = 7$ is recommended. The pseudo-code of the filtering is detailed in Algorithm 24. It uses Algorithm 23 to interpolate the standard deviation corresponding to a given intensity. Algorithm 23 uses Algorithm 22 to get the corresponding standard deviation by a simple affine transformation.

Figure 3 shows the noise curve for the test image *Lena*. The non-filtered curve is drawn with solid lines and the filtered curve (five iterations) with dashed lines, using $D = 7$, $p = 0.5\%$, $w = 8$ and 6 bins. Note that the peak in the blue channel has decreased.

**1.3. Discarding saturated pixels.** When the number of photons measured by the CCD or CMOS detector during the exposure time is too high, its output may get saturated, and therefore underestimated. When the signal saturates the output of the CCD or CMOS detector, the measured variance in the saturated areas of the image is zero. Figure 4 shows an image with some saturated pixels. If the saturated pixels are taken into account when measuring the noise, the noise curve is no more reliable. Since the intensity of the saturated pixels is much higher than the intensity of most of the pixels in the image, there is usually a large gap between the values of normal non saturated pixels and the saturated ones and the noise curve will interpolate the standard deviation values. Of course, the information given by the noise curve inside this gap is
<figure><img src="image_100.png" /><figcaption>FIGURE 3. (a): test image *Lena* used to compare the noise curves with and without filtering; (b): noise curve for *Lena*. The non-filtered curve is drawn with solid lines and the filtered curve (five iterations) with dashed lines, using *D* = 7, *p* = 0.5%, *w* = 8 and 6 bins. Note that the peak in the blue channel is corrected after the filtering.</figcaption></figure>

<figure><img src="image_101.png" /><figcaption>FIGURE 4. Image with saturated pixels.</figcaption></figure>

not correct at all. In general, the strategy used to discard saturated pixels is to avoid the blocks
that contain a group of four connected exactly equal pixels, in any of the channels. This is useful
not only to discard saturated pixels, but also to avoid processing blocks whose pixels have suffered
other types of alterations that can be detected by finding these special blocks. For example, JPEG
encoding with a high compression factor sets to zero the value of the high-frequency coefficients in
many of the 8×8 blocks. Among other undesired effects like lower frequency artifacts and blocking
patterns, it can also create smooth zones, since high-frequency coefficients of the DCT of the block
were set to zero by the JPEG encoder. In natural images that have not been highly compressed,
the probability of finding a set of four connected pixels sharing exactly the same value is very low,
because of noise and texture. Nevertheless, we were unable to reproduce this problem using the
PCA algorithm, which does not seem to be affected by saturated pixels in raw images. However,
the option is left available to the user in the demo. The pseudo-code can be found in Algorithm
25. This algorithm checks if a pixel belong to a group of 2 × 2 pixels whose intensity is the same (in practice, with a difference between them below $\epsilon = 10^{-3}$). The blocks that contain at least one invalid pixels are not taken into account by the algorithm.
**Algorithm 21** Algorithm classifying blocks by their means.

1: **CLASSIFY_BY_MEAN** - Splits the input elements into disjoint *bins* according to the mean of the elements trying that each bin has the same cardinality. **Input bins**: number of bins. **Input data**: list of input data elements. **Input N**: number of elements/bin. **Input datal**: list of means of the input elements. **Output limits_begin[b]**: the lower intensity bound for bin *b*. **Output limits_end[b]**: the upper intensity bound for bin *b*. **Output data_bins[b]**: list of elements at bin *b*. **Output datal_bins[b]**: list of means of the elements at bin *b*.

2: samples_per_bin = ⌊N/bins⌋

3: limits_begin = zeros(bins)

4: limits_end = zeros(bins)

5: num_elements = zeros(bins)

6: data_bins = array(bins)

7: datal_bins = zeros(bins)

8: buffer = array(N)

9: bufferl = zeros(N)

▷ Sort data by datal

10: indices = argsort(datal, N)

▷ Min and max

11: min_datal = datal[indices[0]]

12: max_datal = datal[indices[N-1]]

13: lim0 = min_datal

14: elements_count = 0

15: bin = 0

16: **for** idx = 0 ... N **do**

17: **if** idx == N **then**

18: finished_loading = true

19: **else**

20: lim1 = datal[indices[idx]]

21: finished_loading = ¬ (bin == bins - 1) ∧ (elements_count ≥ samples_per_bin)

22: **end if**

23: **end for**

24: **if** finished_loading **then**

25: data_bins[bin] ← buffer

26: datal_bins[bin] ← bufferl

▷ Update limits and number of elements of the bin

27: limits_begin[bin] = lim0

28: limits_end[bin] = lim1

29: num_elements[bin] = elements_count

▷ Prepare for the next element

30: lim0 = lim1

31: bin = bin + 1

32: elements_count = 0

33: **else**

▷ Keep loading...

34: buffer[elements_count] = data[indices[idx]]

35: bufferl[elements_count] = datal[indices[idx]]

36: elements_count = elements_count + 1

37: **end if**

38: limits_end[bins-1] = max_datal
**Algorithm 22** Obtain a corresponding standard deviation by an affine transformation.

1: **AFFINE.** Input $(\mu_c, \sigma_c)$: current control point. Input $(\mu_e, \sigma_e)$: endpoint control point. Input $\mu$: intensity of the control points whose standard deviation is wanted. Output $\sigma$: standard deviation attributed to the intensity $\mu$.

2: $\varepsilon = 10^{-6}$

3: **if** $|\mu_c - \mu_e| < \varepsilon$ **then** ▷ Avoid dividing by zero

4: $s = 0$

5: **else**

6: $s = \frac{\sigma_c - \sigma_e}{\mu_c - \mu_e}$

7: **end if**

8: $\sigma = (\mu - \mu_e)s + \sigma_e$

**Algorithm 23** Interpolates an affine standard deviation from the points of the given noise curve.

1: **INTERPOLATION.** Input $(\mu_c, \sigma_c)$: known control points. Input $\mu$: the intensity of the point whose interpolated standard deviation is wanted. Output $\sigma$: the interpolated standard deviation of the point whose intensity is $\mu$.

▷ Find the nearest control point

2: $i = \operatorname{argmin}_i (\mu_c[i] - \mu)$

3: $m = \mu_c[i]$

4: **if** $\mu < m$ **then** ▷ on the right of $\mu$

5: **if** $i = 0$ **then** ▷ Treat boundary

6: $i = 1$

7: $m = \mu_c[i]$

8: **end if**

9: $m_1 = \mu_c[i - 1]$

10: $m_2 = m$

11: $s_1 = \sigma_c[i - 1]$

12: $s_2 = \sigma_c[i]$

13: **else** ▷ on the left of $\mu$

14: $N = \operatorname{len}(\mu_c)$

15: **if** $i \ge N - 1$ **then** ▷ Treat boundary

16: $i = N - 2$

17: $m = \mu_c[i]$

18: **end if**

19: $m_1 = m$

20: $m_2 = \mu_c[i + 1]$

21: $s_1 = \sigma_c[i]$

22: $s_2 = \sigma_c[i + 1]$

23: **end if**

return AFFINE($m_1, s_1, m_2, s_2, \mu$)
**Algorithm 24** Filters a noise curve.

1: **FILTER_CURVE** Input $(\mu_c, \sigma_c)$: list of control points to be filtered. **Input D**: diameter. **Input** allow_up: allow the points to go up and down. Otherwise, they are only allowed to go down. **Output** $\sigma_c^o$: returned list filtered standard deviations

2: $B = \text{len}(\mu_c)$

3: $\sigma_c^o \leftarrow \emptyset$

4: **for** $b = 0 \dots B - 1$ **do**

5: &nbsp;&nbsp;&nbsp;&nbsp;mu_current, std_current = $\mu_c[b], \sigma_c[b]$

6: &nbsp;&nbsp;&nbsp;&nbsp;left = mu_current - $D$

7: &nbsp;&nbsp;&nbsp;&nbsp;right = mu_current + $D$

▷ Adjust the diameter for the points near the boundary

8: **if** left < $\mu_c[0]$ **then**

9: &nbsp;&nbsp;&nbsp;&nbsp;dist = $\mu_c[b] - \mu_c[0]$

10: &nbsp;&nbsp;&nbsp;&nbsp;left = mu_current - dist

11: &nbsp;&nbsp;&nbsp;&nbsp;right = mu_current + dist

12: **else**

13: &nbsp;&nbsp;&nbsp;&nbsp;**if** right > $\mu_c[B - 1]$ **then**

14: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;dist = $\mu_c[B - 1] - \mu_c[b]$

15: &nbsp;&nbsp;&nbsp;&nbsp;left = mu_current - dist

16: &nbsp;&nbsp;&nbsp;&nbsp;right = mu_current + dist

17: **end if**

18: **end if**

▷ Add the interpolated control points inside the interval [left, right]

19: $\text{sum\_window} = 0$

20: $L = 0$

21: **for** $\mu$ = left...right (with step $\Delta = 0.05$) **do**

22: &nbsp;&nbsp;&nbsp;&nbsp;$\text{sum\_window} += \text{INTERPOLATION}(\mu_c, \sigma_c, \mu)$

23: $L += 1$

24: **end for**

25: $\text{std\_new} /= L$

26: **if** allow_up **then**

27: $\text{std\_filtered} = \text{std\_new}$

28: **else**

29: $\text{std\_filtered} = \text{std\_new}$ **if** $\text{std\_new} < \text{std\_current}$ **else** $\text{std\_current}$

30: **end if**

31: $\sigma_c^o \leftarrow \text{std\_filtered}$

32: **end for**

return $\sigma_c^o$
**Algorithm 25** Algorithm for the detection of groups of four equal pixels.

1: **DETECT_EQUAL_Pixels** - Creates a mask of valid pixels. **Input I**: input image. **Input N<sub>x</sub>**: width of **I**. **Input N<sub>y</sub>**: height of **I**. **Input w**: block side. **Input num\_channels**: number of channels of **I**. **Output mask**: mask of VALID/INVALID pixels.

2: $\varepsilon = 10^{-3}$

3: **for** i = 0... N<sub>x</sub> - 1 **do**

4:     **for** j = 0... N<sub>y</sub> - 1 **do**       ▷ Check if the pixel is not too close to the image boundary

5:     **if** i < N<sub>x</sub> - w + 1 ∧ j < N<sub>y</sub> - w + 1 **then**

6:         **for** c = 0 ... num\_channels - 1 **do**

7:             u = **I**.get\_channel(c)       ▷ Look if the 2 × 2 block is constant

8:             pixel\_status = (INVALID if c == 0 else mask[x,y])       ▷ Try to validate pixel

9:         **if** |u[i, j] - u[i + 1, j]| > ε ∨ |u[i + 1, j] - u[i, j + 1]| > ε ∨ |u[i, j + 1] - u[i + 1, j + 1]| > ε
10:     **then**

10:         pixel\_status = VALID

11:     **end if**

12:         mask[i, j] = pixel\_status

13:     **end for**

14:     **else**

15:         mask[i, j] = INVALID

16:     **end if**

17:     **end for**

18: **end for**
# The Ponomarenko et al. method

In the article An Automatic Approach to Lossy Compression of AVIRIS Images N.N. Ponomarenko et al. propose a new method to specifically compress AVIRIS (Airborne Visible/Infrared Imaging Spectrometer) images. As part of the compression algorithm, a noise estimation is performed with a proposed new algorithm based on the computation of the variance of overlapping 8 × 8 blocks. The noise is estimated on the high-frequency orthonormal DCT-II coefficients of the blocks. To avoid the effect of edges and textures, the blocks are sorted according to their energy measured on a set of low-frequency coefficients. The final noise estimation is obtained by computing the median of the variances measured on the high-frequency part of the spectrum of the blocks using only those whose energy (measured on the low-frequencies) is low. A small percentile of the total set of blocks (typically the 0.5%) is used to select those blocks with the lower energy at the low-frequencies. 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

1.1. Notation and Terminology. This chapter prepares the detailed description of the noise estimation algorithm given in Section 1.2 by fixing its notation and terminology.

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

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

* $N_x, N_y$: the width and height of $\tilde{U}$ in pixels.

* $\tilde{U}(x, y)$: the gray-level value of $\tilde{U}$ at pixel $(x, y)$, $x \in [0, N_x - 1]$ and $y \in [0, N_y - 1]$.

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

* $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 + 1)(N_y - w + 1)$.

* $\mathbf{D}_{x,y}$: the result of applying the orthonormal 2D DCT-II to a block $\mathbf{W}(x, y)$. Its coefficients are $\mathbf{D}_{x,y}(i, j)$ and the transform is defined as

$$ \mathbf{D}_{x,y}(i,j) = Q_w(i)Q_w(j) \sum_{n_x=0}^{w-1} \sum_{n_y=0}^{w-1} \mathbf{W}(x+n_x, y+n_y) \cos\left[\frac{\pi}{w}\left(n_x + \frac{1}{2}\right)i\right] \cos\left[\frac{\pi}{w}\left(n_y + \frac{1}{2}\right)j\right], $$
with $x \in [0, N_x - w - 1], y \in [0, N_y - w - 1], i \in [0, w - 1], j \in [0, w - 1]$ and $Q_N(k)$ is a normalization factor
$$Q_N(k) = \begin{cases} \frac{1}{\sqrt{N}}, & k = 0 \\ \sqrt{\frac{2}{N}}, & k \neq 0 \end{cases}.$$

From now on the DCT operator will refer specifically to this definition of the orthonormal 2D DCT-II.

* $\mathbf{W}_m$: a $w \times w$ pixels block in $\tilde{\mathbf{U}}$ 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 - w + 1} \rfloor, x = m - y(N_x - w + 1), m \in \{0, 1, \dots, M - 1\}$.

* $\mathbf{D}_m(i, j) = \text{DCT}(\mathbf{W}_m)$ where $m$ is the index of the block.

* $T$: threshold used by the function $\delta(i, j)$ that labels the coefficients of the transformed blocks $\mathbf{D}_{x,y}(i, j)$ as low-frequency coefficients (see Section 1.2.2 for more details).

## 1.2. The Algorithm.

1.2.1. Step 1: Computing the Set of Transformed Blocks {$\mathbf{D}_m(i, j)$}. From an image $\tilde{\mathbf{U}}$ of width $N_x$ and height $N_y$, corrupted with additive white Gaussian noise of variance $\sigma^2$, it is extracted a set of $M = (N_x - w + 1)(N_y - w + 1)$ (overlapping) $w \times w$ blocks {$\mathbf{W}_m$}, where $m$ is the index of the block, $m \in \{0, 1, \dots, 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, ...). Unlike other methods that pre-filter the image before extracting noise variance information from the blocks [57], the Ponomarenko et al. method measures the variance directly on $\mathbf{W}_m$. The DCT of each of these blocks is computed and gives the set {$\mathbf{D}_m$} of transformed blocks. The DCT coefficients in each block are denoted by $\mathbf{D}_m(i, j)$ where $m$ is the index of the block and $0 \le i, j < w$ is the frequency pair associated to that coefficient.

1.2.2. Step 2: Defining a Function to Label the Low/High Frequency Coefficients. The algorithm labels coefficients of the transformed blocks as belonging to low or medium/high frequencies. A coefficient corresponds to a low frequency if and only if $\delta(i, j) = 1$. If not, it is labeled as belonging to the medium/high frequencies set, where $\delta$ is defined by

$$\delta(i, j) = \begin{cases} 1, & (i + j < T) \land (i + j \neq 0), \rightarrow \text{low frequencies} \\ 0, & (i + j \ge T) \lor (i + j = 0) \rightarrow \text{medium/high frequencies.} \end{cases}$$

where $T$ is a given threshold, and $\land$ and $\lor$ stand for the AND and OR logical operators, namely. Note that this function does not label the mean of the block term $(i + j = 0)$ as a low-frequency.

1.2.3. Step 3: Estimating the Block Empirical Variance only with the Low-Frequency Coefficients. Given the set of transformed blocks {$\mathbf{D}_m$} with $m = 0, 1, \dots, M - 1$ the set of (empirical) variances associated to the low-frequency coefficients of the block $m$ is defined as

$$\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),$$

where $\theta = \sum_{i=0}^{w-1} \sum_{j=0}^{w-1} \delta(i, j)$ is the adequate normalization factor to get a mean.