%  Giacomo Boracchi
%  23/11/2010
%

%% notation

% y  original image
% h  LTI discrete filter
% n  noise
% z  observation


%% Linear Filters

FILTER_SIZE = 11;

% load image  check the filename!!
y=im2double(imread('image_lena512.png'));
figure(1), imshow(y,[]),title('y, original image')

% create averaging filter
h=ones(FILTER_SIZE)/(FILTER_SIZE^2);

% generate "averaged" image
z = conv2(y , h , 'same');

% display
figure(1), imshow(y,[]),title('original')
figure(2), imshow(z,[]),title('observation')
figure(3), bar3(h),title('kernel'), axis tight
figure(4), 
subplot(2,1,1), hist(y(:) , 1000), title('intensities on y'), axis([0 1 0 3000])
subplot(2,1,2), hist(z(:) , 1000), title('intensities on z'), axis([0 1 0 3000])

%% shift filter

h=zeros(25);
h(13,1)=1;

z=conv2(y , h , 'same');

% display
figure(1), imshow(y),title('original')
figure(2), imshow(z),title('observation')
figure(3), bar3(h),title('kernel')

%% Gaussian Filters

h_size = 51;
sigma_gauss = 5;
h = FSPECIAL('gaussian' , h_size , sigma_gauss);
figure(3),bar3(h),title('kernel')

z=conv2(y , h , 'same');

% display
figure(1), imshow(y),title('original')
figure(2), imshow(z),title('observation')
figure(3), bar3(h),title('kernel')


%% impulse response

y = zeros(51);
y(26 , 26) = 1;

z=conv2(y , h , 'same');

% display
figure(1), imshow(y),title('original')
figure(2), bar3(z),title('observation'), axis tight
figure(3), bar3(h),title('kernel'), axis tight



%% Generating Gaussian White Noiseł
y=im2double(imread('image_lena512.png'));
figure(1), imshow(y,[]),title('y, original image')

sigma_noise = 0.1

n = sigma_noise * randn(size(y));
figure(4), imshow(n,[])

% Additive White Gaussian Noise  (AWGN)
z = y + n;
figure(5),imshow(z);

%% denoising via smoothing

% with averaging filters
figure(6),imshow(conv2(z,ones(7)/(7^2)),[])

% with gaussian smoothing 
figure(7),imshow(conv2(z,h),[])




%% Generate salt and pepper noise 
perc=0.05;
coords=rand(round(perc*size(y,1)*size(y,2)),2);
coords=ceil([size(y,1)*coords(:,1),size(y,2)*coords(:,2)])

z=y;

for ii=1:floor(size(coords,1))
    z(coords(ii,1),coords(ii,2))= round(rand(1));% 2*double(ii>size(coords,1)/2)-1;
end

% matlab generated noise    
z = imnoise(y,'salt & pepper',0.05);

% display

%% Nonlinear filters

SIZE =3;
y_hat = conv2(z , ones(SIZE) / SIZE^2 , 'valid');
y_med = medfilt2(z , [SIZE , SIZE]);
figure(1), imshow(y),title('original')
figure(2), imshow(z),title('observation')
figure(3), imshow(y_hat,[]),title('denoised via smoothing')
figure(5), imshow(y_med,[]),title('denoised via median filter')


figure(4), 
subplot(3,1,1), hist(y(:) , 100), title('intensities on y'), axis([0 1 0 3000])
subplot(3,1,2), hist(z(:) , 100), title('intensities on z'), axis([0 1 0 3000])
subplot(3,1,3), hist(y_hat(:) , 100), title('intensities on z'), axis([0 1 0 3000])