%% Blur, 
% observation model
% 
%  z= (h \circledast y) (noise free observation) 
%
%  z= (h \circledast y) +n (noisy observation) 

%% Differentiating Filters
y = im2double(imread('image_lena512.png'));
figure(1), imshow(y),title('original image');

dx=[1 0 -1; 1 0 -1; 1 0 -1];
dy=dx';

% compute gradient components
Gx=conv2(y , dx , 'same');
Gy=conv2(y , dy , 'same');

figure(3),imshow(Gx , []),title('horizontal derivative')
figure(4),imshow(Gy , []),title('vertical derivative')

% Gradient Norm
M=sqrt(Gx.^2 + Gy.^2);

figure(5),imshow(M,[]),title('gradient magnitude')



%%% Nonlinear Filtering
%%% Median Filtering

%% Edges
disp('differentiating filters')
diffx=[1 -1]
diffy=diffx'

%smoothing filters Previtt
sx=ones(2,3);
sy=sx';

disp('derivative filters Previtt')
dx=conv2(sy , diffx)
dy=conv2(sx , diffy)
% conv2(sy,diffx,'valid')

%smoothing filters Sobel
sx=[1 2 1 ; 1 2 1];
sy=sx';

disp(' derivative filters Sobel')
dx=conv2(sy,diffx)
dy=conv2(sx,diffy)


%% Derivative
y_x=conv2(y , dx , 'same');
y_y=conv2(y , dy , 'same');

Grad=zeros(size(y,1),size(y,2),2);
Grad(:,:,1)=y_x;
Grad(:,:,2)=y_y;

% grad_r_c = Grad(r,c,:);

Norm_Grad = sqrt(Grad(: , : , 1) .^ 2 + Grad(: , : , 2) .^ 2);
Dir_Grad=atand(-sign(Grad(:,:,1)).*Grad(:,:,2)...
                              ./(abs(Grad(:,:,1))+eps));

figure(2),imshow(Norm_Grad,[]),title('norm');
figure(3),imshow(Dir_Grad,[]),title('directions');

%% Hard Threshold  - Play on the choice of the Threshold T
% T=3*median(Norm_Grad(:));
% Norm_Grad(Norm_Grad<T)=0;

% Binary Threshold  - Play on the choice of the Threshold T
% T = 5*median(Norm_Grad(:));
[T , eff ] = graythresh(Norm_Grad(:));
[hh , hh_bins] = hist(Norm_Grad(:) , 1000);
figure(7), bar(hh_bins , hh)
hold on
plot([T , T] , [min(hh) , max(hh)] , 'g' ,'LineWidth' , 3)
hold off

Norm_Grad_Binary=double(Norm_Grad> T);
Norm_Grad_Hard=Norm_Grad.*Norm_Grad_Binary;

figure(3),imshow(Norm_Grad_Binary,[]),title('Binary Thresholded gradient norm');
%%
figure(2),imshow(Norm_Grad,[]),title('Thresholded gradient norm');

figure(4),imshow(Norm_Grad_Hard,[]),title(' Hard Thresholded gradient norm');
figure(5),imshow(Norm_Grad_Binary.*Dir_Grad,[])

%% dilate 

% dilate via convolution
% Norm_Grad_mask=double(conv2(Norm_Grad_Binary,ones(3),'same')>0);
Norm_Grad_mask = imerode(Norm_Grad_Binary , ones(3));
Norm_Grad_mask = imdilate(Norm_Grad_mask , ones(5));
figure(3),imshow(Norm_Grad_mask),title('Mask on gradient norm');

% show directions
figure(4),imshow(Dir_Grad.*Norm_Grad_mask,[]),title('Gradient direction'),impixelinfo;