In your expression:
variance = medfilt2(double(I_noise).^2, [N N], 'symmetric') - (medfilt2(double(I_noise), [N N],'symmetric')).^2;
both of the terms on the right hand side of the expression are identical. Subtracting one from the other will obviously give zero.
This is the same for example as the expression:
x = 2 - 2
which of course gives an answer of x = 0;
Also I think you are getting confused about median and average.