function s = snr(noisydata, original)

[m,n] = size(noisydata);

mean_original = mean(mean(original));
tmp           = original - mean_original;
var_original  = sum(sum(tmp.*tmp));

noise      = noisydata - original;
mean_noise = mean(mean(noise));
tmp        = noise - mean_noise;
var_noise  = sum(sum(tmp.*tmp));

if var_noise == 0
    s = 999.99; %% INF. clean image
else
    s = 10 * log10(var_original / var_noise);
end

return