Estimate Variance of Additive White Gaussian Noise (AWGN) Given Multiple Realizations with Different Mean

Question

Given $ N $ images of the same scene, where each image is corrupted by additive white Gaussian noise of the same variance $ {\sigma}^{2} $. How can $ {\sigma}^{2} $ be estimated?

So we have an Image $ X \in \mathbb{R}^{m \times n} $ which we have few realizations noisy variations of it $ {\left\{ {X}_{i} \right\}}_{i = 1}^{N} $:

$$ {X}_{j} = X + {V}_{j} $$

Where $ V \in \mathbb{R}^{m \times n} $ is AWGN such that $ {V}_{j} \left( i , j \right) \sim \mathcal{N} \left( 0, {\sigma}^{2} \right). $

Ex., for $N=2$, let $y$ be the ground truth scene and $x_n$ be the noisy relizations, then $x_1=y+ν_1$ and $x_2=y+ν_2$ where $ ν_1, ν_2 \sim N(0,σ^2) $ and hence $x_1−x_2=ν_1−ν_2:=ν$ where $ν∼N(0,2σ^2)$. Therefor $σ^2$ can be estimated by calculating the variance of the difference of the two realizations, and then dividing by two.

Answer 1

Basically each pixel is a realization so all you need is to work in the 3rd dimension (Though you can also get better by using the Spatial Data).

Method 1

So the trick here is to use the multiple images to estimate the Mean (True value) of each pixel and then calculate the STD on all samples (numRows * numCols * numRealizations).

Assuming we have single channel image. So we pack all give images into a tensor tI with dimensions: numRows * numCols * numRealizations.

We calculate the mean per pixel (Averaging on the 3rd dimension) and then subtract each image from the calculated average image.
Then we're left with many realizations of the noise (Well, noise and the left over from the estimated mean error).

Method 2

Estimate the STD in the 3rd dimension per pixel. Them average all numRows * numPixels estimations.

Method 3

Just as method 2. But since the property which obeys to linear operations is the Variance calculate the variance along the 3rd dimension, average it over all pixels and then take the sqrt() of the average Variance.

Summary

Here is a graph showing the estimated noise STD as a function of the number of realizations:

As can be seen, for an image with many pixels (Say more than 20,000 or so) Method 3 becomes almost perfect.

The full code is available on my StackExchange Signal Processing Q61273 GitHub Repository (Look at the SignalProcessing\Q61273 folder)..

Remarks

1. Method 1 - Basically the Maximum Likelihood Estimator for the problem. It is biased as this is the property of the ML for the Variance / STD.
2. Method 2 - Works nice, really intuitive.
3. Method 3 - Basically the MLE with the non biased estimation for the variance. As can be seen. It is also the best in real world.

Yep, Method 3 is much better only because it divides the data by N while Method 1 divides by N - 1. But it has significance impact on performance.