Noise Estimation / Noise Measurement in Image

Question

I want to estimate the noise in an image.

Let's assume the model of an Image + White Noise. Now I want to estimate the Noise Variance.

My method is to calculate the Local Variance (3*3 up to 21*21 Blocks) of the image and then find areas where the Local Variance is fairly constant (By calculating the Local Variance of the Local Variance Matrix). I assume those areas are "Flat" hence the Variance is almost "Pure" noise.

Yet I don't get constant results.

Is there a better way?

Thanks.

P.S. I can't assume anything about the Image but the independent noise (Which isn't true for real image yet let's assume it).

Answer 1

You can use the following method to estimate the noise variance (this implementation works for grayscale images only):

def estimate_noise(I):    H, W = I.shape    M = [[1, -2, 1],        [-2, 4, -2],        [1, -2, 1]]    sigma = np.sum(np.sum(np.absolute(convolve2d(I, M))))   sigma = sigma * math.sqrt(0.5 * math.pi) / (6 * (W-2) * (H-2))    return sigma 

Reference: J. Immerkær, “Fast Noise Variance Estimation”, Computer Vision and Image Understanding, Vol. 64, No. 2, pp. 300-302, Sep. 1996 [PDF]

Answer 2

The problem of characterizing signal from noise is not easy. From your question, a first try would be to characterize second order statistics: natural images are known to have pixel to pixel correlations that are -by definition- not present in white noise.

In Fourier space the correlation corresponds to the energy spectrum. It is known that for natural images, it decreases as 1/f^2 . To quantify noise, I would therefore recommend to compute the correlation coefficient of the spectrum of your image with both hypothesis (flat and 1/f^2), so that you extract the coefficient.

Some functions to start you up:

import numpy def get_grids(N_X, N_Y):     from numpy import mgrid     return mgrid[-1:1:1j*N_X, -1:1:1j*N_Y]  def frequency_radius(fx, fy):     R2 = fx**2 + fy**2     (N_X, N_Y) = fx.shape     R2[N_X/2, N_Y/2]= numpy.inf      return numpy.sqrt(R2)  def enveloppe_color(fx, fy, alpha=1.0):     # 0.0, 0.5, 1.0, 2.0 are resp. white, pink, red, brown noise     # (see http://en.wikipedia.org/wiki/1/f_noise )     # enveloppe     return 1. / frequency_radius(fx, fy)**alpha #  import scipy image = scipy.lena() N_X, N_Y = image.shape fx, fy = get_grids(N_X, N_Y) pink_spectrum = enveloppe_color(fx, fy)  from scipy.fftpack import fft2 power_spectrum = numpy.abs(fft2(image))**2 

I recommend this wonderful paper for more details.