% Adds noise to a real image 
% Computes the SNR(u0/u)=signal-to-noise-ratio 

% factor coefficient of the noise in u0(i,j)=u(i,j)+k*noise(i,j)
% small k (less noise) => larger SNR
% large k (more noise) => smaller SNR
% remove noise by maximizing the SNR

k=15; 

% read an image in BMP format
u=imread('Lena.bmp'); 
u=double(u);
[M N]=size(u);

% visualize the image u in Matlab
imagesc(u); axis image; axis off; colormap(gray);

% noise of mean 0 and variance 1, multiplied by a factor k
noise=randn(size(u));

% construct a noisy version u0(i,j) of the original image u(i,j)
u0=u+k*noise;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% compute SNR(u0/u) using the variances

% compute mean of u0

mean_u0=0.;

for i=1:M,
  for j=1:N,
  mean_u0=mean_u0+u0(i,j);
  end
end

mean_u0=mean_u0/(M*N);

% compute variance of u0

var_u0=0.;

for i=1:M,
  for j=1:N,
  var_u0=var_u0+(u0(i,j)-mean_u0)^2;
  end
end

var_u0=var_u0/(M*N);

% compute mean of (u0-u)

mean_u0_u=0.;

for i=1:M,
  for j=1:N,
  mean_u0_u=mean_u0_u+(u0(i,j)-u(i,j));
  end
end

mean_u0_u=mean_u0_u/(M*N);

% compute variance of u0-u

var_u0_u=0.;

for i=1:M,
  for j=1:N,
  var_u0_u=var_u0_u+(u0(i,j)-u(i,j)-mean_u0_u)^2;
  end
end

var_u0_u=var_u0_u/(M*N);

SNR=10*log10(var_u0/var_u0_u)

% open a new figure window
figure

% visualize the image in Matlab (re-scaled)
imagesc(u0); axis image; axis off; colormap(gray);

% To save the image into an image format, use "imwrite" but you have first to
% convert it into the uint8 format or others supported by image formats