Given a white Gaussian noise signal Nσ on a sampling grid, its variance σ2can be estimated from a small w×w pixels sample. However, in natural images we observe U˜ = U+Nσ, 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×w samples of 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 aected by the edges and textures than a higher percentile (50%, for example). The 0.5%-percentile is empirically proven to be adequate for most natural, medical, and microscopy images. The Percentile method is adapted to deal with signal-dependent noise, which is realistic with the Poisson noise model obtained by a CCD or CMOS detector in a digital camera.
1. Introduction
The Percentile method [1] assumes the signal-dependent noise model. The output of the Percentile method is a noise curve, that is, a function that relates the intensity of the image with a noise standard deviation. The Percentile method tries to estimate the standard deviationσ 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 aected by the edges and textures than a higher percentile (50%, for example). However, it might happen that for a certain intensity interval all the samples belong to textures and edges. In that case, the variance measured is explained by the geometry of the image and not the by noise. In order to minimize that eect, the noise curves can be ltered as explained in Section 1.2 of Chapter 9.
For a review of several noise estimation methods we refer the reader to the work of Lebrun et al. [1] and Olsen [55].
197
2. Noise Estimation Method
2.1. Notation and Terminology. This chapter prepares the detailed description of the noise estimation algorithm given in Chapter 2.2 by xing its notation and terminology.
• U: the noiseless ideal image.
• U: the noisy input image.˜
• Nx,Ny: the size of U.˜ Nxand Ny are always odd. If they are not, the leftmost column or the bottom row are rst removed fromU.˜
• U˜c: the discrete noisy image U after cropping˜ (s−1)/2 columns and rows from each of the four sides ofU using the function CROP˜ (s−1)/2, wheresis odd. The details of this function are given in algorithm 28. The size of U˜c is therefore (Nx−(s−1))×(Ny− (s−1)) = (Nx−s+ 1)×(Ny−s+ 1).
• R: the operator used to obtain the discrete lter F with supports×s, with sodd.
• F: the discrete lter with supports×sused to high-pass the noisy image, obtained from operator R.
• F(x, y): the value of the discrete lter F at position(x, y), x∈[0, s−1], y∈[0, s−1].
• ∧: logical conjunction (and operator).
• ∨: logical disjunction (or operator).
• ¬: logical negation (not operator).
• U˜f: the cropped high-passed version of U:˜ U˜f :=CROPs−1
0otherwise extends F with zeros,
U˜z(x, y) :=
U˜(x, y), x∈[0, Nx−1]∧y∈[0, Ny−1]
0 otherwise extends U with zeros and the˜
CROPs−1 function removess−1columns or rows at each of the four boundaries of the ltered image to avoid boundary eects. Since the convolution ( ˜U∗F)(x, y) is dened for x∈ [0, Nx+s−2]∧y ∈ [0, Nx+s−2], the cropped imageU˜f(x, y) is dened for x∈[0, Nx−s]∧y ∈[0, Ny−s]. The Fast Fourier Transform (FFT) algorithm [123] is used to speed-up the computation of the convolution; the details are given in Chapter 2.2.1.
• U˜c(x, y): the gray-level intensity of the pixel ofU˜cat columnxand rowy,x∈[0, Nx−s], y∈[0, Ny−s].
• U˜f(x, y): the gray-level intensity of the pixel ofU˜fat columnxand rowy,x∈[0, Nx−s],
• IVn: the index of the element at the positionn in the list V after sorting the elements of V in ascending order.
• p: the (small)p-percentile.
• σˆ2: biased variance at thep-percentile of V.
• σ˜2: nal unbiased variance at thep-percentile of V.
2.2. The Algorithm.
2.2.1. Step 1: Pre-ltering the Input Noisy Image. In order to get rid of deterministic tenden-cies due to signal structure, the image is rst pre-ltered with a high-pass lter F implemented as a discrete stencil with support s×s. This stencil corresponds to a discretization of a certain operator R, typically a dierential operator or a waveform. Convolving the image with such a lter removes smooth variations inside the blocks, which increases the number of blocks where noise dominates and on which the variance estimate will be reliable. Mastis proposed a similar approach [57], where operator F writes as a simple subtraction of the average or median to each 7×7block. For the Percentile method a lter based on the DCT with support7×7is proposed.
Given ans×sblock in the imageU at position˜ (x, y), its orthonormal 2D DCT-II is
Figure 1. From left to right, the discrete stencils associated ∆, ∆∆and ∆∆∆
discrete operators, respectively.
Filter F is made by taking the normalized product of cosines that correspond to the highest frequency[s−1, s−1], that is, The lter F presented here was empirically proven to give the best results. However, other typical dierential operators can be used, like directional derivatives, the ∆ (Laplace) operator, or its iterations ∆∆, ∆∆∆, all implemented as discrete stencils. Figure 1 shows the discrete stencils associated to these operators.
The ltered imageU˜fis obtained by cropping the discrete convolution, CROPs−1
h
( ˜U∗F)(x, y)i . Note that this cropping operation avoids the boundary eects of the convolution. To speed-up the computation, the Fast Fourier Transform (FFT) algorithm is used:
(1) Consider the signal