% dip08_2.m                 gr 2016
%    Nonlinear filters
clear all;

f = double(imread('.\img\cameraman.tif'));
[m,n] = size(f);

 fn = f + randn(m,n)*10;       % additive Gaussian noise
% fn = f .* ((rand(m,n)-.5)*.4+1);    % uniform multiplicative noise in [0.8 , 1.2]
% fn = f;  z = rand(m,n);   fn(z<.005)=0;  fn(z>1-.005)=255;    % impulse noise, prob. 2*.005
 z = rand(m,n);   fn(z<.005)=0;  fn(z>1-.005)=255;      % mixed impulse-other noise

fn(fn<0)=0; fn(fn>255)=255;  % note: actual noise variance will be less than expected

w = [1 1 1; 1 1 1; 1 1 1] / 9;

return;
%%
g1 = filter2(w,fn);              % averaging filter (arithmetic mean)
g2 = exp(filter2(w,log(fn+.01)));        % geometric mean
 imshow([f, fn; g1, g2],[0 255]);  %%%%%%%%%%%%%%%%%%%%%%%%%%%%
% note PSNR = Peak SNR = 10log(256^2/MSE) [dB]
 PSNR_noisy = 10 * log10( 256^2 / (mean(mean((fn-f).^2))) )
 PSNR_g1 = 10 * log10( 256^2 / (mean(mean((g1-f).^2))) )
 PSNR_g2 = 10 * log10( 256^2 / (mean(mean((g2-f).^2))) )
[mean(f(:))  mean(fn(:))  mean(g1(:))  mean(g2(:))]
g22 = g2 - mean(g2(:)) + mean(fn(:));  % PSNR improves setting the right luminance
 PSNR_g22 = 10 * log10( 256^2 / (mean(mean((g22-f).^2))) )

%%
g3 = medfilt2(fn, [3 3]);                 % median filter on a 3*3 support
dum(:,:,1) = fn;                % median filter on a 3*3 plus-shaped support
dum(:,:,2) = [zeros(1,n); fn(1:m-1,:)];
dum(:,:,3) = [fn(2:m,:); zeros(1,n)];
dum(:,:,4) = [zeros(m,1)  fn(:,1:n-1)];
dum(:,:,5) = [fn(:,2:n)  zeros(m,1)];
g4 = median(dum,3);
 imshow([f, fn; g3, g4],[0 255]);   %%%%%%%%%%%%%%%%%%%%%%%%%%%%
 PSNR_noisy = 10 * log10( 256^2 / (mean(mean((fn-f).^2))) )
 PSNR_g3 = 10 * log10( 256^2 / (mean(mean((g3-f).^2))) )
 PSNR_g4 = 10 * log10( 256^2 / (mean(mean((g4-f).^2))) )
 %%
dum(:,:,6) = fn;       % center-weighted plus-shape median filter, weight 3
dum(:,:,7) = fn;
g5 = median(dum,3);
 imshow([f, fn; g3, g5],[0 255]);   %%%%%%%%%%%%%%%%%%%%%%%%%%%%
 PSNR_g5 = 10 * log10( 256^2 / (mean(mean((g5-f).^2))) )
 %%
g6 = fn;
for i=2:m-1;   for j=2:n-1;            % alpha-trimmed, 3*3, d=4 or d=6
   dum = fn(i-1:i+1,j-1:j+1);
   dum = sort(dum(:));
   g6(i,j) = mean(dum(3:7));
 %  g6(i,j) = mean(dum(4:6));
   end; end;
 imshow([f, fn; g3, g6],[0 255]);   %%%%%%%%%%%%%%%%%%%%%%%%%%%%
 PSNR_g6 = 10 * log10( 256^2 / (mean(mean((g6-f).^2))) )

return;

% Note: the geom. mean filter spreads dark areas
%   --> it makes pepper (salt) noise more (less) visible

