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

% image dimensions
M=100;
N=100;

% 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=20; 

% construct an original synthetic image u(i,j)
for i=1:M,
  for j=1:N,
     u(i,j)=125;
     if 40<i,
      if i<90,
      if 40<j,
      if j<90,
      u(i,j)=50;
      end
      end
      end
     end

     d=(i-40)*(i-40)+(j-40)*(j-40);
     if d<25*25,
       u(i,j)=200;
     end
  end
end

% 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