Why does imnoise add different noise than randn?

I don't understand what imnoise is doing differently than when you add noise to an image with randn. Here is my code

I create an 8bit image with 8bit pixel depth. Then add gaussian noise manually and with imnoise. I remove the pixel values from the noisy images and use fitdist to find sigma. I also have a control where I use imnoise on zeros. I'd expect the noise added from all methods to have the same mu and sigma.

What's happening here?

EDIT: for typos

clear;

n=256;
mu=0;
sigma=2;
var=sigma^2;

% 8bit image with 8bit pixel depth.
sample = randsample(n,n*n,true);
I = reshape(sample,n,n);

% add gaussian noise
A = I + sigma*randn(n,n)+mu;
B = imnoise(rescale(I),"gaussian",mu,var);
C = imnoise(zeros(n,n),"gaussian",mu,var);

% fit to normal distribution
A_dist = fitdist(reshape(A-I,n*n,1),'normal');
B_dist = fitdist(reshape(B-I,n*n,1),'normal');
B_dist_prime = fitdist(reshape(B,n*n,1),'normal');
C_dist = fitdist(reshape(C,n*n,1),'normal');

A_dist.sigma

ans = 1.9942

A_dist.mu

ans = 0.0029

B_dist.sigma

ans = 73.8421

B_dist.mu

ans = -128.2028

B_dist_prime.sigma

ans = 0.4665

B_dist_prime.mu

ans = 0.5010

C_dist.sigma

ans = 0.4564

C_dist.mu

ans = 0.4018

Answer 1

[DGM]

You're seeing two things.

First, your image is improperly-scaled for its class (1-256 for class 'double'), but the image you're feeding to imnoise() is forced to unit scale. All floating point images fed to imnoise() are assumed to be unit-scale, so B and C are also unit-scale. B-I will result in nonsense.

Second, the distributions are all appearing wonky because you're seeing severe truncation. Internally, we're working in unit-scale. You're specifying gaussian noise with a sigma of 2, so a lot of that noise pushes data outside that interval. In the case of a null image, at least half the pixels are pushed outside [0 1]. The last thing that happens is the image is clamped to unit-scale and then cast/rescaled back to the input class.

If we stick with unit-scale and use moderate variance, the results are more comparable. Consider the example.

n=256;
mu=0;
gauvar = 0.01;
sigma = sqrt(gauvar)

sigma = 0.1000

% an image which is scaled properly for its class
I = rand(n,n);

% add gaussian noise
noise = sigma*randn(n,n)+mu;
A = imclamp(I + noise); % the output is clamped!
B = imnoise(rescale(I),"gaussian",mu,gauvar);
C = imnoise(zeros(n,n),"gaussian",mu,gauvar); % truncation is heavily asymmetric

% fit to normal distribution
ndist = fitdist(reshape(noise,n*n,1),'normal');
A_dist = fitdist(reshape(A-I,n*n,1),'normal');
B_dist = fitdist(reshape(B-I,n*n,1),'normal');
C_dist = fitdist(reshape(C,n*n,1),'normal');

[ndist.sigma ndist.mu]

ans = 1×2

0.0995 -0.0004

[A_dist.sigma A_dist.mu]

ans = 1×2

0.0941 -0.0002

[B_dist.sigma B_dist.mu]

ans = 1×2

0.0943 0.0000

[C_dist.sigma C_dist.mu] % still heavily truncated on one side

ans = 1×2

0.0580 0.0395

If you want to see exactly what's going on, you can check toolbox/images/images/+images/+internal/algimnoise.m, but that's basically it.

[Jonathan Sullivan]

Thanks!