# Multiscale Noise Estimation and Removal for Digital Images

Miguel Colom

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

## To cite this version:

Miguel Colom. Multiscale Noise Estimation and Removal for Digital Images. Image Processing [eess.IV]. Universitat de les Illes Balears (Espagne), 2014. English. NNT: tel-04124122

HAL Id: tel-04124122

https://hal.science/tel-04124122

Submitted on 19 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

**Universitat**
de les Illes Balears

# DOCTORAL THESIS

**Defense date: July 7, 2014**

MULTISCALE NOISE ESTIMATION AND
REMOVAL FOR DIGITAL IMAGES

Miguel Colom Barco

**Thesis Supervisors:**

Antoni Buades Capó and Jean-Michel Morel

**Jury:**

Bartomeu Coll Vicens (president),

Josep Lluís Lisani Roca (secretary)

Gabriele Facciolo Furlan (committee member)

Andrés Almansa (committee member)

Enric Meinhardt Llopis (committee member)

Bernard Rougé (invited)
A mi abuela, Baba.
A todos mis amigos en París.
This version of the document was compiled on 2023/06/19 at 18:08:34.

The last updated version of this document can always be downloaded from:

http://mcolom.info/download/thesis/thesis_miguel_colom.pdf

Author's contact

Email: miguel & mcolom.info (please change the ampersand by the @ sign).

Website: http://mcolom.info
Acknowledgements

Starting this thesis at the UIB has been quite a great adventure, but to be honest: also much fun and plenty of special moments with so many good friends I made along this adventure.

In the UIB I have to thank Joan and Cati for their warm welcome and friendship (you'll be the next ones writing your dissertations, after me. Good luck with it!). Also, Jonathan y Thibaud whom were doing a research internship.

I spent some days at CMLA to collaborate with the team at ENS-Cachan, so it was a great pleasure to have met Ariane, Aude, Barbara, Benjamin, Bruno, Carlo, Claire, Enric, Gabriele, Irène, Ives, José, Julie, Julien, Lara, Laurent, Loïc, Marc, Marie, Martín, Mauricio, Morgan, Nicola, Nicolas C., Nicolas L., Pauline, Rachel, Rafa, Samy, Saad, Tristan, Yohann, and Zhongwei. And also to Nick Chriss, for having improved my English with his classes, and all the secretary team: Véronique Almadovar, Micheline Brunetti, Sandra Doucet, Virginie Pauchont, and Carine Saint-Prix. And definitely, to everybody I've ever met at CMLA.

I have also to thank DxO-Labs, Centre National d'Etudes Spatiales, European Research Council, Office of Naval Research, Spanish Government, Direction Générale de l'Armement, Fondation Mathématique Jacques Hadamard, and Agence Nationale de la Recherche for funding some parts of my research.

I'm thankful to Vladimir V. Lukin, Alessandro Foi, and Benoit Vozel for their useful comments and interaction. It's been really helpful. Also, to Tomeu Coll for kindly looking at the papers I published in IPOL, suggesting changes and noticing mistakes. And to José Luis Lisani, for his thorough inspection of the papers I sent to IPOL and many useful comments, as part of the editorial process before publishing.

For sure, I really thank everybody I met at Colegio de España in Paris (Cité Universitaire) when I arrived to Paris to start a master before starting with the PhD. With not doubt, that's been one of the most fun and enriching experiences in my life! The list of friends I made there is too long to write their names in a single page.

Also, to the friends I made in Paris after I left Cité Universitaire, specially to Ceci, Jaime, Vicente, Cristina and little Luz. And all my friends and family in Palma, who understand my scientific adventures and always give me their support.

To my PhD supervisors Jean-Michel Morel and Toni Buades, for the weekly scientific discussions which gave rise to this thesis, for the projects I've been involved in, for the hours they've spent correcting all the drafts of my articles, and for being always so close. Also, to all the members of the commission that kindly accepted to evaluate this thesis.

And to Gloria, because you make every single moment the best.
# Contents

<table><tr><td>Summary</td><td>9</td></tr><tr><td>Preface</td><td>13</td></tr><tr><td>Bibliography of the thesis</td><td>25</td></tr><tr><td><b>Part 1. NOISE ESTIMATION</b></td><td><b>27</b></td></tr><tr><td>Chapter 1. The homoscedastic noise model</td><td>29</td></tr><tr><td>1. Introduction</td><td>29</td></tr><tr><td>2. Noise models</td><td>33</td></tr><tr><td>3. Review of homoscedastic block-based noise estimators</td><td>36</td></tr><tr><td>Chapter 2. Signal-dependent noise estimation</td><td>65</td></tr><tr><td>1. Introduction</td><td>65</td></tr><tr><td>2. State-of-the-art in white noise estimation</td><td>67</td></tr><tr><td>3. Non-parametric noise ground-truth curve</td><td>69</td></tr><tr><td>4. Non-parametric signal-dependent noise estimation</td><td>71</td></tr><tr><td>5. Cross-validation of several methods. Discussion</td><td>74</td></tr><tr><td>6. Conclusion</td><td>79</td></tr><tr><td>Chapter 3. The noise throughout the camera processing pipeline</td><td>81</td></tr><tr><td>1. The noise curves at each step of the camera pipeline</td><td>81</td></tr><tr><td>2. Overlapping of noise curves with different exposure times</td><td>95</td></tr><tr><td>3. Mean ground-truth curves</td><td>96</td></tr><tr><td>4. Comparison of the autocorrelation functions at different scales</td><td>97</td></tr><tr><td>5. Conclusions</td><td>102</td></tr><tr><td>Chapter 4. Multiscale estimation of intensity and frequency dependent noise</td><td>105</td></tr><tr><td>1. Introduction</td><td>105</td></tr><tr><td>2. Blind noise estimation principles</td><td>107</td></tr><tr><td>3. Noise estimation algorithm</td><td>110</td></tr><tr><td>4. Comparison</td><td>112</td></tr><tr><td>5. Validation with ground truth JPEG noise</td><td>116</td></tr><tr><td><b>Part 2. PATCH DENOISING</b></td><td><b>127</b></td></tr></table>
Chapter 5. Bayesian patch-based methods 129

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

Chapter 6. Generic tools for noise reduction 133

1. Aggregation of estimates 133

2. Iteration and “oracle” filters 134

3. Dealing with color images 136

4. Trying all generic tools on an example 136

Chapter 7. Detailed analysis of the Non-Local Means and the Non-local Bayes methods 141

1. Non-local means 141

2. Non-local Bayesian denoising 145

Chapter 8. The “Noise Clinic”: a universal denoiser 147

1. Introduction 147

2. A generalized nonlocal Bayesian algorithm 149

3. Obtaining the covariance matrix of noise patches 151

4. The multiscale algorithm 154

5. Validation 157

6. Results 161

7. Discussion 167

Part 3. REPRODUCIBLE RESEARCH CONTRIBUTIONS 177

Chapter 9. How to adapt homoscedastic noise estimators to signal-dependent noise 179

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

Chapter 10. The Ponomarenko et al. method 189

1. Noise Estimation Method 189

2. Online Demo 192

Chapter 11. The Percentile method 195

1. Introduction 195

2. Noise Estimation Method 196

3. Optimal Parameters 203

4. Online demo 203

Chapter 12. The PCA method 205

1. Noise estimation method 205

2. Principal components on natural images 206

3. Algorithm 207

4. Optimizing the PCA computation 207

5. Online demo 209

6. Appendix: proof of Theorem 2 210
CONTENTS

7

<table><tr><td>Chapter 13. Evaluation of the adapted methods</td><td>215</td></tr><tr><td>1. Evaluation with simulated white Gaussian noise</td><td>215</td></tr><tr><td>2. Evaluation comparing the noise curve of the raw image with the ground truth</td><td>218</td></tr><tr><td>3. Evaluation of the multiscale coherence of the result</td><td>224</td></tr><tr><td>4. Online demo</td><td>226</td></tr><tr><td>5. Complexity analysis of the algorithms</td><td>228</td></tr><tr><td>6. Conclusion</td><td>230</td></tr><tr><td>Final conclusion</td><td>233</td></tr><tr><td>Bibliography</td><td>235</td></tr></table>
<figure><img src="image_4.png" /><figcaption>FIGURE 1. Organization of the thesis.</figcaption></figure>
# Summary

## English

Any image, digital or analogic, contains not only information from the scene being photographed but also external interferences known as *noise*. The resulting image is the combination of the ideal image without noise with noise itself.

The “ideal image without noise” is a mathematic abstraction and it is not available in reality. Thus, it is needed methods that given only the degraded image are capable to properly characterize noise. This characterization using the noisy image is known as *blind noise estimation* since it does not use any additional information out the the noisy image.

Once noise has been properly characterized, the next step is to obtain a version of the image which is as close as possible to the ideal image. This process is known as *blind denoising*, since the ideal image is not available. Denoising methods exploit the property of *autosimilarity* of the small blocks that form the image to infer the geometry of the blocks of the ideal image. Denoising is a process guided by previous noise estimation.

Given that both noise estimation and denoising are performed *blindly*, it is important that noise characterization is as complete as possible. In this thesis several techniques for noise estimation are discussed, from the simplest which just consider homoscedastic noise, through those which consider the Poissonian model, to finally the new technique that we propose to obtain a complex noise model that depends on both intensity and frequency.

Regarding denoising, this thesis is mainly focuses on Bayesian techniques. The thesis finally reaches with the presentation of the *Noise Clinic*, the tool which we propose for automatic noise estimation and denoising. The Noise Clinic combines the automatic estimation of a complex noise model with its elimination at each of the scales of the image. This allows to restore a large typology of images, including those compressed with JPEG.

## Català

Qualsevol imatge, ja sigui digital o analògica, conté no només informació de l'escena fotografiada, sinó també interferències externes conegudes com *renou*. La imatge resultant és la combinació de la imatge ideal sense renou, amb el renou mateix.

La “imatge ideal sense renou” és una abstracció matemàtica i no està disponible a la realitat. Per tant, cal utilitzar mètodes que, donada únicament la imatge deteriorada, siguin capaços de caracteritzar adequadament el renou. Aquesta caracterització a partir de la imatge amb renou es
coneix com estimació a cegues del renou, atés que no s'utilitza cap altra informació addicional a part de la imatge amb renou.

Un cop caracteritzat adequadament el renou, el següent pas és obtenir una versió de la imatge que sigui tan fidel com sigui possible a la imatge ideal. Aquest procés es coneix com eliminació de renou a cegues, ja que la imatge ideal no està disponible. Els mètodes d'eliminació de renou aprofiten la propietat d' autosimilaritat dels petits blocs que componen la imatge per inferir la geometria dels blocs de la imatge ideal. L'eliminació de renou és un procés guiat per l'estimació de renou prèvia.

Atès que tant l'estimació com l'eliminació de renou es realitzen a cegues, és important que la caracterització del renou sigui tan completa com sigui possible. En aquesta tesi es discuteixen en detall les diverses tècniques per a l'estimació de renou, des de les més simples que únicament consideren renou homoscedàstic, passant per les que consideren el model poissonià de renou, fins a finalment la nova tècnica que proposem per obtenir un model de renou complex, que depèn tant de la intensitat com de la freqüència.

Pel que fa a l'eliminació de renou, aquesta tesi se centra especialment en les tècniques basades en el model bayesià. La tesi culmina amb la presentació de la *Noise Clinic*, l'eina que proposem per a l'estimació i eliminació automàtiques del renou. La Noise Clinic combina l'estimació automàtica d'un model de renou complex amb la seva eliminació en cadascuna de les escales de la imatge. Això permet restaurar una tipologia extensa d'imatges, incloent les comprimides amb JPEG.

## Castellano

Cualquier imagen, ya sea digital o analógica, contiene no solamente información de la escena fotografiada, sino también interferencias externas conocidas como *ruido*. La imagen resultante es la combinación de la imagen ideal sin ruido, con el propio ruido.

La "imagen ideal sin ruido" es una abstracción matemática y no está disponible en la realidad. Por lo tanto, es necesario utilizar métodos que, dada únicamente la imagen deteriorada, sean capaces de caracterizar adecuadamente el ruido. Esta caracterización a partir de la imagen ruidosa se conoce como *estimación a ciegas del ruido*, ya que no se utiliza ninguna otra información adicional aparte de la imagen ruidosa.

Una vez caracterizado adecuadamente el ruido, el siguiente paso es obtener una versión de la imagen que sea tan fiel como sea posible a la imagen ideal. Este proceso se conoce como *eliminación de ruido a ciegas*, ya que la imagen ideal no está disponible. Los métodos de eliminación de ruido aprovechan la propiedad de *autosimilaridad* de los pequeños bloques que componen la imagen para inferir la geometría de los bloques de la imagen ideal. La eliminación de ruido es un proceso guiado por la estimación de ruido previa.

Dado que tanto la estimación como la eliminación de ruido se realizan a ciegas, es importante que la caracterización del ruido sea tan completa como sea posible. En esta tesis se discuten en detalle las diversas técnicas para la estimación de ruido, desde las más simples que únicamente
consideran ruido homoscedástico, pasando por las que consideran el modelo poissoniano de ruido, hasta finalmente la nueva técnica que proponemos para obtener un modelo de ruido complejo, que depende tanto de la intensidad como de la frecuencia.

En cuanto a la eliminación de ruido, esta tesis se centra especialmente en las técnicas basadas en el modelo bayesiano. La tesis culmina con la presentación de la *Noise Clinic*, la herramienta que proponemos para la estimación y eliminación automáticas del ruido. La Noise Clinic combina la estimación automática de un modelo de ruido complejo con su eliminación en cada una de las escalas de la imagen. Esto permite restaurar una tipología extensa de imágenes, incluyendo las comprimidas con JPEG.
<figure><img src="image_5.png" /></figure>
# Preface

Digital images are matrices of regularly spaced pixels, each containing a photon count. This photon count is a stochastic process due to the physical quantum nature of light. It follows that all images are noisy. Ever since digital images exist, numerical methods have been proposed to improve the signal to noise ratio. Such “denoising” methods require a noise model and an image model. This thesis addresses the definition of noise models and their estimation from the digital image themselves. It also develops the main application which we call “blind denoising”, namely the fully automatic noise detection and removal. This is done in the framework of state of the art denoising algorithms which are mostly patch-based. For this reason, the thesis also presents a synthetic theory of patch-based methods.

This thesis is divided into three parts:

(1) Noise Estimation

(2) Patch Denoising

(3) Reproducible Research Contributions

## Part 1: NOISE ESTIMATION

In the first part of the thesis, Noise Estimation, we discuss several strategies to estimate the noise. The simplest experimental procedure to evaluate a noise estimation strategy just consists on simulating white Gaussian noise and adding it to a noise-free image (or to an image which is supposed to contain a very small or negligible noise). Then, the noise is estimated with some *homoscedastic*<fnref n="1" /> noise estimator and its variance is obtained (in this context, “homoscedastic” means that the variance of the noise does not vary depending on the intensity or the frequencies in the noisy image). In Chapter 1 several homoscedastic noise methods are presented and discussed. This work by Lebrun, Colom, Buades, and Morel was published in the Acta Numerica journal with the article *Secrets of image denoising cuisine* [1].

Unfortunately, the simple homoscedastic noise model is not useful to estimate the noise in real digital noise images. Indeed, the very first image acquired by the camera at the focal plane (the *raw* image) exhibits a noise that depends on the intensity. This noise, called the *photon noise* can be modeled with a Poisson distribution, for which the variance is an increasing function of the expectation (the mean intensity), also called a *noise curve*. This is related to the physical quantum

<footnote n="1">¹Here *homoscedastic* noise refers to a set of random variables with the same finite variance regardless their mean.</footnote>
nature of light, for which the emission of individual photons by any body is a Poisson random
process. Therefore, even in the raw image the noise is signal-dependent: the amount of noise
increases with the intensity of the underlying ideal image. This invalidates for real applications
the use of classic methods (discussed in Chapter 1) that only estimate a global variance of the
noise for the whole image. However, most of the noise estimation methods are patch-based and in
consequence they can be easily adapted to estimate signal-dependent noise. Although the exact
distribution for photon noise is the Poisson distribution, when the exposure time $\lambda$ is large enough
($\lambda > 1000$), the Poisson distribution of parameter $\lambda$ can be approximated with small error by
a Gaussian distribution with $\mu = \sigma^2 = \lambda$. Therefore, many signal-dependent noise estimation
methods assume that the noise is white and Gaussian for each intensity level [2, 3, 4, 5, 6].
However, the assumption that states that the variance is linear with the intensity is false in general,
since the saturation of the detectors at the most dark and bright pixels of the image gives a nonlinear
function of the variance according to the intensity. Even if the noise function under saturation can
be predicted quite accurately, as shown by Foi in the article *Practical Poissonian-Gaussian Noise
Modeling and Fitting for Single-Image Raw-Data* [7], the noise does not need to follow at all the
linear model and might follow any other model, as demonstrated by Boie and Cox in *An Analysis
of Camera Noise* [8]. In general, the information about the noise model that corrupted an image,
the characteristics of the detector and the exact transformations it suffered until the final image
was formed, are unknown. Therefore, in that situation the only solution is to obtain a noise model
directly from the noisy image, without assuming any prefixed model or parameters (*nonparametric
estimation*). In Chapter 2 we present to procedure to adapt any patch-based homoscedastic noise
estimator to obtain a signal-dependent noise estimation. In continuation we give a nonparametric
method that overcomes the state of the art in signal-dependent noise estimation for raw images. An
extensive cross-validation procedure is described to compare this new method with state-of-the-art
parametric methods and with laboratory calibration methods giving a reliable ground-truth, even
for nonlinear detectors. The procedure to obtain a ground-truth is described in detail. This work
was published by Colom, Buades, and Morel in the Journal of the Optical Society of America A,
with the article *Nonparametric Noise Estimation Method for Raw Images* [9].

In Chapter 3 we study in detail the characteristics of the noise through each step of the camera
processing chain, namely:

(1) The raw image acquisition at the focal plane of the camera.

(2) Demosaicing, the obtain a color image.

(3) White balance, the compensate the different gains of the detector at each channel in order
to get realistic colors.

(4) Gamma correction, to increase the dynamics of the image and therefore enhance the
visualization of dark pixels.

(5) JPEG compression, to reduce the size of the file that will finally contain the data of the
image by lossy compression.

The aim is to understand how each of these transformations affects the noise curve obtained
with a signal-dependent noise estimator. We used two different cameras (Canon EOS 30D and
<figure><img src="image_6.png" /><figcaption>FIGURE 2. Typical color spots or stains that can be observed in a JPEG image, caused mainly by two different steps of the camera processing. First, demosaicing correlates noise (thus creating low frequency noise that looks like color spots) and afterwards the gamma correction step increases the energy of the pixels, specially the darkest. As a result, the noise is converted in color spots which are clearly visible all over the image, specially at the darkest zones.</figcaption></figure>

Nikon D80), two ISO speeds (1250 and 1600), and four exposure times (1/30s, 1/250s, 1/400s, 1/640s). We identify and discuss the sources of the global perturbation that we observe as “noise” (dark noise, photon noise, readout noise, shot noise, and electronic noise), explain each of the steps in the camera processing pipeline, and for each of the step, discuss the obtained noise curves. Chapter 3 explains the origin of the kind of noise that is observed at the final JPEG image: small color spots, especially at the darkest zones. Figure 2 shows the typical color spots or stains that can be observed in a JPEG image, caused mainly by two different steps of the camera processing. First, demosaicing correlates noise (thus creating low frequency noise that looks like color spots) and afterwards the gamma correction step increases the energy of the pixels, specially at the darkest. As a result, the noise is converted in color spots which are clearly visible all over the image, specially at the darkest zones.

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 making it 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.

Figure 3 shows the noise curves obtained with a signal-dependent noise estimation along all the processing chain: raw image, demosaicing, white balance, gamma (tone curve) correction and JPEG compression, using a Canon EOS 30D camera. In solid lines, the temporal estimation (ground-truth) and in dashed lines the spatial estimation.

In fact, if the noise estimator assumes that the noise is signal-dependent but does not take into account that the noise depends on the frequency, the noise is strongly subestimated, as can be
<figure><img src="image_7.png" /><figcaption>FIGURE 3. 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 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. In solid lines, the temporal estimation (ground-truth) and in dashed lines the spatial estimation.</figcaption></figure>

observed in Figure 3, where the spatial estimation given by the algorithm (dashed lines) is really underestimating the noise (solid lines, the ground truth from the temporal estimation). 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 on the intensity, frequency and scale from the noisy image itself. Chapter 3 concludes that any noise estimation algorithm that is intended for real images must consider the noise not only signal-dependent, but also frequency-dependent (SFD noise). The noise must be estimated according both the intensity and the frequency.

Chapter 4 presents our proposed method to estimate SFD noise. We present a new semi-distance to measure the likeliness between two patches, called the *sparse distance* and depict a new algorithm that overcomes the current state of the art algorithm for frequency-dependent noise. The algorithm is validated using both simulations and observations of denoising results on real images. This work by Colom, Lebrun, Buades, and Morel is submitted to the IEEE Transactions On Image Processing with the article *Multiscale Estimation of Intensity and Frequency Dependent Noise* [10].
## Part 2: PATCH DENOISING

The second part of the thesis, Patch Denoising, discusses the obvious application of noise estimation: using the noise model obtained from the noisy image itself to remove the noise from it and obtain a new version of the image for which the noise has been removed (or at least, minimized). Of course, details, textures, and edges must be preserved. This process is called denoising.

Chapter 5 discusses the Bayesian patch-based method, which gives 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 $u$ with dimension $\kappa \times \kappa$, and $\tilde{P}$ an observed noisy version of $P$, the model gives by assuming a Gaussian model and the independence of noise pixel values

$$ \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}$, the goal is to deduce $P$ by maximizing $\mathbb{P}(P|\tilde{P})$. Using Bayes' rule, the last conditional probability can be written as

$$ \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$. If we assume that patches $Q$ similar to a given patch $P$ also follow the Gaussian model with empirical sample covariance matrix $C_P$ and empirical sample mean $\bar{P}$, then

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

Using the classical Bayesian equivalence of problems,

$$ \max_P \mathbb{P}(P|\tilde{P}) \Leftrightarrow \max_P \mathbb{P}(\tilde{P}|P)\mathbb{P}(P) $$

and after some calculus, Chapter 5 concludes that a restored patch $\hat{P}_1$ can be obtained from the observed patch $\tilde{P}$ by the one step estimation

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

Chapter 6 discusses and gives also the detailed algorithmic descriptions of three generic tools used to denoise an image:

* Aggregation

* “Oracle”

* Color transform

The idea behind aggregation techniques is to 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$. 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).

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. Another type of aggregation technique, like in the works of Raphan, Van De Ville, and Deledalle, considers the risk estimate rather than the variance to locally attribute more weight to the estimators with small risk [11, 12, 13].

The use of an “oracle” is another technique used to improve the denoising results. 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.

A color transform is also a useful tool to avoid artifacts introduced by denoising. Indeed, most denoising algorithms treat each image channel independently, which may introduce color artifacts easily noticeable by the eye. Denoising the image in a different color space avoids this problem. For example, passing from the RGB to the YUV colorspace. Finally, Chapter 6 shows an example (Figure 1) where it can be clearly identified the effect of each of these tools when performing denoising.

Chapter 7 presents and discusses in detail two patch-based denoising methods, the classic NL-means method which uses similar patches to denoise by aggregation and the NL-Bayes method (Chapter 7), which uses an oracle step to obtain a list of similar patches for each patch it denoises. A cleaned version of the patch is obtained from the improved sample covariance matrix of the similar patches thanks to the oracle. Detailed algorithmic descriptions for both methods are given in Algorithm 18 and Algorithm 19, namely.

Chapter 8 develops the project which we named *Noise Clinic*, that involves applying the multi-scale intensity-frequency estimation method in Chapter 4 to obtain a noise model from the patches of the noisy image itself (nonparametric estimation) to denoise it without any prior information (blind denoising) using a Bayesian patch-based method. The method uses a *multiscale* strategy that allows to estimate accurately strongly correlated low-frequency noise, as long as the image is large enough.

Indeed, 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 effectuated by uncontrolled hardware and software (and sometimes by chemical means). For such images, recent progress in noise estimation permits to estimate from a single image a noise model which is simultaneously signal and frequency dependent. Chapter 8 proposes a multiscale denoising algorithm adapted to this broad noise model. 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. The consistency of this algorithm is also verified on simulated distorted images. This algorithm is finally compared to the unique state of the art previous blind denoising method [14] by Portilla.

Figure 4 shows a 7 February 1966 picture of the Honolulu Conference, with General Earle Wheeler, Secretary of State Dean Rusk, and President Lyndon B. Johnson depicted in Camp Smith, Hawaii<fnref n="2" /> (a), a detail in the noisy image (b), and the same detail in the denoised image (c). It can be observed that the noise has been removed, while keeping the details of the image. Figure 5 shows the noise curves obtained when denoising the image using two scales with the Noise Clinic. On the left, the mean of the standard deviations of the high (a) and low (b) frequencies at the first scale. On the right, the mean of the standard deviations of the high (b) and low (d) frequencies at the second scale. It can be observed that the noise at the first scale is lower than the noise at the second scale, and that within any of the scales the noise at the low frequencies is higher than the noise at the high frequencies. This is because the noise has mainly energy in its low frequencies and is highly correlated. Figure 6 shows a detail of difference between the noisy image and the image denoised at each scale. In scale #1 (b), the noise has a particular spatial structure, since low-frequency noise is detected at the second scale. On scale #0 (a), the noise is less correlated, but still mainly low-frequency noise. To improve the visualization of the images, the histogram of the difference image has been equalized.

## Part 3: REPRODUCIBLE RESEARCH CONTRIBUTIONS

The third part of the thesis presents three of the reproducible research contributions of this dis- sertation:

* how to adapt block-based homoscedastic noise estimators to measure intensity-dependent noise;

* how to filter the obtained noise curves to cancel undue oscillations;

* how saturated pixels distort the shape of noise curve and how to avoid them.

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 ac- curate 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.

Section 1.1 of Chapter 9 discusses how to adapt most of the patch-based noise estimation methods to measure intensity-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

<footnote n="2">²From the National Archives and Records Administration. The National Archives and Records Administration provides images depicting American and global history which are public domain or licensed under a free license.</footnote>
<figure><img src="image_8.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(b)</figcaption></figure>

<caption>(c)</caption>

FIGURE 4. A 7 February 1966 picture of the Honolulu Conference, with General Earle Wheeler, Secretary of State Dean Rusk, and President Lyndon B. Johnson depicted in Camp Smith, Hawaii (a), a detail in the noisy image (b), and the same detail in the denoised image (c). It can be observed that the noise has been removed, while keeping the details of the image.

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. The procedure to obtain the noise curve is discussed and the description is given in Algorithm 21.

The noise curve obtained 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
<figure><img src="image_10.png" /><figcaption>FIGURE 5. Noise curves obtained when denoising the image with two scales. On the left, the mean of the standard deviations of the high (a) and low (b) frequencies at the first scale. On the right, the mean of the standard deviations of the high (b) and low (d) frequencies at the second scale.</figcaption></figure>

be overestimated. To solve it, the noise curve obtained can be filtered. The pseudocode of the filtering is detailed in Algorithm 24.

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. If saturated pixels are taken into account when measuring the noise, the noise curve is no more reliable. Section 1.3 of Chapter 9 presents the strategy proposed to discard saturated pixels, that consists on rejecting the blocks that contain a group of four connected exactly equal pixels, in any of the channels. The pseudocode can be found in Algorithm 25.

These three general tools (adaptation to signal-dependent noise, noise curve filtering, and avoiding saturated pixels in the estimation) were applied to the

* Ponomarenko et al. method;
<figure><img src="image_11.png" /><figcaption>FIGURE 6. Detail of difference between the noisy image and the image denoised at each scale. In scale #1 (b), the noise has a particular spatial structure, since low-frequency noise is detected at the second scale. On scale #0 (a), the noise is less correlated, but still mainly low-frequency noise. To improve the visualization of the images, the histogram of the difference image has been equalized.</figcaption></figure>

• Percentile method;

• PCA method.

Chapter 10 discusses and analyzes in deep detail the Ponomarenko et al. method inside this new framework. In the article *An Automatic Approach to Lossy Compression of AVIRIS Images* [15, 16] N. N. Ponomarenko, V. V. Lukin, M. S. Zriakhov, A. Kaarna, and J. T. Astola 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.

Chapter 11 discusses and analyzes in depth the Percentile method [17] using the presented framework. Given a white Gaussian noise signal $N_\sigma$ on a sampling grid, its variance $\sigma^2$ can be estimated from a small $w \times w$ pixels sample. However, in natural images we observe $\tilde{U} = U + N_\sigma$, the combination of the geometry of the scene that is photographed and the added noise. In this case, estimating directly the standard deviation of the noise from $w \times w$ samples of $\tilde{U}$ is not reliable since the measured standard deviation is not explained just by the noise but also from
the geometry of **U**. The Percentile method tries to estimate the standard deviation σ from w × w
blocks of a high-passed version of **Ũ** by a small p-percentile of these standard deviations. The idea
behind is that edges and textures in a block of the image increase the observed standard deviation
but they never make it decrease. Therefore, a small percentile (0.5%, for example) in the list of
standard deviations of the blocks is less likely to be affected by the edges and textures than a
higher percentile (50%, for example). The 0.5%-percentile is empirically proven to be adequate for
most natural, medical and microscopy images.

Chapter 12 discusses and analyzes in deep detail the PCA method from the article *Image Noise Level Estimation by Principal Component Analysis* [18, 3], where S. Pyatykh, J. Hesser, and L. Zheng propose a new method to estimate the variance of the noise in an image from the eigenvalues of the covariance matrix of the overlapping blocks of the noisy image. Instead of using all patches of the noisy image, the authors propose an iterative strategy to adaptively chose the optimal set containing the patches with lowest variance. The method is analyzed inside the presented framework.

For these three methods, the following tests were performed in Chapter 13 for the Ponomarenko
et al., Percentile and PCA methods:

* Tests on simulated white Gaussian noise using the noise-free images. In this case we took seven bins to classify the blocks according to their means (see Section 1.1 of Chapter 9).

* Tests on a set of real raw images obtained by a Canon EOS 30D camera (see Figure 2). The procedure explained in Section 1.1 of Chapter 9 was used to get a noise curve. The results were compared to the ground-truth noise curve of the camera.

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

We show the results obtained with the set of images for the Ponomarenko et al., Percentile and PCA methods.

Of course, estimating homoscedastic noise is not enough to characterize the signal-dependent
noise in real digital images. In order to evaluate the accuracy of the three methods with signal-
dependent noise, we obtained a ground-truth noise curve of the Canon EOS 30D camera with the
procedure explain in Chapter 2. Then, we obtained the noise curves for each of the raw images
in 2 and compared them with the ground-truth noise curve of the camera. Figure 7 shows the
ground-truth noise curve for raw images for the Canon EOS 30D camera, with ISO speed 1600.

To evaluate the accuracy of each method, we compared the obtained noise curve obtained by
each method with the ground-truth.

After evaluating all three methods with the three tests (simulated homoscedastic noise, intensity-
dependent noise, and multiscale coherence), the conclusion is that the strategy followed by the
modified Ponomarenko et al. method is the best and can be considered as the state of the art
in intensity-dependent noise estimation. The Percentile method gives slightly worse results than
<figure><img src="image_12.png" /><figcaption>FIGURE 7. Ground-truth noise curve for the Canon EOS 30D camera, for ISO speed 1600.</figcaption></figure>

Ponomarenko et al. method and PCA is in general worse and requires at least the double of samples than Ponomarenko and Percentile to achieve similar results.

Section 5 of Chapter 13 discusses about the algorithmic complexity the operations needed to adapt homoscedastic noise estimation methods into intensity-dependent, the noise filtering, and the detection of blocks with saturated pixels. These operations can be considered a common framework where almost all noise estimation algorithms can be put inside and estimate intensity-dependent noise. The algorithmic complexity analysis of the Ponomarenko et al., Percentile, and PCA methods within the common framework is also given.

Chapter 6 presents the conclusions of the research presented in this thesis.
# Bibliography of the thesis

The following research articles were published or submitted during this thesis:

## Published

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

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

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

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

## To appear

M. Colom, G. Facciolo, M. Lebrun, N. Pierazzo, M. Rais, Y. Wang, and J.M. Morel, "A mathematical perspective of image denoising", *International Congress of Mathematicians*, 2014, To appear

## Submitted

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

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

M. Colom and A. Buades, "Analysis and Extension of the PCA Method, Estimating a Noise Curve from a Single Image", *Image Processing On Line*, 2014, Submitted

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

M. Lebrun, M. Colom, and J.M. Morel, "Multiscale image blind denoising", *IEEE Transactions On Image Processing*, 2014, Submitted
<figure><img src="image_13.png" /></figure>
Part 1

NOISE ESTIMATION
In this first part of the thesis, we address the problem of the different noise models and how ade-
quate they are to estimate noise in different scenarios. From the simplest models for white Gaussian
noise to the most complex noise model which are able to estimate noise in JPEG-compressed im-
ages. All the noise estimation models described in this thesis (and thus, all denoising methods),
are block-based.

Chapter 1 discusses the homoscedastic white Gaussian noise model, the simplest model. It is not directly applicable to real images, where the physics of light make the variance of the noise depend on the intensity. However, it is possible to adapt the noise estimator under this model to measure signal-dependent noise, as explained in Chapter 2 and hence the interest of considering this model. In this chapter we evaluate several homoscedastic noise estimation methods and also introduce the objective of noise estimation which is, obviously, denoising. The classic NL-means denoiser is shown as an example, and three generic tools for denoising are discussed.

In Chapter 2 we discuss the adaptation of block-based noise estimators to measure signal-
dependent noise, a model that is consistent with the photon noise (modeled with the Poisson dis-
tribution). The concept of “noise curve” is introduced and a non-parametric estimating the noise
curve directly from a single raw image is described. We also detail the procedure that permits to
obtain a ground-truth noise curve.

In order to check the validity of the signal-dependent noise model, Chapter 3 shows the noise curves at each step of the camera processing pipeline (from the initial raw to the final JPEG image) and compares them with the corresponding ground-truth curve. The signal-dependent noise model is valid for raw images, but is not enough for correlated noise, as shown in Chapter 3 This noise curves along the step of the camera pipeline show that, indeed, the signal-dependent model is not sufficient after demosaicing.

Chapter 4 discusses a new estimation method able to measure the noise even in JPEG-encoded
images, where the noise model is complex, with low-frequency noise, and highly correlated.
# The homoscedastic noise model

In this chapter, our main goal is to review not less than 13 blind homoscedastic white noise estimation methods using a single image. These methods will be discussed within the context of block-based denoising. Indeed, the main idea behind denoising techniques is to exploit the self-similarity property of natural images, to find similar patches whose aggregation gives a denoised version of them. The measure of the distance between patches depends on the variance of the noise and therefore noise estimation is a needed step for blind denoising. We will show this using a classic denoiser (NL-means) as an example. The performance of all methods depends on three generic tools: color transform, aggregation, and an “oracle” step, which will be presented in this chapter.

In this chapter we will only focus on homoscedastic noise estimation, that is, assuming that the noise variance is fixed and does not depend on the intensity or the frequency of the underlying noise-free image. Since adapting block-based noise estimators to signal and frequency dependent noise is relatively easy, it is justified to analyze first homoscedastic noise estimators. Later in Chapter 2 we discuss the adaptation of block-based noise estimator to signal-dependent noise and in Chapter 4 the adaptation to signal and frequency dependent noise.

## 1. Introduction

Most digital images and movies are currently obtained by a CCD or CMOS detector. The value $\tilde{u}(\mathbf{i})$ observed by a sensor at each pixel $\mathbf{i}$ is a Poisson random variable whose mean $u(\mathbf{i})$ would be the ideal image. The difference between the observed image and the ideal image $\tilde{u}(\mathbf{i}) - u(\mathbf{i}) = n(\mathbf{i})$ is called “shot noise”. On a motionless scene with constant lighting, $u(\mathbf{i})$ can be approached by simply accumulating photons for a long exposure time, and by taking the temporal average of this photon count, as illustrated in Figure 1.

Accumulating photon impacts on a surface is therefore the essence of photography. The first Nicéphore Niépce photograph in 1826 [23] was obtained after an eight hours exposure. The problem of a long exposure is the variation of the scene due to changes in light, camera motion, and incidental motions of parts of the scene. The more these variations can be compensated, the longer the exposure can be, and the more the noise can be reduced. If a camera is set to a long exposure time, the photograph risks motion blur. If it is taken with short exposure, the image is dark, and enhancing it reveals the noise.

A possible solution is to take a burst of images, each with short-exposure time, and to average them after registration. This technique, illustrated in Figure 1, was evaluated recently in a paper
that proposes fusing bursts of images taken by cameras [24]. This paper shows that the noise reduction by this method is almost perfect: fusing *m* images reduces the noise by a $\sqrt{m}$ factor.

It is not always possible to accumulate photons. There are obstacles to this accumulation in astronomy, biological imaging and medical imaging. In day to day images, the scene is moving, which limits the exposure time. The main limitations to any imaging system are therefore the noise and the blur. In this review, experiments will be conducted on photographs of scenes taken by normal cameras. Nevertheless, the image denoising problem is a common denominator of all imaging systems.

A naive view of the denoising problem would be: how to estimate the ideal image, namely the mean $u(\mathbf{i})$, given only one sample $\tilde{u}(\mathbf{i})$ of the Poisson variable? The best estimate of this mean is of course this unique sample $\tilde{u}(\mathbf{i})$. Getting back a better estimate of $u(\mathbf{i})$ by observing only $\tilde{u}(\mathbf{i})$ is impossible. Getting a better estimate by using also the rest of the image is obviously an ill-posed problem. Indeed, each pixel receives photons coming from different sources.

Nevertheless, a glimpse of a solution comes from image formation theory. A well-sampled image $u$ is band-limited [25]. Thus, it seems possible to restore the band-limited image $u$ from its degraded samples $\tilde{u}$, as was proposed in 1966 in [26]. This classic Wiener-Fourier method consists in multiplying the Fourier transform by optimal coefficients to attenuate the noise. It results in a convolution of the image with a low-pass kernel.

From a stochastic viewpoint, the band-limitedness of $u$ also implies that values $\tilde{u}(\mathbf{j})$ at neighboring pixels $\mathbf{j}$ of a pixel $\mathbf{i}$ are positively correlated with $\tilde{u}(\mathbf{i})$. Thus, these values can be taken into account to obtain a better estimate of $u(\mathbf{i})$. These values being nondeterministic, Bayesian approaches are relevant and have been proposed as early as 1972 in [27].

In short, there are two complementary early approaches to denoising, the Fourier method, and the Bayesian estimation. The Fourier method has been extended in the past thirty years to other linear space-frequency transforms such as the windowed DCT [28] or the many wavelet transforms [29].

Being first parametric and limited to rather restrictive Markov random field models [30], the Bayesian method are becoming non-parametric. The idea for the recent non parametric Markovian estimation methods is a now famous algorithm to synthesize textures from examples [31]. The underlying Markovian assumption is that, in a textured image, the stochastic model for a given pixel $\mathbf{i}$ can be predicted from a local image neighborhood $P$ of $\mathbf{i}$, which we shall call “patch”.

The assumption for recreating new textures from samples is that there are enough pixels $\mathbf{j}$ similar to $\mathbf{i}$ in a texture image $\tilde{u}$ to recreate a new but similar texture $u$. The construction of $u$ is done by nonparametric sampling, amounting to an iterative copy-paste process. Let us assume that we already know the values of $u$ on a patch $P$ surrounding partially an unknown pixel $\mathbf{i}$. The Efros-Leung [31] algorithm looks for the patches $\tilde{P}$ in $\tilde{u}$ with the same shape as $P$ and resembling $P$. Then a value $u(\mathbf{i})$ is sorted among the values predicted by $\tilde{u}$ at the pixels resembling $\mathbf{j}$. Indeed, these values form a histogram approximating the law of $u(\mathbf{i})$. This algorithm goes back
<figure><img src="image_14.png" /><figcaption>FIGURE 1. From left to right: (a) one long-exposure image (time=0.4 s, ISO=100), one of 16 short-exposure images (time=1/40 s, ISO=1600) and their average after registration. The long exposure image is blurry due to camera motion. (b) The middle short-exposure image is noisy. (c) The third image is about four times less noisy, being the result of averaging 16 short-exposure images. From [24].</figcaption></figure>

to Shannon's theory of communication [25], where it was used for the first time to synthesize a probabilistically correct text from a sample.

As was proposed in [32], an adaptation of the above synthesis principle yields an image denoising algorithm. The observed image is the noisy image $\tilde{u}$. The reconstructed image is the denoised image $\hat{u}$. The patch is a square centered at $i$, and the sorting yielding $u(i)$ is replaced by a weighted average of values at all pixels $\tilde{u}(j)$ similar to $i$. This simple change leads to the "non-local means" algorithm, which can therefore be sketched in a few rows.

<caption>**Algorithm 1 Non-local means algorithm**</caption>

1: **Input:** noisy image $\tilde{u}$, $\sigma$ noise standard deviation. **Output:** denoised image $\hat{u}$.
2: Set parameter $\kappa \times \kappa$: dimension of patches.
3: Set parameter $\lambda \times \lambda$: dimension of research zone in which similar patches are searched.
4: Set parameter $C$.
5: **for** each pixel **i** **do**
6: &nbsp;&nbsp;&nbsp;&nbsp;Select a square reference sub-image (or "patch") $\tilde{P}$ around **i**, of size $\kappa \times \kappa$.
7: &nbsp;&nbsp;&nbsp;&nbsp;Call $\hat{P}$ the denoised version of $\tilde{P}$ obtained as a weighted average of the patches $\tilde{Q}$ in a square neighborhood of **i** of size $\lambda \times \lambda$. The weights in the average are proportional to $w(\tilde{P}, \tilde{Q}) = e^{-\frac{d^2(\tilde{P}, \tilde{Q})}{2C\sigma^2}}$ where $d(\tilde{P}, \tilde{Q})$ is the Euclidean distance between patches $\tilde{P}$ and $\tilde{Q}$.
8: **end for**
9: Aggregation: recover a final denoised value $\hat{u}(i)$ at each pixel $i$ by averaging all values at $i$ of all denoised patches $\hat{Q}$ containing $i$.

It was also proved in [32] that the algorithm gave the best possible mean square estimation if the image was modeled as an infinite stationary ergodic spatial process (see Chapter 7 for an exact statement). The algorithm was called "non-local" because it uses patches $\tilde{Q}$ that are far away from $\tilde{P}$, and even patches taken from other images. NL-means was not the state of the art denoising method when it was proposed. As we shall see in the comparison Section 4 of
Chapter 4, the 2003 Portilla et al. [33] algorithm has a better PSNR performance. But quality criteria show that NL-means creates less artifacts than wavelet based methods. This may explain why patch-based denoising methods have flourished ever since. By now, 1500 papers have been published on nonlocal image processing. Patch-based methods seem to achieve the best results in denoising. Furthermore, the quality of denoised images has become excellent for moderate noise levels. Patch-based image restoration methods are used in many commercial software.

An exciting recent paper in this exploration of nonlocal methods raises the following claim [34]: For natural images, the recent patch-based denoising methods might well be close to optimality. The authors use a set of 20000 images containing about $10^{10}$ patches. This paper provides a second answer to the question of absolute limits raised in [35], “Is denoising dead?”. The Cramer-Rao type lower bounds on the attainable RMSE performance given in [35] are actually more optimistic: they allow for the possibility of a significant increase in denoising performance. The two types of performance bounds considered in [34] and [35] address roughly the same class of patch-based algorithms. It is interesting to see that these same authors propose denoising methods that actually approach these bounds, as we shall see in Chapter 7.

The denoising method proposed in [34] is actually based on NL-means (algorithm 1), with the adequate parameter $C$ to account for a Bayesian linear minimum mean square estimation (LMMSE) estimation of the noisy patch given a database of known patches. The only and important difference is that the similar patches $Q$ are found on a database of $10^{10}$ patches, instead of on the image itself. Furthermore, by a simple mathematical argument and intensive simulations on the patch space, the authors are able to approach the best average estimation error which will ever be attained by any patch-based denoising algorithm [34].

These optimal bounds are nonetheless obtained on a somewhat restrictive definition of patch-based methods. A patch-based algorithm is understood as an algorithm that denoises each pixel by using the knowledge of: a) the patch surrounding it, and b) the probability density of all existing patches in the world. It turns out that state of the art patch-based denoising algorithms use more information taken in the image than just the patch. For example, most algorithms use the obvious but powerful trick to denoise all patches, and then to aggregate the estimation of all patches containing a given pixel to denoise it better. Conversely, these algorithms generally use much less information than a universal empirical law for patches. Nevertheless, the observation that at least one algorithm, BM3D [36] might be arguably very close to the best predicted estimation error is enlightening. Furthermore, doubling the size of the patch used in the [34] paper would be enough to cover the aggregation step. The difficulty is to get a faithful empirical law for $16 \times 16$ patches.

The “convergence” of all algorithms to optimality will be corroborated here by the thorough comparison of nine recent algorithms (Section 4 of Chapter 4). These state of the art algorithms seem to attain a very similar qualitative and quantitative performance. Although they initially seem to rely on different principles, our final discussion will argue that these methods are equivalent.

Image restoration theory cannot be reduced to an axiomatic system, as the statistics of images are still a widely unexplored continent. Therefore, a complete theory, or a single final algorithm
closing the problem are not possible. The problem is not fully formalized because there is no
rigorous image model. Notwithstanding this limitation, rational recipes shared by all methods can
be given, and the methods can be shown to rely on only very few principles. More precisely, this
chapter will present the following recipes, and compare them whenever possible:

* three techniques that improve every denoising method (Chapter 6);

* a complete review on several families of homoscedastic noise estimation techniques (Section 3);

Nevertheless, this convergence of results and techniques leaves several crucial issues unsolved. (This is fortunate, as no researcher likes finished problems.) With one exception, (the BLS-GSM algorithm [37, 14, 38]), state of the art denoising algorithms are not multiscale. High noises and small noises also remain unexplored.

In a broader perspective, the success of image denoising marks the discovery and exploration of one of the first densely sampled high dimensional probability laws ever (numerically) accessible to mankind: the “patch space”. For 8 × 88 patches, by applying a local PCA to the patches surrounding a given patch, one can deduce that this space has a dozen significant dimensions (the others being very thin). Exploring its structure, as was initiated in [39], seems to be the first step toward the statistical exploration of images. But, as we shall see, this local analysis of the *patch space* already enables state of the art image denoising.

Most denoising and noise estimation algorithms commented here are available at the *Image Processing on Line* (IPOL) journal, http://www.ipol.im/. Each algorithm is given a complete description, the corresponding source code, and can be run online on arbitrary images.

## 2. Noise models

The main source of noise in digital images is the so-called *shot noise*, which is inherent to photon counting. The value $\tilde{u}(i)$ observed by a sensor at each pixel $i$ is a Poisson random variable whose mean would be the ideal image. The standard deviation of this Poisson distribution is equal to the square root of the number of incoming photons $\tilde{u}(i)$ in the pixel captor $i$ during the exposure time. This noise adds up to a thermal noise and to an electronic noise which are approximately additive and white.

For sufficiently large values of $\tilde{u}(i)$, ($\tilde{u}(i) > 1000$), the normal distribution $N(\tilde{u}(i), \sqrt{\tilde{u}(i)})$ with mean $\tilde{u}(i)$ and standard deviation $\sqrt{\tilde{u}(i)}$ is an excellent approximation to the Poisson distribution. If $\tilde{u}(i)$ is larger than 10, then the normal distribution still is a good approximation if an appropriate continuity correction is performed, namely $P(\tilde{u}(i) \le a) \approx P(\tilde{u}(i) \le a + 0.5)$, where $a$ is any non-negative integer.

As a rule of thumb, the noise model is relatively easy to estimate when the raw image comes directly from the imaging system, in which case the noise model is known and only a few parameters must be estimated. For this, efficient methods are described by Foi et al. [40, 7] for Poisson and Gaussian noise.
Nevertheless, the pixel value is signal dependent, since its mean and variance depend on $\tilde{u}(i)$. To get back to the classic “white additive Gaussian noise” used in most researches on image denoising, a variance-stabilizing transformation can be applied: when a variable is Poisson distributed with parameter $\tilde{u}(i)$, its square root is approximately normally distributed with expected value of about $\sqrt{\tilde{u}(i)}$ and variance of about 1/4. Under this transformation, the convergence to normality is faster than for the untransformed variable [41]. The most classic VST is the Anscombe transform [42] which has the form $f(u_0) = b\sqrt{u_0 + c}$.

The denoising procedure with the standard variance stabilizing transformation (VST) procedure follows three steps:

(1) apply VST to approximate homoscedasticity;

(2) denoise the transformed data;

(3) apply an inverse VST.

Note that the inverse VST is not just an algebraic inverse of the VST, and must be optimized to avoid bias [43].

Consider any additive signal dependent noisy image, obtained for example by the Gaussian approximation of a Poisson variable explained above. Under this approximation, the noisy image satisfies $\tilde{u} \approx \tilde{u} + g(\tilde{u})n$ where $n \approx \mathcal{N}(0, 1)$. We can search for a function $f$ such that $f(\tilde{u})$ has uniform standard deviation,

$$f(\tilde{u}) \approx f(\tilde{u}) + f'(\tilde{u})g(\tilde{u})n.$$

Forcing the noise term to be constant, $f'(\tilde{u})g(\tilde{u}) = c$, we get

$$f'(\tilde{u}) = \frac{c}{g(\tilde{u})},$$

and integrating

$$f(\tilde{u}) = \int_{0}^{\tilde{u}} \frac{c}{g(t)} dt.$$

When a linear variance noise model is taken, this transformation gives back an Anscombe transform. Most classic denoising algorithms can also be adapted to signal dependent noise. This requires varying the denoising parameters at each pixel, depending on the observed value $\tilde{u}(i)$. Several denoising methods indeed deal directly with the Poisson noise. Wavelet-based denoising methods [44] and [45] propose to adapt the transform threshold to the local noise level of the Poisson process. Lefkimmiatis et al. [46] have explored a Bayesian approach without applying a VST. Deledalle et al., [47] argue that for high noise level it is better to adapt NL-means than to apply a VST. These authors proposed to replace the Euclidean distance between patches by a likelihood estimation taking into account the noise model. This distance can be adapted to each noise model such as the Poisson, the Laplace or the Gamma noise [48], and to more complex (speckle) noise occurring in radar (SAR) imagery [49].

Nonetheless, dealing with a white uniform Gaussian noise makes the discussion on denoising algorithms far easier. The recent papers on the Anscombe transform [43] (for low count Poisson noise) and [50] (for Rician noise) argue that, when combined with suitable forward and inverse VST transformations, algorithms designed for homoscedastic Gaussian noise work just as well as
ad-hoc algorithms signal-dependent noise models. This explains why in the rest of this chapter the noise is assumed uniform, white and Gaussian, having previously applied, if necessary, a VST to the noisy image. This also implies that we deal with *raw* images, namely images as close as possible to the direct camera output before processing. Most reflex cameras, and many compact cameras nowadays give access to this raw image.

But there is definitely a need to denoise current image formats, which have undergone unknown alterations. For example, the JPEG-encoded images given by a camera contain a noise that has been altered by a complex chain of algorithms, ending with lossy compression. Noise in such images cannot be removed by the current state of the art denoising algorithms without a specific adaptation. The key is to have a decent noise model. For this reason, it is important to be able to estimate the noise from the noisy image itself, without assuming any noise model and without relying or trusting any prior information, such as the model of the captor or the kinds of transformations that the image might have undergone.

Compared to the denoising literature, research on noise estimation is a poor cousin. Few papers are dedicated to this topic. Among the recent papers one can mention [51], which argues that images are scale invariant and therefore noise can be estimated by a deviation from this assumption. Unfortunately this method is not easily extended to estimate scale dependent or signal dependent noise, like the one observed in most digital images in compressed format. As a rule of thumb, the noise model is relatively easy to estimate when the raw image comes directly from the imaging system, in which case the noise model is known and only a few parameters must be estimated. For this, efficient methods are described in [7], [40] for Poisson and Gaussian noise.

In this chapter we will focus on methods that allow for local, signal and scale dependent noise. Later on Chapter 3 it will be shown why only considering signal and scale dependent noise is not enough and Chapter 4 will present our proposed algorithm for signal, frequency, and scale dependent noise.

One cannot denoise an image without knowing its noise model. It might be argued that the noise model comes with the knowledge of the imaging device. Nevertheless, the majority of images dealt with by the public or by scientists have lost this information. This loss is caused by format changes of all kinds, which may include resampling, denoising, contrast changes and compression. All of these operations change the noise model and make it *signal and scale dependent*.

The question that arises is why so many researchers are working so hard on denoising models, if their corpus of noisy images is so ill-informed. It is common practice among image processing researchers to add the noise themselves to noise-free images to demonstrate the performance of a method. This proceeding permits to reliably evaluate the denoising performance, based on a controlled ground truth. Nevertheless the denoising performance may, after all, critically depend on how well we are able to estimate the noise. Most world images are actually encoded with lossy JPEG formats. Thus, noise is partly removed by the compression itself. Furthermore, this removal is scale dependent. For example, the JPEG 1985 format divides the image into a disjoint set of 8×8 pixels blocks, computes their DCT, quantizes the coefficients and the small ones are replaced
by zero. This implies that JPEG performs a frequency dependent threshold, equivalent to a basic Wiener filter. The same is true for JPEG 2000 (based on the wavelet transform).

In addition, the Poisson noise of a raw image is signal dependent. The typical image processing operations, demosaicking, white balance and tone curve (contrast change) alter this signal-dependency in a way which depends on the image itself.

In short:

* the noise model is different for each image;

* the noise is signal dependent;

* the noise is scale dependent;

* the knowledge of each dependence is crucial to denoise properly any given image which is not raw, and for which the camera model is available.

Thus, estimating JPEG noise is a complex and risky procedure, as well explained in [52] and [53]. It is argued in [54] that noise can be estimated by involving a denoising algorithm. Again, this procedure is probably too risky for noise and scale dependent signal.

A review and comparison of classic homoscedastic noise estimation methods is presented in Section 3.

## 3. Review of homoscedastic block-based noise estimators

This chapter presents a comparison and discussion of classic block-based noise estimation methods, most of them from [55]. To evaluate them, ten noise-free images are used in the tests. Noiseless images are obtained by taking snapshots with a reflex camera of scenes under good lighting conditions and with a low ISO level. To reduce further the noise level, the average of each block of 5 × 5 pixels is computed, reducing therefore the noise by a 5 factor. Since the images are RGB, the mean of the three channels can be computed, reducing the noise by a further $\sqrt{3}$ factor. The noise is therefore reduced by a $5\sqrt{3} \approx 8.66$ factor, and the images can be considered noise-free. The set of images can be seen in Figure 3. Some of them present large flat regions with little texture, while others do not contain any flat regions and are highly textured. The size of each noise-free test image is 704 × 469 pixels.

In this chapter we only evaluate the performance of the estimator using simulated homoscedastic (fixed variance) noise. The noise is added to the noise-free images and then the estimation given by each noise estimator is compared to the value of the simulated noise. Seven noise levels were applied to these noise-free images: $\sigma \in \{1, 2, 5, 10, 20, 50, 80\}$. Figure 2 shows the result of adding white homoscedastic Gaussian noise with $\sigma \in \{1, 2, 5, 10, 20, 50, 80\}$ to the noise-free image traffic.

These images, from left to right and from top to bottom will be referred to as *bag*, *building1*, *computer*, *dice*, *flowers2*, *hose*, *leaves*, *lawn*, *stairs*, and *traffic*. For the uniform-noise tests, five noise levels were applied to these noise-free images: $\sigma \in \{1, 2, 5, 10, 20, 50, 80\}$. Figure 2 shows the result of adding white Gaussian noise with $\sigma \in \{1, 2, 5, 10, 20, 50, 80\}$ to the noise-free image *traffic*.
<figure><img src="image_15.png" /><figcaption>FIGURE 2. Result of adding white homoscedastic Gaussian noise with $\sigma \in \{2, 5, 10, 20, 50, 80\}$ to the noise-free image traffic. It may need a zoom in to perceive the noise for $\sigma = 2$ and $\sigma = 5$.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>FIGURE 3. Noise-free images used in the tests. From left to right and from top to bottom: *bag*, *building1*, *computer*, *dice*, *flowers2*, *hose*, *leaves*, *lawn*, *stairs* and *traffic*.</figcaption></figure>
<table><caption>FIGURE 4. 3 × 3 window centered at pixel $u_5$.</caption><tr><td>u1</td><td>u2</td><td>u3</td></tr><tr><td>u4</td><td>u5</td><td>u6</td></tr><tr><td>u7</td><td>u8</td><td>u9</td></tr></table>

**3.1. The Average method.** The “Average method” [55] tries to minimize the effect of edges and textures in the image by subtracting a low-pass filtered version of the image from the noisy image: $U = I - f(I)$. In the Average method, function $f$ is the mean of the pixels centered at some given pixel with a $w \times w$ pixels window ($w = 3$). Since the blocks only contain 9 pixels, the sample variance is biased and must be corrected. Indeed, if pixel $u_5$ is in the center of the block (Figure 4), the computed sample variance is $S_9^2(u_5 - \sum_{i=1}^9 u_i) = S_9^2(u_5 - \frac{1}{9}u_5 - \frac{1}{9}\sum_{i=1, i \neq 5}^9 u_i) = S_9^2(\frac{8}{9}u_5) + \frac{1}{81}S_9^2(\sum_{i=1, i \neq 5}^9 u_i) = \frac{64}{81}\sigma^2 + \frac{8}{81}\sigma^2 = \frac{8}{9}\sigma^2$. Here, $S_9^2$ is the sample variance. That is, $S_9^2(B) = \frac{1}{8}\sum_{i=1}^9(B_i - \bar{B})$, where $\bar{B} = \frac{1}{9}\sum_{i=1}^9 B_i$ and $B_i$ are pixel intensities of the block $B$. To avoid that the remaining edges and textures affect the noise estimation, blocks where the magnitude of the intensity gradient is above a certain threshold are discarded. The gradient is computed at the central pixel of each block. For example, if $u_5$ is the central pixel (Figure 4), the magnitude of the gradient would be calculated as $\sqrt{(u_6 - u_5)^2 + (u_8 - u_5)^2}$. The threshold is obtained from the normalized cumulative histogram of the magnitude of the intensity gradients. The threshold is the value in the cumulative histogram corresponding to the $p$ percent of the blocks, with $p = 1\%$. Because of the fact that only blocks with smallest variance are taken into account, the estimate is biased and has to be corrected by a factor learned on noise in order to get the final $\tilde{\sigma}$ estimate. With $p = 1\%$, the empirical correction factor is $k = 1.1609$. The standard deviation of all blocks whose intensity gradient magnitude is below the threshold is computed. The mean of these standard deviations, multiplied by the correction factor, is estimated as $\tilde{\sigma}$. The pseudo-code for the Average method is in Algorithm 2. The results for white Gaussian noise are given in Table 1.
**Algorithm 2** Average noise estimation algorithm.

1: **AVERAGE** - Return the standard deviation of the image noise. **Input** **U**: input image. **Input** **p**: percent of pixels to compute (1% typically). **Output** $\tilde{\sigma}$: estimated standard deviation of the noise in the image.

2: $w = 3$.

3: $k = 1.1608823968593502$.

4: $g[x, y] = \sqrt{(U[x+1, y] - U[x, y])^2 + U[x, y+1] - U[x, y]}, \forall (x, y) \in \text{supp}(I)$. ▷ Compute gradient magnitude

5: $h_c, h_l = \text{HISTOGRAM}(g)$. ▷ Compute the histogram of $g$. Here, $h_c$ is a list containing the number of elements in each bin and $h_l$ a list containing the limits $h_l$ of each bin.

6: $acc = \text{ACC}(h_c, h_l)$. ▷ Accumulated histogram of $g$

7: $i = i : \text{acc}[i] = \frac{p}{100}$. ▷ Index of threshold bin

8: $t = \frac{h_l[i]+h_l[i+1]}{2}$. ▷ Get threshold

9: $B \leftarrow$ all $w \times w$ overlapped blocks in U.

10: $F \leftarrow \text{mean}(B)$. ▷ Mean of all overlapped blocks in U

11: $D = C(B) - F$. ▷ Filter image: the mean of the block is subtracted from the central pixel $C(B)$ in the block $B$.

12: $S = 0$. ▷ List of standard deviations

13: **for** all $w \times w$ blocks $B_D \in D$, where $(x, y)$ is the center of the block **do**

14: **if** $g[x, y] < t$ **then**

15: $S \leftarrow \sqrt[9]{8S_{w^2-1}^2(B_D)}$

16: **end if**

17: **end for**

return $\tilde{\sigma} = k \times \text{mean}(S)$. ▷ $k$ is the correction factor

<table><caption>TABLE 1. Average method results with simulated homoscedastic white Gaussian noise.</caption><thead><tr><th>Image / &sigma;&#771;</th><th>&sigma; = 1</th><th>&sigma; = 2</th><th>&sigma; = 5</th><th>&sigma; = 10</th><th>&sigma; = 20</th><th>&sigma; = 50</th><th>&sigma; = 80</th></tr></thead><tbody><tr><td>bag</td><td>6.00</td><td>6.75</td><td>9.18</td><td>13.56</td><td>22.45</td><td>51.51</td><td>80.79</td></tr><tr><td>building1</td><td>1.99</td><td>3.25</td><td>6.77</td><td>12.23</td><td>22.12</td><td>51.62</td><td>81.17</td></tr><tr><td>computer</td><td>3.51</td><td>4.40</td><td>7.68</td><td>12.79</td><td>22.58</td><td>51.96</td><td>81.56</td></tr><tr><td>dice</td><td>1.30</td><td>2.27</td><td>5.26</td><td>10.27</td><td>20.27</td><td>50.35</td><td>80.61</td></tr><tr><td>flowers2</td><td>1.69</td><td>2.71</td><td>5.74</td><td>10.63</td><td>20.57</td><td>50.67</td><td>80.37</td></tr><tr><td>hose</td><td>2.46</td><td>3.31</td><td>6.05</td><td>10.80</td><td>20.60</td><td>50.25</td><td>80.37</td></tr><tr><td>leaves</td><td>4.91</td><td>5.58</td><td>7.98</td><td>12.28</td><td>21.58</td><td>50.89</td><td>80.48</td></tr><tr><td>lawn</td><td>5.13</td><td>5.79</td><td>8.05</td><td>12.30</td><td>21.58</td><td>50.96</td><td>80.78</td></tr><tr><td>stairs</td><td>2.37</td><td>3.33</td><td>6.08</td><td>10.75</td><td>20.48</td><td>50.39</td><td>80.39</td></tr><tr><td>traffic</td><td>4.04</td><td>5.26</td><td>8.63</td><td>13.74</td><td>23.37</td><td>52.30</td><td>81.44</td></tr><tr><td>constant</td><td>1.00</td><td>2.01</td><td>5.00</td><td>10.01</td><td>20.04</td><td>50.06</td><td>80.52</td></tr></tbody></table>
<table><caption>TABLE 2. Median method results with simulated homoscedastic white Gaussian noise.</caption><thead><tr><th>Image / ̃σ</th><th>σ = 1</th><th>σ = 2</th><th>σ = 5</th><th>σ = 10</th><th>σ = 20</th><th>σ = 50</th><th>σ = 80</th></tr></thead><tbody><tr><td>bag</td><td>6.39</td><td>7.18</td><td>9.67</td><td>14.03</td><td>22.92</td><td>52.01</td><td>81.46</td></tr><tr><td>building1</td><td>2.03</td><td>3.33</td><td>6.92</td><td>12.51</td><td>22.46</td><td>52.11</td><td>81.80</td></tr><tr><td>computer</td><td>3.63</td><td>4.56</td><td>8.00</td><td>13.22</td><td>23.06</td><td>52.61</td><td>82.34</td></tr><tr><td>dice</td><td>1.33</td><td>2.31</td><td>5.34</td><td>10.38</td><td>20.42</td><td>50.68</td><td>81.23</td></tr><tr><td>flowers2</td><td>1.78</td><td>2.85</td><td>5.97</td><td>10.90</td><td>20.87</td><td>51.12</td><td>80.89</td></tr><tr><td>hose</td><td>2.56</td><td>3.42</td><td>6.18</td><td>10.94</td><td>20.82</td><td>50.68</td><td>81.00</td></tr><tr><td>leaves</td><td>5.13</td><td>5.83</td><td>8.23</td><td>12.53</td><td>21.86</td><td>51.34</td><td>81.10</td></tr><tr><td>lawn</td><td>5.37</td><td>6.04</td><td>8.28</td><td>12.55</td><td>21.90</td><td>51.39</td><td>81.53</td></tr><tr><td>stairs</td><td>2.46</td><td>3.45</td><td>6.22</td><td>10.90</td><td>20.70</td><td>50.79</td><td>81.05</td></tr><tr><td>traffic</td><td>4.21</td><td>5.58</td><td>9.14</td><td>14.26</td><td>23.98</td><td>53.04</td><td>82.22</td></tr><tr><td>constant</td><td>1.01</td><td>2.02</td><td>5.03</td><td>10.07</td><td>20.17</td><td>50.36</td><td>81.04</td></tr></tbody></table>

**3.2. Median method.** The Median method [55] is a variant of the Average (Section 3.1). In this method, function *f* to filter the image is the median. The same correction factor used for the Average method is used here. The white Gaussian noise estimates are given in Table 2. The only difference with the Average method in the algorithm description is that the mean of the block is replaced by the median of the block.
<table><caption>TABLE 3. Block method results with simulated homoscedastic white Gaussian noise.</caption><thead><tr><th>Image / ̂σ</th><th>σ = 1</th><th>σ = 2</th><th>σ = 5</th><th>σ = 10</th><th>σ = 20</th><th>σ = 50</th><th>σ = 80</th></tr></thead><tbody><tr><td>bag</td><td>3.51</td><td>4.20</td><td>7.07</td><td>12.36</td><td>22.10</td><td>51.87</td><td>81.00</td></tr><tr><td>building1</td><td>1.23</td><td>2.20</td><td>5.28</td><td>10.50</td><td>20.58</td><td>51.04</td><td>81.02</td></tr><tr><td>computer</td><td>1.42</td><td>2.42</td><td>5.54</td><td>10.62</td><td>20.82</td><td>51.24</td><td>81.63</td></tr><tr><td>dice</td><td>1.20</td><td>2.16</td><td>5.11</td><td>10.09</td><td>20.13</td><td>50.21</td><td>80.04</td></tr><tr><td>flowers2</td><td>1.23</td><td>2.23</td><td>5.34</td><td>10.47</td><td>20.53</td><td>50.52</td><td>80.26</td></tr><tr><td>hose</td><td>1.61</td><td>2.59</td><td>5.61</td><td>10.58</td><td>20.45</td><td>50.40</td><td>79.97</td></tr><tr><td>leaves</td><td>2.96</td><td>3.65</td><td>6.54</td><td>11.53</td><td>21.33</td><td>50.75</td><td>80.81</td></tr><tr><td>lawn</td><td>2.43</td><td>3.29</td><td>6.35</td><td>11.48</td><td>21.14</td><td>50.55</td><td>80.75</td></tr><tr><td>stairs</td><td>1.37</td><td>2.44</td><td>5.51</td><td>10.59</td><td>20.39</td><td>50.31</td><td>79.76</td></tr><tr><td>traffic</td><td>1.32</td><td>2.33</td><td>5.55</td><td>10.77</td><td>21.23</td><td>51.34</td><td>81.31</td></tr><tr><td>constant</td><td>1.00</td><td>1.99</td><td>5.01</td><td>9.93</td><td>20.00</td><td>49.90</td><td>80.11</td></tr></tbody></table>

**3.3. Block method.** The Block method [56, 57] measures the noise from the set of 7 × 7 pixels blocks with minimal variance. The size of this set of blocks is chosen as a fixed small percentage $p$ of the total number of pixels in the image. The direct estimation $\hat{\sigma}^2$ is the mean of the set of sample variances under the $p$-percentile. Of course, this produces a biased result, because only the blocks with the minimal variance are chosen. Therefore, to obtain $\tilde{\sigma}^2$, the direct sampled variance estimation $\hat{\sigma}^2$ has to be multiplied by a correction factor. This factor is $k = 1.8350$ for $p = 1\%$. The pseudo-code for the Block method is in Algorithm 3 and its results with white Gaussian noise in Table 3. Two different block sized where tests: $w = 11$ (Table 4) and $w = 15$ (Table 5), but increasing the block size does not give, in general, better results. The computation time, however, increases significantly with the block size. The correction factors for $w = 11$ is $k = 1.4429$ and for $w = 15$ it is $k = 1.2994$.

**Algorithm 3 Block algorithm.**

1: **BLOCK** - Return the standard deviation of the image noise. **Input U**: input image. **Input p**: percentage of blocks to compute (1% typically). **Output $\tilde{\sigma}$**: estimated standard deviation of the noise in the image.

2: $w = 7$.

3: $k = 1.8350297625014702$.

4: $B \leftarrow$ all $w \times w$ blocks in U.

5: $V = 0$ ▷ List of sample variances

6: **for** all $w \times w$ blocks $B_D \in B$ **do**

7: $V \leftarrow S_{w^2-1}^2(B_D)$. ▷ Add sample variances to list $V$

8: **end for**

9: $\tilde{\sigma} = k\sqrt{\text{Mean}(V[0...|V|_{100}])}$. ▷ $k$ is the correction factor
<table><caption>TABLE 4. Block method results with simulated homoscedastic white Gaussian noise (w = 11).</caption><thead><tr><th>Image / ̃σ</th><th>σ = 1</th><th>σ = 2</th><th>σ = 5</th><th>σ = 10</th><th>σ = 20</th><th>σ = 50</th><th>σ = 80</th></tr></thead><tbody><tr><td>bag</td><td>4.58</td><td>5.01</td><td>7.31</td><td>12.46</td><td>22.19</td><td>52.03</td><td>81.33</td></tr><tr><td>building1</td><td>1.23</td><td>2.17</td><td>5.22</td><td>10.35</td><td>20.43</td><td>50.94</td><td>80.86</td></tr><tr><td>computer</td><td>1.42</td><td>2.41</td><td>5.48</td><td>10.59</td><td>20.80</td><td>50.94</td><td>81.62</td></tr><tr><td>dice</td><td>1.21</td><td>2.18</td><td>5.12</td><td>10.12</td><td>20.14</td><td>50.22</td><td>80.20</td></tr><tr><td>flowers2</td><td>1.20</td><td>2.19</td><td>5.29</td><td>10.42</td><td>20.60</td><td>50.59</td><td>80.35</td></tr><tr><td>hose</td><td>1.70</td><td>2.60</td><td>5.63</td><td>10.57</td><td>20.54</td><td>50.46</td><td>80.23</td></tr><tr><td>leaves</td><td>3.37</td><td>3.91</td><td>6.54</td><td>11.50</td><td>21.33</td><td>50.97</td><td>80.88</td></tr><tr><td>lawn</td><td>2.77</td><td>3.46</td><td>6.25</td><td>11.36</td><td>21.13</td><td>50.78</td><td>80.72</td></tr><tr><td>stairs</td><td>1.31</td><td>2.37</td><td>5.43</td><td>10.54</td><td>20.43</td><td>50.47</td><td>79.82</td></tr><tr><td>traffic</td><td>1.31</td><td>2.28</td><td>5.42</td><td>10.61</td><td>21.05</td><td>51.37</td><td>81.17</td></tr><tr><td>constant</td><td>1.00</td><td>2.00</td><td>5.01</td><td>9.93</td><td>20.01</td><td>49.96</td><td>80.03</td></tr></tbody></table>

<table><caption>TABLE 5. Block method results with simulated homoscedastic white Gaussian noise (w = 15).</caption><thead><tr><th>Image / ̃σ</th><th>σ = 1</th><th>σ = 2</th><th>σ = 5</th><th>σ = 10</th><th>σ = 20</th><th>σ = 50</th><th>σ = 80</th></tr></thead><tbody><tr><td>bag</td><td>5.67</td><td>5.98</td><td>7.88</td><td>12.65</td><td>22.30</td><td>52.20</td><td>81.43</td></tr><tr><td>building1</td><td>1.24</td><td>2.16</td><td>5.20</td><td>10.33</td><td>20.36</td><td>50.90</td><td>81.13</td></tr><tr><td>computer</td><td>1.43</td><td>2.40</td><td>5.45</td><td>10.52</td><td>20.70</td><td>50.95</td><td>81.40</td></tr><tr><td>dice</td><td>1.22</td><td>2.19</td><td>5.14</td><td>10.14</td><td>20.12</td><td>50.21</td><td>80.23</td></tr><tr><td>flowers2</td><td>1.19</td><td>2.18</td><td>5.26</td><td>10.34</td><td>20.55</td><td>50.66</td><td>80.50</td></tr><tr><td>hose</td><td>1.79</td><td>2.64</td><td>5.61</td><td>10.56</td><td>20.59</td><td>50.35</td><td>80.21</td></tr><tr><td>leaves</td><td>3.79</td><td>4.26</td><td>6.66</td><td>11.46</td><td>21.23</td><td>50.99</td><td>80.90</td></tr><tr><td>lawn</td><td>2.99</td><td>3.58</td><td>6.20</td><td>11.26</td><td>21.18</td><td>50.83</td><td>80.52</td></tr><tr><td>stairs</td><td>1.28</td><td>2.33</td><td>5.38</td><td>10.53</td><td>20.48</td><td>50.50</td><td>79.91</td></tr><tr><td>traffic</td><td>1.34</td><td>2.28</td><td>5.36</td><td>10.54</td><td>20.94</td><td>51.32</td><td>81.27</td></tr><tr><td>constant</td><td>1.00</td><td>2.00</td><td>5.01</td><td>9.97</td><td>19.97</td><td>49.99</td><td>80.08</td></tr></tbody></table>
**3.4. Gradient method.** In the Gradient method [58, 59], the magnitude of the intensity gradient is computed using a 3 × 3 least squares fit for all possible pixels of the image (excluding those pixels in the edges of the image, where the gradient can not be computed). For the 3 × 3 least squares fit, both the Sobel and Prewitt operators were tested. The Sobel operator uses the following matrices for the horizontal and vertical gradient components:

$$S_x = \begin{pmatrix} -1 & 0 & +1 \\ -2 & 0 & +2 \\ -1 & 0 & +1 \end{pmatrix}, S_y = \begin{pmatrix} -1 & -2 & -1 \\ 0 & 0 & 0 \\ +1 & +2 & +1 \end{pmatrix}$$

The Prewitt operator uses these matrices:

$$P_x = \begin{pmatrix} -1 & 0 & +1 \\ -2 & 0 & +2 \\ -1 & 0 & +1 \end{pmatrix}, P_y = \begin{pmatrix} -1 & -2 & -1 \\ 0 & 0 & 0 \\ +1 & +2 & +1 \end{pmatrix}$$

The same results were obtained computing the gradient with the Sobel or the Prewitt operators. Because for a white Gaussian noise the horizontal and vertical magnitudes in the intensity gradient are normally distributed (with mean $\mu = 0$ and variance $\sigma^2$), the magnitude of the noise intensity gradient is Rayleigh distributed. Indeed, if $X$ and $Y$ are normal random variables with mean $\mu = 0$ and variance $\sigma_R^2$, then $R = \sqrt{X^2 + Y^2}$ is Rayleigh distributed, with parameter $\sigma_R$. Note that $\sigma_R$ is not the variance of the noise, but the variance of the directional gradient components.

It is important that blocks do not overlap, because the random variables $X$ and $Y$ must be independent in order to obtain a Rayleigh distribution $R$ with $\sigma$ parameter. A histogram of the magnitude of the intensity gradient is made. As explained in Section 1 of Chapter 11, edges and textures affect the right part of the histogram mainly. If the image is not dominated by edges and textures, the $\sigma_R$ parameter can be estimated by finding the mode (the value at which the distribution is maximal) of the Rayleigh distribution. However, this is only possible if edges and textures are moderate, and therefore not affecting the part of the histogram containing the mode. Before finding the mode of the distribution, the histogram is filtered at least four times. To low-pass filter the histogram, the same procedure used to filter the noise curves explained in Chapter 1.2 was applied.

Once the mode $\sigma_R$ is found in the filtered histogram, the noise standard deviation $\tilde{\sigma}$ can be obtained by multiplying the mode by a correction factor $f = 0.7607$. Unfortunately, trying to fit a Rayleigh distribution to the normalized histogram is a difficult task, because edges and textures affect the histogram, changing the shape of the distribution. When normalizing the histogram, these undesired values in the tail change the shape of the normalized distribution. Therefore, even if a Rayleigh if fitted onto the steep rising portion of the affected histogram, it would not be fitting the noise image distribution. Figure 5 illustrates this problem. In green, the ground-truth distribution (pure noise). In blue, the distribution obtained with a noisy image with several edges and textures. In red, the Rayleigh distribution fitted to the steep rising part (20%) of the histogram. On the left, fitting a Rayleigh distribution to the steep rising portion histogram of the
<figure><img src="image_17.png" /><figcaption>FIGURE 5. Left: fitting a Rayleigh distribution to the steep rising portion (20%) histogram of the variance of pure-noise ($\sigma = 5$) patches. The fit (red) is very close to the ground truth (green) and therefore the obtained $\hat{\sigma}$ is accurate. Right: trying to fit a Rayleigh distribution to the steep rising portion (20%) histogram of the variance of the *leaves* image (many edges and textures). The obtained parameter $\sigma$ does not correspond to that of the noise (ground truth distribution).</figcaption></figure>

variance of pure-noise ($\sigma = 5$) patches. The fit (red) is computed over the ground truth (green). Thus the obtained $\hat{\sigma}$ is accurate. On the right, trying to fit a Rayleigh distribution to the steep rising portion histogram of the variance of the *leaves* image (many edges and textures). The obtained parameter $\sigma$ does not correspond to that of the noise (ground truth distribution). The pseudo-code for the algorithm is presented in Algorithm 4 and the results for white Gaussian noise in Table 6.

**Algorithm 4** Gradient noise estimation algorithm.

1: **GRADIENT** - Return the standard deviation of the image noise. **Input U**: input image. **Output $\tilde{\sigma}$**: estimated standard deviation of the noise in the image.

2: $f = 0.76071735110441019$.

▷ Correction factor

3: **for** $[x, y] \in [1, N_x - 2] \times [1, N_y - 2]$ **do**

4: $g_x[x, y] = (U[x+1, y-1]+2U[x+1, y]+U[x+1, y+1]-U[x-1, y-1]-2U[x-1, y]-U[x-1, y+1])/8$

5: $g_y[x, y] = (U[x-1, y-1]+2U[x, y-1]+U[x+1, y-1]-U[x-1, y+1]-2U[x, y+1]-U[x+1, y+1])/8$.

6: $g[x, y] = \sqrt{g_x^2 + g_y^2}$

▷ Compute $3 \times 3$ least squares fit gradient magnitude

7: **end for**

8: $H = \text{HISTOGRAM}(g)$.

▷ Compute the histogram of $g$. Here, $H$ is a list with the number of

9: $k = \operatorname{argmax}_k(H[k])$.

10: $\text{mode} = \text{median}(\text{bin}_k)$.

▷ Median of the elements in bin $\#k$

11: $\tilde{\sigma} = f \times \text{mode}$.
<table><caption>TABLE 6. Gradient method results with simulated homoscedastic white Gaussian noise.</caption><thead><tr><th>Image / σ̃</th><th>σ = 1</th><th>σ = 2</th><th>σ = 5</th><th>σ = 10</th><th>σ = 20</th><th>σ = 50</th><th>σ = 80</th></tr></thead><tbody><tr><td>bag</td><td>5.04</td><td>5.42</td><td>7.47</td><td>12.58</td><td>20.64</td><td>47.99</td><td>74.67</td></tr><tr><td>building1</td><td>0.76</td><td>2.78</td><td>6.20</td><td>11.14</td><td>20.91</td><td>51.57</td><td>77.98</td></tr><tr><td>computer</td><td>1.02</td><td>1.09</td><td>6.75</td><td>11.67</td><td>20.52</td><td>49.16</td><td>80.61</td></tr><tr><td>dice</td><td>0.92</td><td>2.10</td><td>4.97</td><td>9.71</td><td>19.94</td><td>51.17</td><td>80.50</td></tr><tr><td>flowers2</td><td>0.74</td><td>2.73</td><td>6.35</td><td>9.83</td><td>19.43</td><td>51.44</td><td>83.17</td></tr><tr><td>hose</td><td>2.14</td><td>3.12</td><td>5.81</td><td>10.68</td><td>20.39</td><td>48.95</td><td>75.72</td></tr><tr><td>leaves</td><td>4.26</td><td>5.35</td><td>6.77</td><td>12.09</td><td>19.19</td><td>51.54</td><td>82.97</td></tr><tr><td>lawn</td><td>4.06</td><td>5.05</td><td>7.07</td><td>11.81</td><td>21.49</td><td>49.06</td><td>76.83</td></tr><tr><td>stairs</td><td>2.21</td><td>3.19</td><td>5.81</td><td>10.52</td><td>20.37</td><td>50.67</td><td>83.89</td></tr><tr><td>traffic</td><td>0.96</td><td>1.06</td><td>6.62</td><td>11.33</td><td>23.15</td><td>49.29</td><td>80.19</td></tr><tr><td>constant</td><td>0.98</td><td>2.00</td><td>5.31</td><td>9.43</td><td>17.72</td><td>50.89</td><td>80.32</td></tr></tbody></table>
<table><caption>TABLE 7. Summary of the Pyramid method rules.</caption><thead><tr><th>Rule #</th><th>Condition</th><th>Computation of ̃σ</th></tr></thead><tbody><tr><td>1</td><td>q[l, 4] &lt; 1</td><td>̃σ = 0</td></tr><tr><td>2</td><td>l<sub>0</sub> = n</td><td>̃σ = &radic;v[n]</td></tr><tr><td>3, 4</td><td>l<sub>0</sub> = 3 or l<sub>0</sub> = 4</td><td>WARNING: CAN NOT ESTIMATE</td></tr><tr><td>5.1</td><td>ρ<sub>0</sub> &lt; ρ ≤ ρ<sub>1</sub></td><td>̃σ = &radic;v[3]</td></tr><tr><td>5.2</td><td>ρ<sub>1</sub> &lt; ρ ≤ ρ<sub>2</sub></td><td>̃σ = &radic;δv[3] + (1 − δ)v[4]</td></tr><tr><td>5.3</td><td>ρ<sub>2</sub> &lt; ρ ≤ ρ<sub>3</sub></td><td>̃σ = &radic;δv[4] + (1 − δ)v[5]</td></tr><tr><td>5.4</td><td>ρ<sub>3</sub> &lt; ρ ≤ ρ<sub>4</sub></td><td>̃σ = &radic;δv[5] + (1 − δ)v[6]</td></tr><tr><td>6.1</td><td>T &lt; α[l<sub>0</sub> − 1] ≤ 0</td><td>̃σ = &radic;δv[l<sub>0</sub> − 2] + (1 − δ)v[l<sub>0</sub> − 1]</td></tr><tr><td>6.2</td><td>−1 &lt; α[l<sub>0</sub>] ≤ −0.5</td><td>̃σ = &radic;δv[l<sub>0</sub> − 2] + (1 − δ)v[l<sub>0</sub> − 1]</td></tr><tr><td>6.3</td><td>−0.5 &lt; α[l<sub>0</sub>] ≤ T</td><td>̃σ = &radic;((1 + δ)v[l<sub>0</sub> + 1] + (1 − δ)v[l<sub>0</sub>])/2</td></tr></tbody></table>

**3.5. Pyramid method.** The Pyramid method [6] considers $n$ levels in the input image, $l = 2, 3, \dots, n$, where $N = 2^n$ is the size of the image. At each level, the variance of all non-overlapping blocks of size $2^l \times 2^l$ is computed. To compute this variance, the unbiased estimator for the population variance is used: $S_{B-1}^2(X) = \frac{1}{B-1} \sum_{i=1}^B (X - \bar{X})$, where $\bar{X} = \frac{1}{B} \sum_{i=1}^B x_i$ and $B$ the number of elements in $X$. For $l = 2, 3, \dots, n-1$, the four smallest values for each level $l$ are stored. A slippage test in levels $l = 2, 3, \dots, n-2$ is done and the set $v(l)$ is defined as the mean of the smallest values (between one and four values are used, depending on the case). Finally, the contribution to $v(l)$ (signal) and that of the noise is obtained by analyzing the shape of $\alpha(l) = \frac{v(l-1)}{v(l)} - \beta(l)$. The function $\beta$ is experimentally determined by the authors to be $\beta(l) = 1 - \frac{2^{6-l}}{10}$. The method proposes several conditional rules that, if held, give an estimate of the noise. It can happen that the method determines that the noise is poorly separated from the signal and no result is obtained. However, the method can detect this situation and prevents giving a bad result. The results for the Pyramid method for white Gaussian noise are given in Table 9. Code NE means that the method considered that noise and signal were not well separated and therefore it did not give any estimation. Rule 6.3 was applied when no rule held.

Function $\alpha$ is computed at each level $l$ as $\alpha[l] = \frac{v[l-1]}{v[l]} - \beta[l]$. $l_u$ is the minimal argument $l$ that gets $\alpha[l] < 0$ for $l = 3, \dots, n$ and $l_0$ the minimal argument $l$ that gets $\sum_{i=l_u}^l \alpha[l] < T$ with $l = 3, \dots, n-2$, where $T = -0.1$. The pyramid method proposes several rules as conditions that if held, provide a way to get the noise standard deviation estimation. This rules are summarized in Table 7. Some of the rules interpolate to get the value of $\delta$. The interpolation function is: INTERPOLATION($a, \rho, b$) = $\frac{\rho-b}{a-b}$. Rules 5.2, 5.3, 5.4, 6.1, 6.2 and 6.3 use interpolation. The interpolation parameters to get $\delta$ according to the rule are given in Table 8. Rules 5.1, 5.2, 5.3 and 5.4 use the following variables: $\rho = \alpha[5] + \alpha[6]$, $\rho_0 = -2$, $\rho_1 = -1.5$, $\rho_2 = -1$, $\rho_3 = -0.5$, $\rho_4 = T$, where $T = -0.1$.
**Algorithm 5** Pyramid noise estimation algorithm.

1: **PYRAMID** - Return the standard deviation of the image noise. **Input U**: input image. **Output $\tilde{\sigma}$**: estimated standard deviation of the noise in the image.

2: $n = \frac{\log(N)}{\log(2)}$. ▷ Image is $N \times N$ and $N = 2^n$.

3: $\beta[l] = 1 - 0.1 \times 2^{6-l}$ for $l = 1...n$.

4: $Q = \emptyset$. i-order statistic $q[l, i]$, $i \in \{1, 2, 3, 4\}$

5: $r = \emptyset$. ▷ Thresholds $r[l, i]$, $l \in \{1, 2, ..., n\}$, $i \in \{1, 2, 3\}$

6: $v = \emptyset$. ▷ Variances $v[l]$, $l \in \{1, 2, ..., n\}$

7: **for** $l = 1$ to $n$ **do**

8: $K = 4^{n-l}$, $c = 2^l$, $V = \emptyset$. ▷ List of variances in level $l$

9: **for** $x = 0$ to $N$ step $c$ **do**

10: **for** $y = 0$ to $N$ step $c$ **do**

11: $v = \text{Var}(u[x..., x + c - 1, y..., y + c - 1])$. ▷ Sample variances

12: $V \leftarrow v$. ▷ Store sample variance

13: **end for**

14: **end for**

15: $q[l, 1], q[l, 2], q[l, 3], q[l, 4] = \text{Min}(V)$. ▷ Find 1, 2, 3 and 4 order statistics of $V$

16: $r[l, 1] = (q[l, 2] - q[l, 1])/(q[l, 4] - q[l, 1])$

17: $r[l, 2] = (q[l, 3] - q[l, 2])/(q[l, 4] - q[l, 2])$

18: $r[l, 3] = (q[l, 4] - q[l, 3])/(q[l, 4] - q[l, 2])$

19: $v[l] = \text{Mean}(q[l, :])$ if $r[l, 1] \le 0.5$.

20: $v[l] = (q[l, 2] + q[l, 3] + q[l, 4])/3$ if $r[l, 2] \le 0.7$.

21: $v[l] = (q[l, 3] + q[l, 4])/2$ if $r[l, 3] \le 0.7$.

22: $v[l] = q[l, 4]$ otherwise.

23: **end for**

24: Compute $\alpha$ and $l_u$

25: **PYRAMID_APPLY_RULES**.

<table><caption>TABLE 8. Pyramid method interpolation parameters.</caption><thead><tr><th>Rule #</th><th>Interpolation</th></tr></thead><tbody><tr><td>5.2</td><td>$\delta = \text{INTERPOLATION}(\rho_1, \rho, \rho_2)$</td></tr><tr><td>5.3</td><td>$\delta = \text{INTERPOLATION}(\rho_2, \rho, \rho_3)$</td></tr><tr><td>5.4</td><td>$\delta = \text{INTERPOLATION}(\rho_3, \rho, \rho_4)$</td></tr><tr><td>6.1</td><td>$\delta = \text{INTERPOLATION}(T, \alpha[l_0 - 1], 0)$</td></tr><tr><td>6.2</td><td>$\delta = \text{INTERPOLATION}(-1, \alpha[l_0], -0.5)$</td></tr><tr><td>6.3</td><td>$\delta = \text{INTERPOLATION}(-0.5, \alpha[l_0], T)$</td></tr></tbody></table>
<table><caption>TABLE 9. Pyramid method results with simulated homoscedastic white Gaussian noise.</caption><thead><tr><th>Image / ̃σ</th><th>σ = 1</th><th>σ = 2</th><th>σ = 5</th><th>σ = 10</th><th>σ = 20</th><th>σ = 50</th><th>σ = 80</th></tr></thead><tbody><tr><td>bag</td><td>NE</td><td>NE</td><td>NE</td><td>19.40</td><td>26.37</td><td>49.03</td><td>73.02</td></tr><tr><td>building1</td><td>0.00</td><td>1.96</td><td>4.44</td><td>8.92</td><td>21.12</td><td>62.97</td><td>83.91</td></tr><tr><td>computer</td><td>0.00</td><td>1.82</td><td>4.38</td><td>10.08</td><td>20.63</td><td>69.31</td><td>75.54</td></tr><tr><td>dice</td><td>0.00</td><td>0.00</td><td>4.32</td><td>10.88</td><td>29.48</td><td>43.42</td><td>75.34</td></tr><tr><td>flowers2</td><td>0.00</td><td>1.58</td><td>5.28</td><td>9.80</td><td>23.04</td><td>64.19</td><td>84.65</td></tr><tr><td>hose</td><td>0.00</td><td>NE</td><td>6.30</td><td>12.37</td><td>19.46</td><td>46.34</td><td>79.86</td></tr><tr><td>leaves</td><td>3.25</td><td>3.63</td><td>6.72</td><td>13.28</td><td>28.53</td><td>50.15</td><td>83.29</td></tr><tr><td>lawn</td><td>3.13</td><td>4.39</td><td>6.66</td><td>10.18</td><td>27.44</td><td>49.68</td><td>74.18</td></tr><tr><td>stairs</td><td>0.00</td><td>1.87</td><td>5.37</td><td>13.28</td><td>19.25</td><td>48.51</td><td>79.78</td></tr><tr><td>traffic</td><td>0.00</td><td>NE</td><td>5.01</td><td>10.37</td><td>25.69</td><td>52.96</td><td>67.94</td></tr><tr><td>constant</td><td>0.00</td><td>0.00</td><td>4.61</td><td>7.27</td><td>19.39</td><td>36.29</td><td>77.57</td></tr></tbody></table>
<figure><img src="image_18.png" /><figcaption>FIGURE 6. Problem of the *Scatter* method when trying to figure out the standard deviation of the noise from the peak of the histogram in images dominated by textures and edges. On the left, the normalized histograms of standard deviations for pure noise and *bag* images. On the right, the normalized standard deviation histograms for pure noise and *dice* images. The noise standard deviation is $\sigma = 5$ for both images.</figcaption></figure>

**3.6. The Scatter method.** The Scatter method [60] is intended not only for additive noise, but also for multiplicative noise. It uses the Hough transform to have a plot of the averages $M$ and variances $V$ of all 8 × 8 blocks in the image. Then, it fits a straight line to the main cluster, using the Hough transform. The slope and intercept of this line are the multiplicative and additive components of the noise. In the case of signal independent white Gaussian noise, the Hough transform boils down to the largest peak of the standard deviations histogram. Therefore the standard deviation is figured out by the largest peak, which is enough to compare this method with the other ones.

In the white Gaussian noise test, it was used the largest peak in the histogram to figure out $\tilde{\sigma}$ (Algorithm 7). When the peak of the histogram technique is used, $\hat{\sigma}$ must be corrected by a factor $k = 1.0229$ to get the final estimation $\tilde{\sigma}$. If the image is dominated by edges and textures, the Scatter method fails to get a reliable estimation because the large peak can be produced by them and not by the noise.

Figure 6 illustrates this problem. On the left, it shows the normalized histogram of standard deviations corresponding to pure noise and the one corresponding to the *bag* image, both with $\sigma = 5$. On the right, the normalized histograms of a pure noise and of the *dice* image. The image *bag* has many edges and textures, while the *dice* image has not. The estimate in *bag* is wrong, because the variance retrieved from the blocks corresponds to the variations of the signal, rather than noise. The estimate in *dice* is accurate, because this image contains large flat regions and little texture (although a small bump can be observed on the right part of the histogram). When the peak of the histogram is used to get $\hat{\sigma}$, a correction factor $k = 1.0229$ must be used to get $\tilde{\sigma}$. The pseudo-code for the Scatter method for signal-dependent noise is Algorithm 6. However, here
the method is evaluated for white Gaussian noise (Algorithm 7) and its results are shown in Table 10.

**Algorithm 6** Scatter signal-dependent noise estimation algorithm.

1: **SCATTER** - Return the standard deviation of the image noise. **Input U:** input image. **Output $\tilde{\sigma}$:** estimated standard deviation of the noise in the image.

2: $w = 8$.

3: $B \leftarrow$ all $w \times w$ block in $U$.

4: $M \leftarrow \text{Mean}(B)$. ▷ Obtain the mean of each $w \times w$ block.

5: $V \leftarrow \text{Var}(B)$. ▷ Obtain the variance of each $w \times w$ block.

6: Construct scatter plot of $V$ vs. $M^2$.

7: Fit a straight line fitted to the main cluster, using the Hough transform.

8: Model the fitted line as $Y = mU + a$. ▷ $m$ is the multiplicative and $a$ the additive noise component.

**Algorithm 7** Scatter white Gaussian noise estimation algorithm.

1: **SCATTER\_PEAK** - Return the standard deviation of the image noise. **Input U:** input image. **Output $\tilde{\sigma}$:** estimated standard deviation of the noise in the image.

2: $w = 8$.

3: $k = 1.0229490535653729$.

4: $B \leftarrow$ all $w \times w$ block in $U$.

5: $M \leftarrow \text{Mean}(B)$. ▷ Obtain the mean of each $w \times w$ block

6: $V \leftarrow \text{Var}(B)$. ▷ Obtain the variance of each $w \times w$ block

7: $H = \text{HISTOGRAM}(V)$. ▷ Compute the histogram of $V$. Here, $H$ is a list with the number of elements in each bin.

8: $m = \text{argmax}_i(H[i])$. ▷ Look for the bin with most elements (mode of the histogram).

9: $\hat{\sigma} = \text{Mean}[H(m)]$.

return $\tilde{\sigma} = k \times \hat{\sigma}$. ▷ Corrected value

These experiments show clearly that block variances are not at all reliable to estimate low noises. See how all images having textures lead to a strong over estimation of the noise. They therefore amply justify the use high pass filters to clean up the blocks from slow variations due to the signal and give prominence to noise components.
<table><caption>TABLE 10. Scatter method results with simulated homoscedastic white Gaussian noise.</caption><thead><tr><th>Image / σ̃</th><th>σ = 1</th><th>σ = 2</th><th>σ = 5</th><th>σ = 10</th><th>σ = 20</th><th>σ = 50</th><th>σ = 80</th></tr></thead><tbody><tr><td>bag</td><td>13.52</td><td>13.50</td><td>14.43</td><td>16.87</td><td>23.75</td><td>52.54</td><td>82.71</td></tr><tr><td>building1</td><td>1.34</td><td>2.05</td><td>5.32</td><td>10.47</td><td>20.53</td><td>51.57</td><td>81.38</td></tr><tr><td>computer</td><td>1.54</td><td>2.32</td><td>5.31</td><td>10.65</td><td>20.63</td><td>51.58</td><td>83.29</td></tr><tr><td>dice</td><td>1.32</td><td>2.20</td><td>5.18</td><td>10.19</td><td>20.38</td><td>50.86</td><td>80.96</td></tr><tr><td>flowers2</td><td>1.29</td><td>2.30</td><td>5.26</td><td>10.51</td><td>20.53</td><td>51.18</td><td>80.14</td></tr><tr><td>hose</td><td>2.88</td><td>2.98</td><td>5.63</td><td>10.84</td><td>21.15</td><td>51.03</td><td>80.12</td></tr><tr><td>leaves</td><td>6.85</td><td>6.59</td><td>8.01</td><td>12.05</td><td>21.59</td><td>51.67</td><td>82.07</td></tr><tr><td>lawn</td><td>6.52</td><td>6.71</td><td>8.07</td><td>11.96</td><td>21.47</td><td>51.42</td><td>80.72</td></tr><tr><td>stairs</td><td>1.62</td><td>2.49</td><td>5.61</td><td>10.68</td><td>20.76</td><td>50.71</td><td>80.56</td></tr><tr><td>traffic</td><td>1.18</td><td>2.20</td><td>5.28</td><td>10.59</td><td>20.97</td><td>51.93</td><td>83.44</td></tr><tr><td>constant</td><td>1.01</td><td>2.00</td><td>4.97</td><td>9.94</td><td>20.05</td><td>50.83</td><td>80.36</td></tr></tbody></table>
**3.7. The mean of DCT high-frequency coefficients method.** It is well known that the energy in high-frequency coefficients of a DCT image block corresponds mainly to the noise and edges, while the medium and low frequency coefficients contribute to the geometry of the image. This method computes the variance of the noise from the high-frequency coefficients.

Since the orthonormal DCT-II is a isometry, the variance can be computed both at the spatial domain or the frequential domain. However, computing it at the frequential domain has the advantage that it decorrelates the signal and compacts its energy on a few coefficients, thus making it possible to estimate the noise using only high-frequency DCT-II coefficients.

A correction factor $k = 1.0061$ is needed to correct the obtained standard deviation.

The estimated $\tilde{\sigma}^2$ is obtained as the mean of the variance of coefficient at frequencies [6, 7], [7, 6], and [7, 7]. To compute the variance, the unbiased estimator for the population variance is used: $S_{N-1}^2(X) = \frac{1}{N-1} \sum_{i=1}^N (X - \bar{X})$, where $\bar{X} = \frac{1}{N} \sum_{i=1}^N x_i$ and $N$ the number of elements in $X$. The pseudo-code for the DCT method is given in Algorithm 8 and the results for white Gaussian noise in Table 11.

<caption>**Algorithm 8 DCT noise estimation algorithm.**</caption>

1: **DCT** - Return the standard deviation of the image noise. **Input** $U$: input image. **Output** $\tilde{\sigma}$: estimated standard deviation of the noise in the image.
2: $w = 8$.
3: $k = 1.0061426171829178$. ▷ Correction factor
4: $B \leftarrow$ all $w \times w$ block in $U$.
5: $N$ = number of blocks in $B$.
6: $C = 2$D-DCT($B$). ▷ Compute the 2D orthonormal DCT of each block in $B$
7: $C_1 \leftarrow C[6, 7]$. ▷ Obtain all coefficients at position [6, 7]
8: $C_2 \leftarrow C[7, 6]$. ▷ Obtain all coefficients at position [7, 6]
9: $C_3 \leftarrow C[7, 7]$. ▷ Obtain all coefficients at position [7, 7]
10: $\hat{\sigma} = \sqrt{\text{Mean}[\text{Var}(C_1), \text{Var}(C_2), \text{Var}(C_3)]}$. return $\tilde{\sigma} = k \times \hat{\sigma}$.
<table><caption>TABLE 11. The Mean of DCT high-frequency coefficients method results with uniform noise.</caption><thead><tr><th>Image / σ̃</th><th>σ = 1</th><th>σ = 2</th><th>σ = 5</th><th>σ = 10</th><th>σ = 20</th><th>σ = 50</th><th>σ = 80</th></tr></thead><tbody><tr><td>bag</td><td>2.59</td><td>3.10</td><td>5.58</td><td>10.31</td><td>20.30</td><td>50.24</td><td>80.66</td></tr><tr><td>building1</td><td>3.78</td><td>4.18</td><td>6.24</td><td>10.71</td><td>20.40</td><td>50.40</td><td>80.56</td></tr><tr><td>computer</td><td>3.51</td><td>3.94</td><td>6.06</td><td>10.57</td><td>20.47</td><td>50.62</td><td>80.52</td></tr><tr><td>dice</td><td>1.24</td><td>2.14</td><td>5.10</td><td>10.03</td><td>20.25</td><td>50.23</td><td>80.53</td></tr><tr><td>flowers2</td><td>1.47</td><td>2.28</td><td>5.16</td><td>10.10</td><td>20.05</td><td>50.27</td><td>79.97</td></tr><tr><td>hose</td><td>1.42</td><td>2.23</td><td>5.15</td><td>10.12</td><td>20.12</td><td>50.29</td><td>81.18</td></tr><tr><td>leaves</td><td>2.80</td><td>3.27</td><td>5.64</td><td>10.40</td><td>20.32</td><td>50.04</td><td>81.12</td></tr><tr><td>lawn</td><td>3.06</td><td>3.53</td><td>5.79</td><td>10.44</td><td>20.37</td><td>50.35</td><td>80.24</td></tr><tr><td>stairs</td><td>1.46</td><td>2.28</td><td>5.15</td><td>10.15</td><td>20.11</td><td>50.18</td><td>80.90</td></tr><tr><td>traffic</td><td>4.84</td><td>5.13</td><td>6.90</td><td>11.01</td><td>20.54</td><td>50.63</td><td>80.93</td></tr><tr><td>constant</td><td>1.00</td><td>2.02</td><td>5.04</td><td>10.04</td><td>20.15</td><td>49.84</td><td>80.36</td></tr></tbody></table>
<table><caption>TABLE 12. Mean of DCT high-frequency coefficients method with MAD estimator results with white Gaussian noise.</caption><thead><tr><th>Image / ˜σ</th><th>σ = 1</th><th>σ = 2</th><th>σ = 5</th><th>σ = 10</th><th>σ = 20</th><th>σ = 50</th><th>σ = 80</th></tr></thead><tbody><tr><td>bag</td><td>1.90</td><td>2.71</td><td>5.47</td><td>10.23</td><td>20.20</td><td>50.10</td><td>80.30</td></tr><tr><td>building1</td><td>1.71</td><td>2.77</td><td>5.75</td><td>10.59</td><td>20.23</td><td>50.13</td><td>80.29</td></tr><tr><td>computer</td><td>1.72</td><td>2.68</td><td>5.60</td><td>10.38</td><td>20.39</td><td>50.44</td><td>80.15</td></tr><tr><td>dice</td><td>1.16</td><td>2.10</td><td>5.08</td><td>9.98</td><td>20.17</td><td>49.92</td><td>80.28</td></tr><tr><td>flowers2</td><td>1.19</td><td>2.14</td><td>5.12</td><td>10.04</td><td>19.98</td><td>50.08</td><td>79.74</td></tr><tr><td>hose</td><td>1.32</td><td>2.19</td><td>5.11</td><td>10.06</td><td>20.03</td><td>50.07</td><td>81.21</td></tr><tr><td>leaves</td><td>2.35</td><td>3.00</td><td>5.57</td><td>10.36</td><td>20.26</td><td>49.98</td><td>81.14</td></tr><tr><td>lawn</td><td>2.60</td><td>3.23</td><td>5.70</td><td>10.41</td><td>20.30</td><td>50.22</td><td>79.87</td></tr><tr><td>stairs</td><td>1.39</td><td>2.25</td><td>5.13</td><td>10.11</td><td>19.94</td><td>50.07</td><td>80.46</td></tr><tr><td>traffic</td><td>2.25</td><td>3.23</td><td>6.06</td><td>10.74</td><td>20.43</td><td>50.63</td><td>80.48</td></tr><tr><td>constant</td><td>1.00</td><td>2.01</td><td>5.00</td><td>10.05</td><td>19.95</td><td>49.49</td><td>79.58</td></tr></tbody></table>

### 3.8. The Mean of DCT high-frequency coefficients method with MAD estimator.

We know of no particular reference for this method, except that in view of the other ones, it is simply necessary to test it. The Median of Absolute Deviations (MAD) estimator [61] is a robust estimator of the standard deviation that can be used when the data present outliers. The MAD estimator of X is defined as MAD(X) = median(|X - median(X)|). If X follows a normal distribution, MAD is a biased estimator of the standard deviation and must be multiplied by k = 1.4865 to get back ˜σ. This method is essentially the same as in Section 3.7, but using the MAD estimator instead of directly computing the sampled variance. The pseudo-code for the DCT-MAD method is given in Algorithm 9 and the results for white Gaussian noise in Table 12.

<caption>**Algorithm 9** DCT-MAD noise estimation algorithm.</caption>

1: **DCT** - Return the standard deviation of the image noise. **Input** **U**: input image. **Output** ˜σ: estimated standard deviation of the noise in the image.
2: w = 8.
3: k = 1.486495980422939.
4: B ← all w × w block in **U**.
5: C = 2D-DCT(B).
6: C₁ ← C[6, 7].
7: C₂ ← C[7, 6].
8: C₃ ← C[6, 7].
9: ̂σ = Mean[MAD(C₁), MAD(C₂), MAD(C₃)].
return ˜σ = k × ̂σ.
**3.9. Deconvolution method (E.I.N.V.).** The Estimation of Image Noise Variance method [5] (E.I.N.V) first computes the finite differences derivative of the image, then makes a histogram of the local standard deviations of these high pass filtered values, and finally it evaluates the histogram iteratively in order to converge to the noise standard deviation estimate.

**Step 1: suppression of the original image x**

Given an image $y(m, n)$ corrupted with Gaussian noise $w(m, n)$, $y(m, n) = x(m, n)+w(m, n) \forall (m, n) \in$ the domain of the ideal image $x$:

First, the horizontal derivative is computed:

$$y_1(m, n) = \frac{1}{\sqrt{2}} [y(m + 1, n) - y(m, n)]$$

Then, the vertical derivative of the horizontal derivative $y_1$:

$$y_2(m, n) = \frac{1}{\sqrt{2}} [y_1(m, n + 1) - y_1(m, n)].$$

**Step 2: computing the histogram of local standard deviations** Using a window of size $L \times L$ pixels, with $L = 2K + 1$, $K \in \mathbb{N}$, the variance is estimated by the classic unbiased estimator

$$\hat{\sigma}^2(m, n) = \frac{1}{N_L - 1} \left( \left[ \sum_{i=-K}^{K} \sum_{j=-K}^{K} y_2^2(m + i, n + j) \right] - N_L \hat{\mu}^2(m, n) \right).$$

$N_L$ is the number of pixels in the window ($L^2$).

The local mean is computed as:

$$\hat{\mu}(m, n) = \frac{1}{N_L} \sum_{i=-K}^{K} \sum_{j=-K}^{K} y_2(m + i, n + j).$$

The authors explain that they tested several different window sizes and found that the best results are achieved with the minimal possible window size, $L = 3$.

The histogram $h(k)$ with entries $\hat{\sigma}(m, n)$ is defined as:

$$h(k) = \begin{cases} \#\{(m, n) : k - \frac{1}{2} \le \alpha \hat{\sigma}(m, n) < k + \frac{1}{2}\}, & \text{if } k = 1, \dots, k_{\max}; \\ 2\#\{(m, n) : 0 \le \alpha \hat{\sigma}(m, n) < \frac{1}{2}\}, & \text{if } k = 0, \end{cases}$$

where “#” means the cardinal of the set.

The value $\alpha = 1$ was considered adequate by the authors for most of the images.

**Step 3: Evaluation of the histogram**

An averaging of the values $\hat{\sigma}^2(m, n)$ can be performed by computing the mean square value of the histogram,

$$s_l^2 = \frac{\sum_{k=0}^{k_{\max}} k^2 h(k)}{\sum_{k=0}^{k_{\max}} h(k)}.$$

The value $s_l^2/\alpha^2$ is an initial global estimate for the noise variance $\sigma_w^2$.

The estimation can be improved iteratively thanks to the prior knowledge of the shape which the histogram should have. The histogram should exhibit a rapid descent for large values of $k$, but for a natural image, the descent is significantly slower. The reason they give for this fact is that the histogram is the convolution of a noisy uniform image with the histogram of the ideal image.
<table><caption>TABLE 13. Method for the Estimation of Image Noise Variance (deconvolution method) results with white Gaussian noise.</caption><thead><tr><th>Image / ˜σ</th><th>σ = 1</th><th>σ = 2</th><th>σ = 5</th><th>σ = 10</th><th>σ = 20</th><th>σ = 50</th><th>σ = 80</th></tr></thead><tbody><tr><td>bag</td><td>3.20</td><td>3.82</td><td>6.45</td><td>11.12</td><td>20.65</td><td>50.35</td><td>80.59</td></tr><tr><td>building1</td><td>1.58</td><td>2.53</td><td>5.79</td><td>11.14</td><td>20.92</td><td>50.57</td><td>80.46</td></tr><tr><td>computer</td><td>2.02</td><td>2.87</td><td>5.99</td><td>11.09</td><td>21.05</td><td>50.76</td><td>80.68</td></tr><tr><td>dice</td><td>1.07</td><td>2.12</td><td>5.11</td><td>10.08</td><td>20.12</td><td>50.25</td><td>80.63</td></tr><tr><td>flowers2</td><td>1.12</td><td>2.15</td><td>5.16</td><td>10.15</td><td>20.17</td><td>50.36</td><td>80.33</td></tr><tr><td>hose</td><td>1.68</td><td>2.56</td><td>5.38</td><td>10.24</td><td>20.25</td><td>50.08</td><td>80.63</td></tr><tr><td>leaves</td><td>3.60</td><td>4.12</td><td>6.45</td><td>11.01</td><td>20.64</td><td>50.26</td><td>80.65</td></tr><tr><td>lawn</td><td>3.51</td><td>4.03</td><td>6.39</td><td>10.95</td><td>20.60</td><td>50.33</td><td>80.51</td></tr><tr><td>stairs</td><td>1.84</td><td>2.67</td><td>5.42</td><td>10.29</td><td>20.14</td><td>50.36</td><td>80.48</td></tr><tr><td>traffic</td><td>2.14</td><td>2.98</td><td>6.14</td><td>11.41</td><td>21.21</td><td>50.79</td><td>80.64</td></tr><tr><td>constant</td><td>1.00</td><td>2.00</td><td>5.04</td><td>10.00</td><td>20.03</td><td>50.28</td><td>80.40</td></tr></tbody></table>

Therefore, it is a deconvolution problem. The method fades-out the histogram by a descending weighting function $g(k)$. The following function is used:

$$g_l(k) = \begin{cases} 1 & \text{if } k \le s_l \\ \frac{1}{2} \left[ 1 - \cos \left( \frac{\beta - \frac{k}{s_l}}{1-\beta} \right) \right] & \text{if } s_l < k < \beta s_l \\ 0 & \text{if } k \ge \beta s_l. \end{cases}$$

The value of $\beta$ must be obtained experimentally. The authors of the method found that $\beta = 2.12$ is adequate. This value was used in all tests for this method. An improved value of the square mean $s^2$ can be obtained iteratively by

$$s_{l+1}^2 = \frac{\sum_{k=0}^{k_{max}} k^2 g_l(k) h(k)}{\sum_{k=0}^{k_{max}} g_l(k) h(k)}.$$

A fixed value $l_{max} = 4$ is enough. The initial value $s_l$ if taken from $s_l^2 = \frac{\sum_{k=0}^{k_{max}} k^2 h(k)}{\sum_{k=0}^{k_{max}} h(k)}$. Finally, the estimate $\tilde{\sigma}^2 = \frac{s_{l_{max}}}{\alpha^2}$. It was found that a correction factor $F = 1.0141$ is needed to get an unbiased estimation. The pseudo-code for this method is given in Algorithm 10 and the results for white Gaussian noise in Table 13.

This is comparable to the other ones, in that it computes the average of a weighted quantile of the variance histogram. The chosen differential operator, the second cross derivative, is judicious. The results are bad for low noise. They show that either the high pass filter has not a large enough order, or the quantile is not low enough.
**Algorithm 10** Estimation of Image Noise Variance (deconvolution method) algorithm.

1: **EINV** - Return the standard deviation of the image noise. **Input U**: input image. **Output $\tilde{\sigma}$**: estimated standard deviation of the noise in the image.

2: $F = 1.0140974819332935$.

▷ Correction factor

3: $y_1(m, n) = \frac{1}{\sqrt{2}}[y(m + 1, n) - y(m, n)]$.

4: $y_2(m, n) = \frac{1}{\sqrt{2}}[y_1(m, n + 1) - y_1(m, n)]$.

5: $\hat{\sigma}^2(m, n) = \frac{1}{N_L - 1} \left( \left[ \sum_{i=-K}^{K} \sum_{j=-K}^{K} y_2^2(m + i, n + j) \right] - N_L \hat{\mu}^2(m, n) \right)$.

6: Create histogram $h(k)$.

▷ View step 2

7: $Q = \text{ones}(k_{\text{max}} + 1)$.

▷ Array of $k_{\text{max}}$ ones

8: **for** $i = 1...4$ **do**

9: $N = 0$.

▷ Numerator

10: $D = 0$.

▷ Denominator

11: **for** $k = 0, ... k_{\text{max}}$ **do**

12: $N = N + k^2 Q(k) h(k)$

13: $D = D + Q(k)h(k)$

14: **end for**

15: $s_l = \sqrt{N/D}$

16: **for** $k = 0, ..., k_{\text{max}}$ **do**

17: **if** $k < s_l$ **then**

18: $g_l(k) = 1$

19: **end if**

20: **if** $s_l < k$ and $k < \beta s_l$ **then**

21: $g_l(k) = 0.5 \left[ 1 - \cos \left( \frac{\pi(\beta - k/s_l)}{1 - \beta} \right) \right]$

22: **end if**

23: **if** $k \ge \beta s_l$ **then**

24: $g_l(k) = 0$

25: **end if**

26: **end for**

27: $Q = g_l$.

▷ Copy array

28: $\tilde{\sigma} = (s_l/\alpha)F$

29: **end for**

return $\tilde{\sigma}$
$$L_1 = \begin{array}{|c|c|c|} \hline 0 & 1 & 0 \\ \hline 1 & -4 & 1 \\ \hline 0 & 1 & 0 \\ \hline \end{array} \qquad L_2 = \frac{1}{2} \begin{array}{|c|c|c|} \hline 1 & 0 & 1 \\ \hline 0 & -4 & 0 \\ \hline 1 & 0 & 1 \\ \hline \end{array}$$

<caption>FIGURE 7. Laplace-like operators.</caption>

$$N = 2(L_2 - L_1) = \begin{array}{|c|c|c|} \hline 1 & -2 & 1 \\ \hline -2 & 4 & -2 \\ \hline 1 & -2 & 1 \\ \hline \end{array}$$

<caption>FIGURE 8. Difference of Laplacian operators to get the noise estimation operator in the article *Fast Noise Variance Estimation* [4].</caption>

**3.10. Method for Fast Noise Variance Estimation.** The Fast Noise Variance Estimation method [4] tries to avoid the interference of image structures (edges and textures) present on the image when estimating the noise. To do it, it detects these structures using an operator based on the Laplacian, and cancels them. It considers two 3 × 3 Laplacian stencils $L_1$ and $L_2$ (Figure 7) and makes their difference to obtain the noise estimation operator $L = 2(L_2 - L_1)$ (Figure 8).

Operator $L$ cancels edges and structures in the image, but it is a biased estimator of the variance. The correction factor to get an unbiased estimation of the variance is $\frac{1}{36}$. The pseudocode for this method is given in Algorithm 11 and the results for white Gaussian noise in Table 14.

**Algorithm 11 Fast Noise Variance Estimation algorithm.**

1: **FNVE** - Return the standard deviation of the image noise. **Input U**: input image. **Output $\tilde{\sigma}$**: estimated standard deviation of the noise in the image.

2: $L = [[1, -2, 1], [-2, 4, -2], [1, -2, 1]]$.

3: $V = 0$.

4: **for** $[x, y] \in [1, N_x - 2] \times [1, N_y - 2]$ **do**

5: &nbsp;&nbsp;&nbsp;&nbsp;A = 0.

6: &nbsp;&nbsp;&nbsp;&nbsp;**for** $(i, j) \in (-1, \dots, 1, -1, \dots, 1)$ **do**

7: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$A = A + U(x + i, y + j)L(1 + i, 1 + j)$.

8: &nbsp;&nbsp;&nbsp;&nbsp;**end for**

9: $V = V + A^2$.

10: **end for**

$$\text{return } \tilde{\sigma} = \sqrt{\frac{V}{36(N_x - 1)(N_y - 1)}}.$$
<table><caption>TABLE 14. Fast Noise Variance Estimation method results with simulated homoscedastic white Gaussian noise.</caption><thead><tr><th>Image / σ̃</th><th>σ = 1</th><th>σ = 2</th><th>σ = 5</th><th>σ = 10</th><th>σ = 20</th><th>σ = 50</th><th>σ = 80</th></tr></thead><tbody><tr><td>bag</td><td>3.61</td><td>4.00</td><td>6.09</td><td>10.60</td><td>20.29</td><td>49.94</td><td>80.14</td></tr><tr><td>building1</td><td>5.00</td><td>5.30</td><td>7.01</td><td>11.14</td><td>20.57</td><td>50.15</td><td>80.11</td></tr><tr><td>computer</td><td>5.21</td><td>5.50</td><td>7.17</td><td>11.21</td><td>20.66</td><td>50.42</td><td>80.10</td></tr><tr><td>dice</td><td>1.34</td><td>2.19</td><td>5.09</td><td>10.01</td><td>20.05</td><td>49.94</td><td>80.25</td></tr><tr><td>flowers2</td><td>1.79</td><td>2.49</td><td>5.21</td><td>10.09</td><td>20.00</td><td>50.11</td><td>79.81</td></tr><tr><td>hose</td><td>1.76</td><td>2.47</td><td>5.21</td><td>10.11</td><td>20.09</td><td>49.92</td><td>80.33</td></tr><tr><td>leaves</td><td>3.65</td><td>4.03</td><td>6.10</td><td>10.61</td><td>20.36</td><td>49.84</td><td>80.36</td></tr><tr><td>lawn</td><td>3.78</td><td>4.16</td><td>6.18</td><td>10.63</td><td>20.34</td><td>50.03</td><td>80.06</td></tr><tr><td>stairs</td><td>1.81</td><td>2.51</td><td>5.22</td><td>10.13</td><td>20.03</td><td>50.05</td><td>80.25</td></tr><tr><td>traffic</td><td>6.02</td><td>6.25</td><td>7.76</td><td>11.59</td><td>20.80</td><td>50.41</td><td>80.24</td></tr><tr><td>constant</td><td>1.00</td><td>2.00</td><td>5.02</td><td>9.96</td><td>19.92</td><td>50.04</td><td>80.22</td></tr></tbody></table>
<table><caption>TABLE 15. Wavelet DB3 MAD results with white Gaussian noise.</caption><thead><tr><th>Image / ˜σ</th><th>σ = 1</th><th>σ = 2</th><th>σ = 5</th><th>σ = 10</th><th>σ = 20</th><th>σ = 50</th><th>σ = 80</th></tr></thead><tbody><tr><td>bag</td><td>2.37</td><td>3.12</td><td>5.78</td><td>10.48</td><td>20.03</td><td>49.81</td><td>79.58</td></tr><tr><td>building1</td><td>1.77</td><td>2.85</td><td>5.91</td><td>10.71</td><td>20.39</td><td>49.85</td><td>79.52</td></tr><tr><td>computer</td><td>1.99</td><td>2.97</td><td>5.92</td><td>10.79</td><td>20.28</td><td>49.74</td><td>79.90</td></tr><tr><td>dice</td><td>1.11</td><td>2.08</td><td>5.03</td><td>9.96</td><td>19.80</td><td>49.81</td><td>78.97</td></tr><tr><td>flowers2</td><td>1.16</td><td>2.11</td><td>5.05</td><td>9.98</td><td>20.07</td><td>49.33</td><td>79.08</td></tr><tr><td>hose</td><td>1.47</td><td>2.33</td><td>5.17</td><td>10.06</td><td>19.85</td><td>49.67</td><td>79.27</td></tr><tr><td>leaves</td><td>2.89</td><td>3.50</td><td>5.86</td><td>10.44</td><td>20.08</td><td>49.64</td><td>79.76</td></tr><tr><td>lawn</td><td>3.06</td><td>3.68</td><td>6.07</td><td>10.61</td><td>20.19</td><td>49.64</td><td>79.43</td></tr><tr><td>stairs</td><td>1.77</td><td>2.57</td><td>5.31</td><td>10.16</td><td>20.03</td><td>49.13</td><td>79.68</td></tr><tr><td>traffic</td><td>2.25</td><td>3.29</td><td>6.24</td><td>11.01</td><td>20.67</td><td>50.22</td><td>79.66</td></tr><tr><td>Flat image</td><td>0.99</td><td>1.96</td><td>4.92</td><td>9.92</td><td>19.63</td><td>48.85</td><td>78.59</td></tr></tbody></table>

**3.11. Wavelet MAD.** This method consists on applying the MAD estimator to the detail (HH) coefficients of the Discrete Wavelet Transform (DWT) of the image.

The algorithmic description for this method can be found in Algorithm 12 and its results for white Gaussian noise in Table 15. The DB3 Daubechies wavelet were used.

**Algorithm 12 Wavelet-MAD noise estimation algorithm.**

1: **Wavelet-MAD** - Return the standard deviation of the image noise. **Input U**: input image. **Output ˜σ**: estimated standard deviation of the noise in the image.

2: $w = 8$.

3: Obtain multiresolution wavelet coefficients LL, HL, LH, HH of U.

4: $M = \text{median}(\text{HH})$.

$$ \text{return } \tilde{\sigma} = \frac{\text{median}(|\text{HH} - M|)}{0.6745} $$

**3.12. Conclusion.** In this chapter we presented a review of classic homoscedastic noise estimation methods.

Table 16 gives the RMSE values obtained by evaluating these methods with white Gaussian noise, with $\sigma \in \{1, 2, 5, 10, 20, 50, 80\}$ along all the images. This table is the final result of this study of homoscedastic noise estimation methods, and is particularly decisive for the low noise values. Giving the RMSE depending on $\sigma$ for white Gaussian noise is perhaps not quite intuitive, and requires some attention. For example, a RMSE of 1.0 is an excellent estimate for $\sigma = 80$, but is quite bad if $\sigma = 2$. Thus this last table shows that the results are very bad for all methods but the Ponomarenko et al. and the Percentile methods. Still the error is about 8% for low level noise $\sigma = 2, 5$, but it only doubles for much larger noises, which means that it is ridiculously low for high noise. These precisions on the noise estimation are more than enough for applying denoising algorithms.
<table><caption>TABLE 16. Comparison of the averaged RMSE along all test images in Figure 3 using simulated homoscedastic white Gaussian noise.</caption><thead><tr><th>Method</th><th>σ = 1</th><th>σ = 2</th><th>σ = 5</th><th>σ = 10</th><th>σ = 20</th><th>σ = 50</th><th>σ = 80</th></tr></thead><tbody><tr><td>PCA [3]</td><td><strong>0.31</strong></td><td><strong>0.21</strong></td><td><strong>0.13</strong></td><td><strong>0.25</strong></td><td>0.37</td><td>1.17</td><td>2.57</td></tr><tr><td>DCT-MAD</td><td>0.59</td><td>0.85</td><td>0.48</td><td>0.29</td><td>0.12</td><td>0.38</td><td>0.71</td></tr><tr><td>DCT</td><td>0.99</td><td>2.21</td><td>0.71</td><td>0.31</td><td>0.19</td><td>0.46</td><td>0.74</td></tr><tr><td>Ponomarenko et al. [15, 16]</td><td>0.56</td><td>0.39</td><td>0.28</td><td>0.34</td><td>0.49</td><td>1.38</td><td>1.60</td></tr><tr><td>Wavelet MAD</td><td>0.91</td><td>0.91</td><td>0.73</td><td>0.43</td><td><strong>0.05</strong></td><td><strong>0.25</strong></td><td>0.83</td></tr><tr><td>F.N.V.E. [4]</td><td>1.64</td><td>3.12</td><td>1.40</td><td>0.45</td><td>0.17</td><td>0.26</td><td><strong>0.25</strong></td></tr><tr><td>Percentile [17]</td><td>0.58</td><td>0.52</td><td>0.51</td><td>0.47</td><td>0.67</td><td>1.22</td><td>2.07</td></tr><tr><td>E.I.N.V. [5]</td><td>1.36</td><td>0.72</td><td>1.01</td><td>0.67</td><td>0.39</td><td>0.56</td><td>0.59</td></tr><tr><td>Block [56, 57]</td><td>1.47</td><td>0.41</td><td>0.94</td><td>1.05</td><td>1.27</td><td>1.00</td><td>1.05</td></tr><tr><td>Gradient [58, 59]</td><td>2.44</td><td>0.96</td><td>1.45</td><td>1.30</td><td>0.54</td><td>1.17</td><td>1.76</td></tr><tr><td>Scatter [60]</td><td>7.24</td><td>0.58</td><td>1.75</td><td>1.39</td><td>2.23</td><td>1.56</td><td>2.26</td></tr><tr><td>Average [55]</td><td>3.02</td><td>2.15</td><td>2.31</td><td>1.64</td><td>1.43</td><td>1.63</td><td>0.99</td></tr><tr><td>Median [55]</td><td>3.26</td><td>2.35</td><td>2.54</td><td>1.83</td><td>1.81</td><td>2.17</td><td>1.61</td></tr></tbody></table>

In table 16 they appear three additional methods: the Ponomarenko et al., the Percentile, and the PCA methods. These methods have a special interest because their performance is clearly superior to that of the classic method (the case of the Ponomarenko et al. and the Percentile methods), or because they are relatively recent (the PCA method). These methods have been studied with great detail in the third part of the thesis. The Ponomarenko et al. method is detailed in Chapter 10, the Percentile method in Chapter 11, and the PCA method in Chapter 12, and evaluated in Chapter 13.

The method Estimation of Image Noise Variance is referred to as E.I.N.V. [5] and the method Fast Noise Variance Estimation as F.N.V.E [4]. The Pyramid method [6] is not evaluated here, because it was not able to give an estimate for all σ levels. The Ponomarenko et al. method was configured with $p = 0.005$ and $w = 8$.

It can be concluded that for simulated homoscedastic white Gaussian noise, the PCA methods gives better results than its competitors, for low and moderated noises ($\sigma < 10$). For $\sigma > 10$, simpler methods, as F.N.V.E and the Wavelet MAD methods give the best results.

It is important to note that this table only shows the preliminary results of the study presented in this thesis, where the noise is homoscedastic, white, and simulated. As it will be explained in Chapter 2, simulating this kind is simple noise is not realistic at all, since real images are (at least) signal-dependent, and thus the variance of the noise increases with the intensity. Adapting the method to signal-dependent noise implies changes that make this classification change significantly. In fact, the classic method gives inaccurate results when adapted, whereas the Ponomarenko et al. method gives the best results, followed by Percentile and PCA.
In Chapter 3 it will be shown than the signal-dependent noise is not enough (for example, it
does not apply to digital images after demosaicing the raw image at the focal plane). In Chapter
13 it will be shown that after adapting the methods to signal-dependent noise, the PCA method
is not reliable when the image contains high noise, and needs much more samples/bin than the
Ponomarenko et. al and Percentile methods, which perform much accurately.
<figure><img src="image_19.png" /></figure>
# Signal-dependent noise estimation

Optimal denoising works at best on raw images (the image formed at the output of the focal plane, at the CCD or CMOS detector), which display a white signal-dependent noise. The noise model of the raw image is characterized by a function that given the intensity of a pixel in the noisy image returns the corresponding STD; the plot of this function is the noise curve. This chapter develops a non-parametric approach estimating the noise curve directly from a single raw image. An extensive cross-validation procedure is described to compare this new method with state-of-the-art parametric methods and with laboratory calibration methods giving a reliable ground-truth, even for nonlinear detectors.

The signal-dependent noise model is valid for raw images, but when the noise is correlated and thus frequency-dependent (for example, after demosaicing the raw image), the noise model presented here is not enough, as will be shown in Chapter 3. Chapter 4 will discuss a new model able to measure the noise even in JPEG-encoded images.

## 1. Introduction

Most denoising methods assume that the noise in the image is additive, homoscedastic, white, and Gaussian. Homoscedastic means that the variance of the Gaussian noise is fixed and does not depend on the pixel position or value. By “white” noise, we mean that the noise pixel values are independent (look at Section 3 of Chapter 1 for a review of classic homoscedastic white noise estimators). We shall retain this terminology throughout.

The homoscedasticity assumption is not realistic. The photon emission by a body follows a Poisson distribution which can be approximated by a Gaussian distribution when the number of photons is large enough. But the variance of this Gaussian is signal dependent. In the Poisson model [62, 63, 64, 65, 66, 67, 68], an image value $\tilde{U}(x, y)$ at pixel $(x, y)$ is a Poisson variable with variance and mean equal to $U(x, y)$, where $U$ is the ideal noise-free image. The Poisson noise has therefore a standard deviation (STD) equal to $[U(x, y)]^{1/2}$. Thus, an ideal raw image is a white Poisson noise whose mean at each pixel is the noiseless value. Note that this is related to the quantum nature of light and the probability of emitting a photon, independently of the technology used at the CFA (CCD, CMOS). This Poisson noise adds up to a thermal noise and to an electronic noise which are approximately additive and white, making the final noise model not necessarily Poisson distributed, but still white and signal dependent.

Noise estimation is a necessary preliminary step for most image processing and computer vision algorithms [1]. Nevertheless, several other denoising methods propose to deal directly with
Poisson noise. Wavelet-based denoising methods [45, 44] propose to adapt the transform threshold to the local noise level of the Poisson process. Lefkimmiatis et al. [46] have explored a Bayesian approach and Deledalle et al. [47] have adapted the Non-local means algorithm [32] to Poisson noise. These papers assume that no variance stabilizing transform (VST) transforming the signal dependent noise into a nearly homoscedastic noise is accurate enough to transform the Poisson noise into homoscedastic noise. The advantage of VSTs is that they permit the application of a classic denoising algorithm. The VST associated with Poisson noise is often called Anscombe transform [42], but one can attach a VST to any signal dependent noise model [1]. As a matter of fact, papers on the Anscombe transform [43] (for low count Poisson noise) and [50] (for Rician noise) argue that, when combined with suitable forward and inverse variance stabilizing transformations, algorithms designed for signal independent Gaussian noise work just as well as ad-hoc algorithms for Poisson noise models. These considerations confirm the importance of estimating as accurately as possible the noise curves of raw images, since their accurate knowledge is required to compute the VST. In most CCDs and CMOS detectors, the variance of the noise at a pixel is approximated (assuming that all detectors at the CFA are equivalent and thus neglecting the fixed pattern noise) by a simple linear model $\sigma^2 = A + BU$, where $U$ is the expectation of the intensity of this pixel in the noisy image. This model is valid under the assumption mentioned above of a combination of a Poisson with a thermal noise. Yet, this assumption holds only if the signal is not saturated and the photon count large enough. At the darkest pixels, the Poisson distribution of the noise cannot be approximated by a Gaussian and it becomes a *shot* noise. In short, the noise variance does not necessarily follow the linear model in the darkest and brightest image regions. An accurate estimation of the noise at the darkest zones is crucial since subsequent processes in the camera chain (specially, the gamma correction step) are designed to increase the dynamics in the dark zones. If the noise is not removed at the raw image stage, it might end up really augmented at the final stage.

Parametric noise estimation methods try to obtain the parameters that control a noise model (for example, the *A* and *B* parameters of the linear model). Yet, to get a realistic estimation, they have to take into account the effect of the saturation in the darkest and brightest pixels in the final noise curve. To validate the estimation of a noise estimation method, its noise curve must be compared to a ground-truth curve. Such a ground-truth for a particular camera and settings can be obtained by taking a series of photographs of a pattern, that is mostly flat and contains a wide range of gray levels. The temporal variation of the gray level at a given pixel gives an estimate of the noise STD associated with this gray level. However, the series of photographs must be taken under controlled conditions, to ensure that any variation of the intensity of a pixel can be only explained by the noise. In short, it is a heavy procedure (that is, it requires constant lighting, a camera stabilizer to fix its position, and isolation from any kind of electromagnetic source that may introduce electronic noise into the camera) which also needs access to the camera that took the photographs. It also requires the *a priori* knowledge of the form of the camera noise model, which is not granted. This explains why the establishment of a method able to estimate automatically the noise model from a single snapshot is a valid question. Furthermore, if the method can be shown
to be reliable even without any *a priori* model guess, its credibility will be somewhat augmented. In this chapter, we show that it is indeed possible to use a non-parametric estimator to get an accurate noise curve from the noisy image itself, by measuring the variance locally with patch-based methods [69, 70, 71, 3, 56, 72]. This eliminates the need for lab calibration procedure. Indeed, the procedure described uses one or several photographs taken in arbitrary environment and yields a non-parametric noise model as good (for those images) as the one obtained by the heavier ground truth procedure (laboratory calibration). We also examine the question of whether it is better to use a parametric or a non parametric model when dealing with a single or a few photographs. Our conclusion is that the non parametric method gives results comparable to the parametric method, but is somewhat less risky as it does not propagate local estimation errors caused by the presence of texture in the image.

Our plan follows from the above discussion. Since noise estimation is a well-known procedure for white homoscedastic noise, Section 2 will review the literature on white homoscedastic noise estimation and will point out competitive algorithms. Section 3 explains the procedure that should be followed to get a reliable non-parametric noise curve from a series of images, under controlled conditions. Section 4 discusses how homoscedastic white noise estimation algorithms can be adapted to estimate an arbitrary signal-dependent noise curve. Section 5 compares the Root Mean Squared Errors (RMSE) between the non-parametric ground-truth, the STDs from the series of images and two state of the art parametric methods. Finally, Section 6 presents the conclusions, that validate our proposed nonparametric method, but also the use of two state-of-the-art parametric methods.

## 2. State-of-the-art in white noise estimation

Many noise estimation methods share the following features, which can be summarized in two sentences:

* estimate noise in high frequencies, where noise dominates over signal;

* estimate noise in image regions with the least variation, typically the blocks with the smallest STDs.

Thus, these numerous methods [4, 5, 6, 16, 56, 72, 55, 57, 60, 58, 59, 73, 74] proceed roughly as follows:

* they start by applying some high-pass filter, which concentrates the image energy on its edges, while the noise remains spatially homogeneous;

* they compute the energy of many blocks extracted from this high-passed image;

* they estimate the STDs of the blocks;

* to avoid blocks whose STD is mostly explained by the underlying ideal image, a statistic robust to (many) outliers must be applied. The methods therefore prefer the flattest blocks, which belong to a (low) percentile of the STDs of all the blocks.

Note that the power spectral density of a natural image is not homogeneous. Most of the energy corresponding to its geometry is located at the low and medium frequencies, whereas high
frequency coefficients bring little visual information (with the exception of the edges). Conversely, an image can be considered "highly textured" if the energy at the high-frequency coefficients is as high as the energy observed at edges. Thus, high-passing the image before estimating the noise spatially (or equivalently, estimating the noise only at the high-frequency DCT coefficients) is an initial step for many noise estimation algorithms. This enhances the contribution of the noise. Yet, avoiding edges and textures in the estimation remains necessary.

We shall limit ourselves to discuss the method acknowledged as the best estimator for ho-
moscedastic noise in the review [1], the Ponomarenko et al. method [16], along with the two of
the most competitive parametric methods for noise estimation in raw images [7, 75]. We briefly
describe these competitors in the next paragraphs. For a complete review on noise estimation
methods, we refer the reader to [1] and [55].

2.0.0.1. The Ponomarenko et al. approach. The Ponomarenko et al. [16] method is an extension of the previous method [72], based on the analysis of the DCT coefficients. In short, the Ponomarenko et al. method computes the variance of the high-frequency coefficients of a set of blocks whose variance measured at the low frequencies is minimal. We refer the reader to Chapter 10, where this method is analyzed in deep detail.

We now discuss two parametric methods that will be compared here.

2.0.0.2. Practical Poissonian-Gaussian Noise Modeling and Fitting for Single-image Raw-data.
Foi et al. proposed a simple parametric noise model [7] that takes into account the non-linear
response of the CCD due to the saturation of the signal and noise at the darkest and brightest
pixels of the image. The model assumes the well-known normal approximation, for which the
Poisson distribution $\mathcal{P}(\lambda)$ of the noise can in practice be approximated by the normal distribution:
$\mathcal{P}(\lambda) \rightarrow \mathcal{N}(\mu = \lambda, \sigma = \lambda)$. The method has two stages. In the first step it estimates several pairs
intensity/STD that form a scatter plot. In the second step, the observed pairs are used to fit a
global parametric model. Before applying these two steps, the image is preprocessed. First, the
2D-wavelet transform of the image is computed and the wavelet detail coefficients stored. The 1D
Daubechies wavelet and scaling functions are used to create the 2D kernels of the transform. The
STD of the noise is obtained from the detail coefficients of the transformed signal. In order to
be robust against edges, the image is segmented into level sets according to the intensity. Since
the image to be segmented is noisy, the segmentation is done in a low-pass filtered version. With
the selected regions of the image, the intensity of each pair is obtained as the sample mean of the
approximation wavelet coefficients and the estimated variance with the unbiased sample variance
estimator. The last step of the method is to fit the A and B parameters of the linear model of the
variance, for which a maximum-likelihood (ML) fitting is performed. However, since saturation
makes the response of the CCD or CMOS detector non-linear, the method needs to modify the
expectation and variance estimators to take saturation into account. The authors calculated the
new estimators from the distribution of the non-saturated signal and gave the explicit expression for
the expectation and variance estimators under saturation. Finally, these new pairs are incorporated
into the ML fitting in order to get the A and B parameters of the linear model despite the presence
of saturation. The model is able to predict the non-linear response of the CCD or CMOS detector
under saturation, giving explicitly the variance of the clipped noise for any intensity. Therefore, this method will be used as an example where parametric and nonparametric methods are cross-validated (Section 5).

2.0.0.3. Image Informative Maps for Component-wise Estimating Parameters of Signal-dependent Noise. In the paper [75] Uss et al. propose to adapt the use of disjoint informative maps [76] to estimate a parametric signal-dependent Poisson-like noise model. It discriminates between two kinds of non-overlapping blocks (SW – scanning window): those which belong to textures (TI – texture informative) and those that are suitable for noise estimation (NI – noise informative). To describe the textures of a given SW in the image, the 2D fractal Brownian motion (fBm) model is used, since the model is able to characterize a texture with few parameters. The roughness of the texture is obtained from the Hurst exponent in the fBm model. The estimation of the noise is obtained from a limited set of high-frequency coefficients of the DCT transform of the SWs that belong to the NI map. This idea was introduced in the Ponomarenko et al. method [16] and stated as the state of the art technique for noise estimation [1]. The Cramér-Rao Lower Bound (CRLB) is used to decide if a SW belongs to the TI or NI maps, on the texture parameters and the noise STD of the SW. All the SWs in the image are sorted according to increasing CRLB and then compared to a threshold. The SW below the threshold have the lowest CRLB and therefore belong to the NI map. The rest are assumed to be textures and assigned to the TI map. Since the criterion based on the CRLB relies on the (unknown) texture and noise parameters, the method begins with an initial guess for the NI and TI maps by fixing a noise STD and texture level to have an initial and rough CRLB criterion. Then, with the available CRLB criterion better STD and texture levels are computed, allowing for an even better CRLB criterion. The refining loop is iterated until convergence is reached. To estimate signal-dependent noise, the set of SW is partitioned into disjoint intensity sets according to their mean intensity and the method is applied separately to each set in order to get an (intensity-STD) pair. Therefore, this method is coherent with the claim we make in Section 4, which states that any block-based homoscedastic noise estimation method can be easily adapted to deal with signal-dependent noise, just by splitting the whole set of image blocks of the image into disjoint in intensity sets to apply then the homoscedastic version of the method to each of these sets. For example, if the input image has size $N_x \times N_y$, there are $M = (N_x - w + 1)(N_y - w - 1)$ overlapping blocks, that may be distributed into a set of $M/k$ bins, where each bin contains $k$ image blocks/bin whose mean intensity is a part of complete intensity range of the image.

### 3. Non-parametric noise ground-truth curve

Parametric methods fix an *a priori* model for the noise. For example, at the output of the CCD or CMOS detector a good approximation of the Poisson noise is to use the normal distribution approximation, at least when the number of incoming photons is large enough. Therefore, the variance of the noise is equal to the expectation. The noise at the output of the CCD or CMOS detector is Poissonian and therefore its variance is linear with the intensity. Also, thermal and electronic noise are added, and the noisy signal is amplified afterwards. Thus, the variance of the
noise can be modeled as a function of the intensity of the ideal (noise-free) image: $\sigma^2(\mathbf{U}) = A+B\mathbf{U}$. However, since the dynamics of the digital output from the CCD or CMOS detector is limited, the darkest and brightest pixels of the image can get saturated because of the noise, that becomes clipped noise. Because of the saturation, the probability distribution of the noise is no longer a symmetric normal distribution, but a truncated version with different statistics. The variance of the truncated distribution does not coincide with that of the normal distribution. Therefore, any realistic parametric estimation method must take into account that under saturation the simpler linear model is no longer valid. Some methods [7] adapt their expectation and variance estimators in order to take into account the effect of the saturation before fitting the linear function, while others [75] try to fit with polynomials of higher order or transform the image in such a way that the linear model holds.

In any case, the parametric model has to be validated in order to ensure that the curve they provide is indeed a function that accurately relates the intensity of the ideal image with the STD of the added noise. To do it, the estimations of the parametric method must be compared with a ground-truth noise curve. For the construction of the ground-truth the constraint of using just a single image is not needed. Indeed, it can be built from a series of snapshots of a calibration pattern taken from a camera in fixed position. The series must be taken under controlled laboratory conditions that ensure that the temperature and lighting remain constant. Ideally, any two images of the series should be exactly equal in absence of noise. Therefore, any variation between the images is only explained by stochastic light fluctuations (photon noise and shot noise) and the noise generated by the camera itself (dark noise, readout noise and electronic noise). In Section 1.2 of Chapter 3 it is explained in detail how to build the ground-truth curve for a particular camera and ISO speed. The ground-truth noise curves of the Canon EOS 30D and Nikon D80 for ISO speeds 1250 and 1600 are shown.

If $\tilde{\mathbf{U}}_i(\mathbf{x}, \mathbf{y})$ is a pixel of the noisy image $i$ at position $(\mathbf{x}, \mathbf{y})$, the intensity of the ideal image can be approximated by its empirical expectation $\hat{\mu}(\mathbf{x}, \mathbf{y}) = \mathbf{E}[\{\tilde{\mathbf{U}}_i(\mathbf{x}, \mathbf{y})\}]$ for $i = 1, \dots, N$, where $N$ is the number of snapshots in the series. The empirical variance associated to intensity $\hat{\mu}(\mathbf{x}, \mathbf{y})$ is $\hat{\sigma}^2(\mathbf{x}, \mathbf{y}) = \text{Var}[\{\tilde{\mathbf{U}}_i(\mathbf{x}, \mathbf{y})\]}$.

The calibration pattern must be mostly flat and represent a wide range of gray levels. Since the noise curve mainly depends on the ISO sensitivity, a different noise curve is estimated for each ISO level. Series of different exposure times were taken for each ISO in order to get representative information in the whole gray level range. The noise curves for different times of exposure were combined to obtain a single curve. In order to get a ground-truth noise curve, for each exposure time (1/30s, 1/250s, 1/400s, 1/640s) about two hundred pictures of the calibration pattern were taken. Since each 2 × 2 block of the CFA (Color Filter Array) contains one sample of the red channel, two samples of the green channel and one sample of the blue channel, the raw image was resampled as an image with four different color channels of half width and height. Thus, four different noise ground-truth curves were obtained from the series, each corresponding to one of the four channels of the CFA. By splitting the color range into bins (disjoint in intensity intervals) and
computing the median value of the STDs at each bin, a ground-truth is obtained for the camera noise curve given the ISO and exposure times.

Figure 1 shows the noise curves obtained with a Nikon D80 camera with fixed ISO sensitivity of 1250 and 1600 and four exposure times, $t \in \{1/30s, 1/250s, 1/400s, 1/640s\}$. The obtained curves overlap perfectly. Each one treats a different color interval, thus permitting to fuse them into a single noise curve. This fused curve can be observed in the same figure. For each color value, the fused noise estimation is obtained by the median of the available estimations obtained for the different exposure time. The value for each curve is linearly interpolated using the two closest neighbors. Since the noise curve does not depend on the exposure time, these curves overlap (hence the double values). However, this overlap is not perfect because the STD is computed with a finite number of samples and therefore the estimation has some variance that causes a small error centered at the theoretical value. Curve (b) is the mean of all four curves at different exposure times, which cancels their variation around the theoretical value and therefore it can be used finally as a ground-truth for evaluating noise estimation algorithms. Figure 5 displays the approximation of the computed ground-truth values by a linear model with the Nikon D80 camera. Because of the saturation at the darkest zones, the estimated noise in the dark gray level does not follow a linear model. However, using the partial linear model splitting the curve into three parts might be useful to model this kind of curves is the noise model is known in advance. The ground-truth curve obtained with the procedure presented here describes accurately the characteristics of the noise without depending only on the minimal assumption that the noise depends only on the intensity. Parametric methods assume priors for a particular noise model and therefore their accuracy depends on how realistic these assumptions are. Section 4 shows that it is possible to get a reliable noise curve that matches with negligible error the non-parametric ground-truth and Section 5 shows that indeed it is possible to validate parametric methods with the non-parametric ground-truth curve.

## 4. Non-parametric signal-dependent noise estimation

Parametric models are accurate under the condition of prior knowledge about the noise model. For example, the Foi et al. [7] method assumes the linear model $\sigma^2 = A+BU$ for the variance, but with a saturation effect. On the other hand, Uss et al. showed that the measured noise variance cannot always be fitted with a linear function, but with a polynomial of at least second order [75]. However, we were unable to fit a 2nd, 3rd or 4th order polynomial to the saturated noise curve in Figure 5. In order to use a linear function, these authors modify the intensity of the pixels at each SW of the image by a function that nullifies quadratic and higher terms of the noise variance model. After this transformation, the estimation is accurate.

Parametric methods require a validation, by a comparison to ground truth noise curves. The data in the ground-truth must be empirical, in the sense that it does not assume any prior (with the exception that the variance of the noise is a function of the expectation) and simply measures the variance of the noise as-is. As discussed in Section 3, the major problem of the comparison
<figure><img src="image_20.png" /><figcaption>FIGURE 1. Noise ground-truth curves obtained for a Nikon D80 camera with fixed ISOs of 1250 (a) and 1600 (c) and four exposure times, $t \in \{1/30s, 1/250s, 1/400s, 1/640s\}$ using laboratory calibration. Channels $G_1$ and $G_2$ give the same STD. The obtained curves overlap perfectly. Since they cover different color intervals, their fusion yields a complete noise curve (b), (d).</figcaption></figure>

against the ground-truth is that it is different for each camera model and it must be obtained under controlled laboratory conditions.

Our goal here is precisely to show that the laboratory calibration method used to obtain the ground truth can be replaced by a non-parametric method, estimating directly on the image the signal-dependent noise. We adapted the Ponomarenko et al. method [16], since it is scored as the best method in a previous review [1]. Other non-parametric estimation methods could be used as well. For example, in the paper of Liu et al. [52] a flexible eigenfunction representation of the noise level curves was proposed, but it needs a prior segmentation of the noisy image.

We extended the Ponomarenko et al. method [16] to be able to estimate signal-dependent noise (see the details in Chapter 10). To this aim, the means of the blocks are classified into a disjoint union of variable intervals (bins), in such a way that each interval contains a fixed and large enough number of elements. Thus, these intervals are automatically adapted to the image itself instead of prefixed, since the intensity range of each bin depends on the mean intensity of the $w \times w$ blocks in the image. We found that each bin should contain at least 42000 samples,
which seems to be the lowest number permitting a reliable estimation. The value of $p$ for the $p$-quantile of block variances must be small to avoid blocks with large variance, corresponding to edges and textures. In general, if a bin is made of blocks that belong to a flat or smooth zone, we found experimentally that 210 per bin are enough to estimate the variance (using $p = 0.5$, the median). However, with a smaller $p = 0.005$ percentile value, we can discard 99.5% of the blocks with a higher variance and therefore in general all of the blocks affected by edges and textures. By choosing a large value for the bin cardinality, namely 42000, we ensure that the blocks below this low quantile are still numerous enough: $42000 \times 0.005 = 210$, so that they ensure a reliable estimate of the noise variance.

To each bin a list of image blocks is associated, each of them being endowed with a list of STDs. Notice that a bin does not correspond to a spatial region of the image, but only to a set of blocks with similar means.

Another modification with respect to the original method is the procedure to find the best $p$-quantile. The values in the list $\{V^H(i, j)\}$ depend on the choice of $p$-quantile. If $p$ is small, the method becomes more robust to the influence of textures and geometry of the image; but the accuracy of the estimation also decreases with $p$. Our assumption is that the variance measured using $p = 0.5$ (the median) should not be significantly different to the variance obtained with lower values of $p$, unless the image is composed mainly of textures. The proposed iteration to get a robust estimation of the variance, adapting $p$, is as follows. At each bin,

(1) Initially, $p = 0.5$ (the median) and $\Delta p = 0.005$.

(2) Set $S = \{V^H(i, j) | i + j \ge T\}$.

(3) Set $V_s$ to the median of the values of $S$ under the $\Delta p$-quantile.

(4) Set $V_p$ to the median of the values of $S$ under the $p$-quantile.

(5) If $p \ge \Delta p$ and $V_p \le V_s$ then [Set $p = p - \Delta p$ and go to step 3] else END.

This procedure decreases the initial $p$ from the median to a lower value that makes the estimation robust to textures and geometry, if needed.

About the $w \times w$ size of the scanning window, we use the same value that the authors of the original algorithm proposed, and we found that indeed the best results are obtained with $w = 8$ in most natural images. Since the optimal size of the window depends on the density of edges and level of texturization of the noisy image, if it is a priori known that the image is mainly composed by large flat or smooth areas, it is better to use a larger window (up to $21 \times 21$) and to choose a smaller size in the opposite case (but at least $3 \times 3$) to obtain a reliable variance noise estimation. A larger window estimates the noise more accurately (since the sample variance estimator has itself a variance that depends on the number of available samples), but is less likely to contain only data from flat or smooth zones, and more likely to capture edges and geometry. Nevertheless, the proposed method is “blind”, in the sense that no prior information about the characteristics of the image or the noise is available. Adapting the window size is beyond the scope of this chapter and it is left as future work.
To avoid outliers in the estimation, we systematically discard completely saturated blocks. Indeed, when the number of photons counted 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. Indeed, the effect of the saturation must be measured and given by the non-parametric method in the produced noise curve, but the completely saturated pixels have outlier intensities. Figure 2 shows a noise curve obtained by using or avoiding the saturated blocks, where the modified Ponomarenko et al. [16] algorithm was performed with 49 bins. Since the intensity of the saturated pixels is much higher (outlier) than the intensity of the rest of the points, the noise curve is linearly interpolated along the gap in between. This natural image is a normal scene that is useful to illustrate the problem of the saturation. The bike is not illuminated directly by any source of light (only ambient light) and therefore it does not reflect much light, with the exception of a few points at the handlebar that reflect light with enough power to saturate the detectors. The STD equal to zero (measured near intensity 4000) is indeed correct, but all the interpolated points in between are definitely not. The strategy we adopted was to discard the blocks that contain a sub-group of 2 × 2 pixels sharing the same intensity, in any of the channels. It must be noted that this only removes the blocks containing pixels that are completely saturated, but keeps the rest of the blocks, including those where the noise distribution is truncated, but not absolutely saturated. This permits to measure and observe the saturation in the curve, as shown in Figure 5.

5. Cross-validation of several methods. Discussion

In order to compare ground-truth, parametric methods and our non-parametric method, we used a dataset of 20 images obtained with a Nikon D80 camera using ISO 1250 and exposure time 1/640s. In these images the darkest pixels are saturated and therefore the noise curve does not follow the linear variance model. Our dataset contains some views of a room with objects over a table, images of corridors, bookshelves, (Figure 5) stairs, and classrooms inside a building, with different lighting levels. Also, two outdoor images of highly textured images (Figure 4). For each test image, we computed the RMSE between the STDs given by the method and by the ground-truth. The control points are given by the method and the STD of the ground-truth corresponding to that intensity is obtained by linear interpolation between the two nearest intensity control points of the ground-truth curve.

Figure 3 shows the obtained results. In general, the RMSE of the modified Ponomarenko et al. (red curve) method is close to zero, which means that it could be used to establish a (non-parametric) ground truth. The estimations given by Foi et al. (green curve) and Uss et al. (blue curve) are really close to the non-parametric ground-truth, and are therefore also validated by our approach. The Foi et al. method failed to measure the noise correctly when the images are composed mainly by textures (images #19 and #20, see Figure 4), whereas the Uss et al. and the proposed method gave good results in that case. Note that textures cause a localized error in the non-parametric curve (middle), whereas they cause a global error in the parametric curves. Fig 5 shows three examples of images in our dataset (images #8 and #12) where all algorithms
<figure><img src="image_21.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(c)</figcaption></figure>

FIGURE 2. Noise curve obtained when the saturated pixels are avoided in the noise estimation (b) and when they are taken into account (c) by using the Ponomarenko et al. method [16] with 49 bins in an image with saturated pixels (a).

estimated the noise correctly. Foi et al. method (red curve) matches accurately the ground-truth curves (green and blue), since it is designed to predict the shape of the curve under saturation conditions, whereas Uss et al. estimation is overall correct, except in the saturation zone, as expected. As explained in Section 2, the original Ponomarenko method is only able to estimate homoscedastic noise, that is, a value of STD that does not depend on the intensity. However, in Figs. 4 and 5 we show noise curves, that correspond to the modified Ponomarenko method: added bins to get control points in the curve for different intensities and avoid using completely saturated points before the estimation. Of course, the over-estimation caused by a bin where all samples belong to textures can be avoided if more than a single image is available, by estimating the noise in the mosaic made of several different input images.

As shown in Figure 5, the linear model does not hold when the image is saturated. Uss et al. tried to use a second-order polynomial to fit the saturated noise curve. However, we found that a second-order polynomial was not general enough to fit the saturated curves. Foi et al. assumed the linear model, but taking into account the effect of the saturation. However, both methods assume that the noise can be modeled with a linear function when there is no saturation. This is true for most CCDs or CMOS detectors, but the output recorded in the raw file given by the camera might not be a linear function of the intensity [77]. This makes clear the necessity of validating parametric methods, which assume an *a priori* model for the noise. In contrast, the
<figure><img src="image_24.png" /><figcaption>FIGURE 3. RMSE between the methods and the ground-truth for all the 20 images in our dataset. In general, the RMSE of the modified Ponomarenko et al. (red curve) method is close to zero, which means that indeed it can be considered to be a non-parametric ground-truth curve. The estimations given by Foi et al. (green curve) and Uss et al. (blue curve) are really close to the non-parametric ground-truth, and therefore they are also validated by our approach. (a): obtained RMSEs/image, (b): detailed view.</figcaption></figure>

estimates of a non-parametric method rely on the minimal assumption that the signal is a function of the expectation. In general, the best results are obtained with the modified Ponomarenko et al. method.

To decide if a method is valid or not, its RMSE with respect to the non-parametric ground-truth has to be compared to a threshold. We consider that a method is valid if the RMSE between the measured STD and the ground-truth is less or equal to $\Delta\hat{\sigma}_8 = 0.15$ (assuming that the images are encoded with 8 bits. This value was chosen to be as low as possible and, at the same time, consistent with the accuracy of state-of-the-art noise estimation methods. Since the raw images are encoded with 12 bits, the threshold is $\gamma = \Delta\hat{\sigma}_8 \times 16 = 2.4$. The estimation of the Foi et. al method is considered valid in 17 of 20 images whereas the Uss et al. method is validated with all the images in our dataset.

**5.1. Complexity.** The Uss et al. method follows four steps: (1) initialization of the TI map and the polynomial function for the variance, (2) estimate texture and noise variance for each TI and NI SWs and label the SW into NI or TI, (3) update the CRLB and finally, (4) apply the noise estimator to the samples associated to each bin and update the variance polynomial. Steps 2, 3 and 4 are iterated until convergence is reached. The complexity of the Uss and the modified Ponomarenko method is similar and their complexity is linear with the number of pixels in the image. Both imply an estimation of the noise variance at the DCT coefficients in small patches of the image after classifying the them according to their intensity. For its part, the Foi et al. method
<figure><img src="image_25.png" /><figcaption>FIGURE 4. highly textured images that caused small oscillations in the noise curves with the proposed method and wrong results with Foi et al: images #19 (a) and #20 (b) (see the obtained RMSEs in Figure 3). (c), (d): the ground-truth obtained with the series (green), the non-parametric ground-truth (darker blue), the Uss et al. method (brighter blue) and the Foi et al. method (red).</figcaption></figure>

follows these steps in order to obtain the final parametric model: compute the detail wavelet coefficients of the image, segment the image to find homogeneous zones, estimate locally pairs of intensity/variance, and finally the maximum-likelihood fitting of the global parametric model. All steps can be computed quickly, but unlike Uss. and the modified Ponomarenko method, it requires a previous segmentation of the image.

**5.2. Denoising results.** We used the noise curves obtained with the Uss, Foi, and the modified Ponomarenko methods as the input of the NL-Bayes [78] denoising algorithm, after applying a VST to the noisy image. Only the green channel was used. Note that according to our threshold criterion, both the Uss and Foi methods are validated and therefore their denoising results in almost all images in the dataset are very similar. Figure 6 (a) shows details of the results obtained for the image #3 of our dataset, where the Foi et al. method failed to estimate correctly the noise. While Uss and the modified Ponomarenko methods denoise the image properly, the noise at the dark zone (the bag over the table) remains visible. The image (e) is the test image #20 of our dataset, where the Uss and the modified Ponomarenko method give an valid estimation, whereas the Foi method overestimates. All methods gave an increased RMSE for that particular image,
<figure><img src="image_26.png" /><figcaption>FIGURE 5. Examples of images (#8 (a) and #12 (b) in our dataset) where all algorithms estimated the noise correctly. Their noise curves along all the intensity range (c), (d). Detail of the noise curves only within the range of the estimation given by the modified Ponomarenko et al. method (non-parametric ground-truth, green curve), (e), (f). Note that the Foi et al. method (red curve) matches accurately the ground-truth curves (green and blue), since it is designed to predict the shape of the curve under saturation conditions, whereas the Uss et al. estimation is overall correct, except in the saturation zone.</figcaption></figure>
which causes blurred denoised images and loss of fine details. Both the Uss and the modified Ponomarenko methods give similar visual results, whereas the overestimation in the method blurs the image even more.

<figure><img src="image_27.png" /><figcaption>FIGURE 6. Details of the denoising results with the NL-Bayes algorithm using the noise curves obtained from the noisy images (a, e) with the modified Ponomarenko (b, f), Uss et al. (c, g), and Foi et al. (d, h) methods. Image (a) is the test image #3 of our dataset (it is very dark, so we increased brightness for visualization purposes), where the Foi method was unable to give a reliable estimation and thus the noise is not removed at the dark zones and remains visible (the bag over the table). Image (e) is a detail from the test image #20 of our dataset, where the Uss and the modified Ponomarenko method give an valid estimation, whereas the Foi method overestimates.</figcaption></figure>

## 6. Conclusion

We showed that estimating an accurate noise curve from a single raw image is possible and can be done by an adaptation of a non-parametric noise estimator [16]. The only minimal assumption is that the noise STD is a function of the expected signal. Being able to apply a noise estimator (with relatively low complexity compared to denoising algorithms) to each raw image frees the users of a tedious and sometimes impossible camera calibration task. Indeed, noise curves obtained in an optical lab require measurements for each ISO and each optical setup, a heavy and costly procedure. By estimating the noise directly on the raw image, there is no risk of model error or accuracy loss caused by a noise parameter estimation on another camera.

According to the provided RMSE results (see Figure 3), the non-parametric method proposed here exhibits a very stable error (close to RMSE=0.5) when the image is not composed mainly of
textures. However, even if the image is highly textured (see images #19 and #20 in Figure 4), the error is small and similar to the RMSE obtained with the compared state-of-the-art methods.

In general, the estimation given by the proposed method is as reliable as the actual ground-truth obtained from the temporal series of the series of images of a calibration pattern in the laboratory, and matches the best parametric methods.

Signal-dependent noise is a sufficient model when estimating noise in raw images, but this model is not enough, in general. Noise is transformed at each step of the camera processing chain (raw image, demosaicing, white balance, gamma correction, and JPEG encoding). After demo-saicing, the noise becomes correlated (and thus, frequency-dependent), as explained in Chapter 3. Chapter 4 will propose an algorithm for estimating both intensity and frequency dependent noise, which is valid even for JPEG-encoded images.
# The noise throughout the camera processing pipeline

Noise in raw images follow a signal-dependent Poisson model and can be estimated by adapting the state of the art homoscedastic noise estimation algorithms [9, 16]. However, the JPEG standard [79] is preferred in general, since it reduces drastically the size of the file by using lossy and lossless compression. In the raw image the noise follows a Poisson distribution, it is white, and uncorrelated. After JPEG compression, the noise is frequency-dependent, correlated and no longer white. We will show how the noise is affected at each step of the camera processing chain. The noise curves obtained with the Ponomarenko et al. method [15, 16] along the complete camera processing chain (raw image, demosaicing, white balance, gamma correction, and JPEG-encoding) are shown (along the corresponding autocorrelation matrices of the noise) and compared to the temporal estimation (nonparametric ground-truth curve, see Section 2 of Chapter 3).

## 1. The noise curves at each step of the camera pipeline

In this chapter, the effects of the camera processing chain (from the raw to the final JPEG image) on the noise will be discussed, showing detailed noise curves at each step, using different cameras, ISO speeds, and exposure times. Finally, it will be shown that it is possible to obtain a unique noise curve for raw images depending only on the camera model and the ISO that was configured when the picture was taken. For all the experiments in this chapter, two different cameras (Canon EOS 30D and Nikon D80) were used, combining two ISO speeds (ISO 1250, ISO 1600) and four exposure times (1/30s, 1/250s, 1/400s, 1/640s).

### 1.1. Noise at the Color Filter Array.
Each cell of the mosaic acquired by the CCD<fnref n="1" /> or CMOS<fnref n="2" /> detectors of most digital cameras presents a value at each pixel that can be modeled as a Poisson variable whose expectation is the actual “true” color [8]. The random fluctuations of this value around its mean can be considered independent and can be estimated accurately with state-of-the-art homoscedastic noise estimation algorithms [16, 1, 3]. More precisely, the noise measured at the CCD or CMOS detector is the combination of several sources:

* Dark noise, caused by a small electronic current known as the dark current, created at the depletion zone of the semiconductors of the CCD or CMOS detector because of the high intensity electric fields. The dark current (and therefore the dark noise) increases with the temperature. The dark current at the same temperature is a function of each captor at the CCD or CMOS detector, and therefore the additional current does not

<footnote n="1">¹Charge-Coupled Device.</footnote>

<footnote n="2">²Complementary Metal-Oxide Semiconductor.</footnote>
vary between snapshots. This pattern observed by the camera when the CCD or CMOS
detector is not exposed to light is known as the *dark frame*. Since the dark frame pattern
is always the same and nothing more than a constant bias, it can be subtracted after
taking the photograph of a scene.

*   **Photon noise**, due to the physical nature of light. The photons are emitted as quantum of energy with a rate that has some variance. As we mentioned, this process can be modeled with a Poisson distribution.

*   **Readout noise.** Since the charge accumulated at the CCD is really small, it must be processed by an analogic amplifier which adds noise to the measurement.

*   **Shot noise.** The noise of the Poisson photon acquisition can be modeled in good light as an additive Gaussian noise with standard deviation equal to the square root of the expected intensity. Nevertheless, this approximation is no longer true for very low photon income, in which case the Poissonian image is called “shot noise”. Indeed, in a dark scene most of the photosensitive surface receives a few photons whereas the rest of the surface does not. Therefore, after the amplification of the measured charge, it causes that isolated bright dots appear in the output image.

*   **Electronic noise** caused by the absorption of electromagnetic energy by the semiconductors of the camera circuits and the crosstalk phenomenon, among others.

**1.2. The camera pipeline.** Most of the noise evaluation methods assume or require the noise to be uncorrelated and thus they are adequate to estimate directly the noise in the image formed at the CFA<fnref n="3" />, the raw image. Unfortunately, the first transformation applied to the CFA to obtain a color image (by demosaicing the Bayer pattern) correlates the noise. In addition, the lossy compression performed by JPEG encoding makes the noise frequency-dependent. Many methods [16, 1] assume that the variance of the noise can be estimated by measuring the variance only at high-frequency coefficients, supposing that noise is frequency-independent. Yet, JPEG compression quantizes the values of the coefficients depending of the frequency. This energy loss of certain coefficients causes an underestimation of the noise estimation in these methods.

The following steps and transformations are performed in the camera processing pipeline in
order to obtain a final JPEG image from the raw:

(1) Acquire the image data from the CFA (raw image);

(2) interpolate the data from the CFA to obtain a color image from the grayscale mosaic (demosaicing);

(3) adjust the weight of each color channel in order that the picture has the same colors as the photographed scene (white balance);

(4) modify the dynamic range of the image in order to give more importance to the intensities that human observers are more adapted to detect. This is done a non-linear function to

<footnote n="3">³Color Filter Array.</footnote>
<figure><img src="image_28.png" /><figcaption>FIGURE 1. Left: calibration pattern used to obtain the ground-truth of the cameras. Right: photograph of the calibration pattern taken by the Canon EOS 30D camera, raw image. The image is not grayscale because no white balance has been applied at this first step.</figcaption></figure>

each pixel, called **gamma correction**. For some cameras this function cannot be formulated with a simple expression and it is given as a table that maps the input intensities into different values (**tone curve**);

(5) apply lossy compression (quantization of the DCT coefficients).

(6) apply lossless compression (Huffman encoding).

Another objective of this chapter is to build *ground-truth* (GT) noise curve for the noise at the CFA for both the Canon and Nikon cameras, depending on the ISO. A noise curve associates with each observable intensity in each color channel a value for its standard deviation. To evaluate and validate different noise estimation algorithms the GT is an absolute requirement. To build the ground-truth curves, the calibration pattern shown in Figure 1 (left) was printed with a high quality plotter and several (about 500) raw image pictures of the calibration pattern were taken with the Canon and Nikon cameras with the same lighting conditions, fixing the ISO and the exposure times. To minimize the effect of the edges in the image causing fluctuations due to sub-pixel displacements of the still camera, and not to noise, the photographed calibration pattern contains large different gray level rectangles (Figure 1, right). The image has a green tone because the Bayer pattern presents two green pixels for one red and one blue. The variance of a pixel value along different snapshots of the same still scene (Figure 2) can only be explained by the noise and therefore the noise. Thus, the noise curve obtained computing the standard deviation values of the temporal series gives the ground-truth noise curve for that camera. This noise curve obviously also depends on the ISO and on the exposure time (Figure 4). As shown in Section 2, since the standard deviation of the noise is a function of the intensity, it is possible to overlap all the noise curves from different exposure times into a single ground-truth curve combining the curves from all the exposure times.

1.2.1. *Step 1/5: raw image.* The first step to acquire the raw image is to count the number of incident photons over the CFA along the exposure time, using a CCD or CMOS detector. Because of the photoelectric effect, the CCD or CMOS detector is able to accumulate electric charge by the absorption of electromagnetic energy. The camera electronics measures the voltage produced by the accumulated charge during the exposure time, digitizes the values and finally stores them using some proprietary format (NEF for the Canon and CR2 for the Nikon cameras). For both
<figure><img src="image_29.png" /><figcaption>FIGURE 2. Computing the STD of the temporal series with the calibration pattern.</figcaption></figure>

cameras, the intensity of a pixel in one of the four channels (*R*, *G*₁, *G*₂, *B*) is encoded with 12 bits. At this first step, the noise is uncorrelated, that is, the noise at a certain pixel is not related with the noise at any other pixel with the same signal intensity. Figure 3 shows the autocorrelation function of a rectangle in the calibration pattern. This function only presents a peak at [0, 0], that is, the pixels are almost only correlated with themselves and not with any other neighbor pixel. It can also be observed, for the Canon camera (left), a pattern of some pixels at distance 2 and beyond with small autocorrelation, caused by crosstalk between captors of the same type in the CFA.

<figure><img src="image_30.png" /><figcaption>FIGURE 3. Autocorrelation function for a rectangle of the calibration pattern for the Canon (left) and Nikon (right) cameras with ISO 1600 and t=1/250s in the raw image. The noise is mostly not correlated, although it can be observed, for the Canon camera, a pattern of some pixels at distance 2 and beyond with small autocorrelation, caused by crosstalk between captors of the same type in the CFA. On the other hand, the autocorrelation at [0,0] shows that the image obtained with the Nikon camera contains more noise than the image of the Canon camera.</figcaption></figure>

The noise curve obtained by the temporal series (ground-truth) and the noise curve obtained by any accurate noise estimation algorithm coincide, since the noise is Poisson distributed. Figure 4 shows how the curves match in the case of the Ponomarenko et al. method [16], using ISO 1600 and t=1/250s for both the Canon and Nikon cameras.
<figure><img src="image_31.png" /><figcaption>FIGURE 4. Comparison between the temporal series standard deviation (ground-truth) and the noise curve obtained by the Ponomarenko et al. algorithm in raw images. 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. Left: Canon camera. Right: Nikon camera.</figcaption></figure>

Since noise is Poisson distributed, the variance of the noise is directly related to the intensity of the underlying signal. Therefore, the variance of the noise should follow a simple $σ^2 = a + bu$ linear relation, where *u* is the intensity of the ideal noiseless image and *a* *a* is due to the other mentioned sources of uniform white noise explained in Section 1.2, dark, readout, and electronic noise. In the case of raw images, the noise is only a function of ISO, since it controls the intensity multipliers. Therefore, it is possible to combine the curves obtained with different exposure times (and the same ISO) to get a unique ground-truth curve.

However, some cameras pre-process the data acquired at the CCD before writing the raw image. Also, the saturation of the captors at the most dark and bright zones breaks the linearity between the intensity and the variance. Figure 5 shows the square root of the linear model of the variance obtained by the least squares method. In the case of the Nikon camera, the image is dark enough to show the non-linear effect of the saturation in the noise curve. It is therefore preferable and more reliable to estimate the noise curve without assuming *a priori* a noise model. Thus we shall prefer non-parametric estimation methods that give a signal-dependent noise curve from the image itself.

### 1.2.2. Step 2/5: demosaicing.
The image acquired at the CFA is just a mosaic where each captor measures the intensity of the *R*, *G*₁, *G*₂ or *B* channels. Depending on the CFA, the disposition of the pixels associated with each channel may be different. The goal of the demosaicing step is to create a color image from the mosaic using an interpolation algorithm. In our tests we used the Adams-Hamilton [80, 81] because of its simplicity. Indeed, it bases its estimation on the gradients and pixels at distance one and two.
<figure><img src="image_32.png" /><figcaption>FIGURE 5. Linear approximation of the variance ground-truth for the Canon (up) and the Nikon (down) cameras with ISO 1600 (solid line: original values, dashed line: linear approximation). Globally, the linear approximation coincides with the computed ground-truth values. On the right, a detailed view of the estimation for the darkest pixels in the same noise curve. For the Nikon camera, the estimated noise in the dark gray level values does not coincide with the linear approximation, because of saturation.</figcaption></figure>

Figure 6 shows the comparison between the temporal series and the estimation obtained using the Ponomarenko et al. algorithm for both the Canon (left) and Nikon (right) cameras.

It shows that

* each channel has a different noise level.

* the temporal series has more noise than the single image estimation.

Now each channel has a different noise level since the Adams-Hamilton algorithm does not process each channel in the same way. The temporal series and the single image estimation coincide in the raw image but not after the demosaicing. The reason is that the Adams-Hamilton demosaicing does, among other operations, an averaging of pixels at distance one, and therefore reduces the spatial standard deviation that is measured in a single image. However, in the temporal series the standard deviation is computed using pixels at the same location in the temporal series. This explains why the noise in the temporal series is higher than the noise measured using a single image. Since the image processing pipeline is sequential, the temporal noise curves and those measured on a single image will not coincide anymore after the demosaicing step. Figure 6 permits to compare the estimation obtained using the temporal series and using the Ponomarenko et al. algorithm.
<figure><img src="image_33.png" /><figcaption>FIGURE 6. Comparison between the temporal series standard deviation and the noise curve obtained by the Ponomarenko et al. algorithm after demosaicing. 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 curves do not coincide anymore. Left: Canon camera. Right: Nikon camera.</figcaption></figure>

In the raw image, the autocorrelation is close to zero for any location different from [0, 0], meaning that the noise is not spatially correlated. In the demosaiced image the autocorrelation goes down to zero at distance two or more. This means that the pixel noise is correlated with the noise at pixels whose distance is one. Figure 7 shows that the maximum autocorrelation is attained at [0, 0], but the four points at distance 1, {[0, -1], [0, 1], [-1, 0], [1, 0]}, also have a significant autocorrelation. This result is coherent with the Adams-Hamilton algorithm, where pixels at distance higher than one show a very low autocorrelation and can be considered non-correlated in practice.

<figure><img src="image_34.png" /><figcaption>FIGURE 7. Autocorrelation function for a rectangle of the calibration pattern for the Canon (left) and Nikon (right) cameras with ISO 1600 and t=1/250s in the demosaiced image. Each pixel is correlated with four pixels at distance one. Pixels at a distance higher than one show a very low autocorrelation and they can be considered to be non-correlated in practice. Also, the Nikon camera shows a higher autocorrelation because the noise is higher.</figcaption></figure>
1.2.3. *Step 3/5: white balance.* After the demosaicing each pixel of the image has a color, but it has to be corrected in order to show the same colors a human would observe. The reason is that the color filters of the CFA do not have all the same gain. Therefore, taking a picture of a gray object will produce a colored image if the gain for the R (red), G (green) and B (blue) channels is not balanced. For example, the photograph of the calibration pattern in Figure 1 (right) is bluish since the blue channel has more gain than the others. However, apart from correcting the colors, the white balance allows to simulate the kind of light on a scene. The kind of light is well defined using its *color temperature*, that is, the kind of light that would emit a blackbody at that temperature. For example, a candle light is within the range [1000K, 2000K], the fluorescent lamps within [4000K, 5000K], the light of a sunny day within [5000K, 6500K], etc. Because the calibration pattern only contains gray rectangles, the multipliers are chosen in such a way that the resulting rectangles are gray, that is, the same intensity at each color channel. For example, fixing ISO 1600 and exposure time 1/250s, the multipliers for the Canon camera are $M_R = 1.30$, $M_G = 1.00$, $M_B = 1.08$ and for the Nikon they are $M_R = 1.8$, $M_G = 1.00$, $M_B = 1.14$. In both cameras, the green channel is more sensitive that the red or blue channels and therefore the white balance step increases them more. Figure 8 shows the results of the noise estimation after the white balance for both Canon and Nikon cameras, with ISO 1600 and exposure time 1/250s. Since the white balance factors are higher than one, the whole dynamic range of the image is increased. Compare the maximum intensities of the noise curves corresponding to the demosaicing step and the white balance step. The white balance of course increases the noise. For example, the factor $M_R$ is higher than $M_G$ and $M_B$ as it can be seen in the curves. The red is now the most noisy channel (then the blue and finally the green channel, which was not modified by the white balance). Because of the dark, readout, and electronic noises, the minimum intensity registered by the camera is not encoded at zero, but at some fixed value. It can be observed that this value in the Canon camera (near intensity 128) is higher than in the Nikon camera (near intensity 30) and therefore multiplying the red channel by 1.30 is much more noticeable in the Canon curve than in the curve corresponding to the Nikon camera.

1.2.4. *Step 4/5: gamma correction (tone curve).* Human vision perception is not linear with the signal intensity and can be modeled approximately with a power function. However, the charge accumulated by the CCD or CMOS detector is linear with the number of incident photons on the device during the exposure time. Since the information at the darkest zones is invisible to a human observer, it is applied a power function called the *gamma correction* to the linear data captured by the CCD or CMOS detector before attempting lossy compression. The idea is not to only enhance the contrast of the image but also to encode more accurately the information in the dark areas that is invisible in the raw image by enhancing it with the concave gamma function. This gamma function has the form $f_{k,\gamma}(u) = ku^{\gamma}$, where the typical values for $\gamma$ vary from 1.8 to 2.2. However, commercial cameras do not apply this simple formula but use some precomputed tables to simulate the non-linearity, called *tone curves*. In general, the tone curves saturate the signal at the brightest zones of the image. Figure 9 shows the effect of applying the gamma correction
<figure><img src="image_35.png" /><figcaption>FIGURE 8. Comparison between the temporal series standard deviation and the noise curve obtained by the Ponomarenko et al. algorithm after white balance. 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 minimum intensity value is much lower for the Canon than for the Nikon. Therefore, the displacement to the right of the red channel is more noticeable for the Canon camera. The dynamic range is increased and since each channel has different multipliers, the noise now is higher for the red channel, then for the blue and finally for the green, that has not been modified. Left: Canon camera. Right: Nikon camera.</figcaption></figure>

to the images obtained by the Canon and Nikon cameras. Previously, the demosaicing and white balance were applied.

It can be observed that

* the dynamic range is clearly increased, much more than the white balance stretching, specially in the dark areas;

* the noise increases significantly because of the power law function;

* the noise curve function is no longer monotonically increasing. When the intensity is over some threshold the noise saturates the captor and the measured standard deviation decreases instead of increasing (see Figure 9, right). Indeed, the derivative of the gamma function, which is concave, is larger than 1 in the dark part and smaller than 1 in the bright part. It is easily checked that in a first approximation, the noise standard deviation is multiplied by the derivative of the gamma function. Indeed, assuming that the noise $n$ is small with respect to the signal $u$, we have by a first order asymptotic expansion $\gamma(u+n) \approx \gamma(u) + \gamma'(u)n$. Thus, after gamma correction the noise at level $u$ is approximately $\gamma'(u)n$.

1.2.5. *Step 5/5: JPEG compression.* The final step is to encode the image with the JPEG standard in order to reduce the size of the resulting file. The encoding is not done in the *RBG*