Why my calculation of standard deviation of an image is different from the built-in function

Hello,

My code is shown below, and I opened the built-in function. I found it just calculates the square root or variance. So, the method is the same. I don't know why.

function ret = yStdDeva(im)
tic;
[tR, tC] = size(im);

mean1 = mean(im(:));

% Calculate square of difference between pixels and mean
sdpm = (double(im) - mean1).^2;
sum_sdpm = sum(sdpm(:));

% Calculate variance
imvar = (1/(tR*tC   ))*sum_sdpm;

% Calculate Standard deviation
ret = (imvar)^(1/2);
toc;

Answer 1

[Jan]

You simply use a wrong formula: You have to normalize by the number of elements minus 1. In addition there can be a difference between SQRT and ^0.5 due to the different numerical implementations.

d = double(im);
C = d(:) - mean(d(:));
S = sqrt(sum(C .* C, 1) ./ (numel(d) - 1)):

[ZhG]

Thanks Jan.

[ZhG]

Excuse me, but why it is normalzied by 1. Is this specific for image?

[ZhG]

Actually, it is the same after I had modified my program as your answer. I am confused. The result obtained by built-in function is only 9.6938. But a result of 72.0412 was obtained by my program.

[Jan]

I have converted the image to the type double to avoid saturation effects with integer types. So please try my code and not a program modified like my code.

Answer 2

[ZhG]

However, I found where the problem happens. If I calculate standard deviation of an image without converting it into double, it will generate a result of 9.6938. However, it will does return 72.0412 afte I transformed it to double.

Is there anyone know what the problem is?

[Iain]

The problem is wind-up. To see it, try these at command line:

 uint8(255) + uint8(1) - uint8(1) % Get it here
 uint16(255) + uint16(1) - uint16(1) % Don't get it here
 double(10^30) + double(1) - double(10^30) % Get it here
 double(10^30) - double(10^30) + double(1) % Don't get it here

[Jan]

Operations on integer types are affected of the saturaion:

uint8(250) + uint8(250) == uint8(255) !!

[Iain]

Wind-up is a type of saturation, and I use the terms almost interchangably.

Answer 3

[Image Analyst]

I don't see much difference - just a really slight difference (less than a thousandth of a gray level):

im=imread('cameraman.tif');

% Your way:
[tR, tC] = size(im);
mean1 = mean(im(:));
% Calculate square of difference between pixels and mean
sdpm = (double(im) - mean1).^2;
sum_sdpm = sum(sdpm(:));
% Calculate variance
imvar = (1/(tR*tC   ))*sum_sdpm;
% Calculate Standard deviation
ret = (imvar)^(1/2)

% Standard, built-in way:
std(double(im(:)))

In the command window:

ret =
          62.3412396872523
ans =
          62.3417153186086

Why are you wanting to do it yourself anyway? Why not just use std()?

[ZhG]

because std is slower that what i did. but what i did is not useful. 0:)

Answer 4

[ZhG]

thank you all, guys.