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).
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]
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.
Scikit Image has an estimate sigma function that works pretty well:
http://scikit-image.org/docs/dev/api/skimage.restoration.html#skimage.restoration.estimate_sigma
it also works with color images, you just need to set multichannel=True
and average_sigmas=True
:
import cv2
from skimage.restoration import estimate_sigma
def estimate_noise(image_path):
img = cv2.imread(image_path)
return estimate_sigma(img, multichannel=True, average_sigmas=True)
High numbers mean low noise.