%%%%%%%%%%%%% The main.m file  %%%%%%%%%%%%%%%
function [I_temp,NVI]=colcannydiffsigma_ozel(img,Sigmax1,Sigmax2,Sigmay1,Sigmay2)
[x,map]=imread(img,'bmp');               

% The algorithm parameters:
% 1. Parameters of edge detecting filters:
%    X-axis direction filter:
Nx1=10;Nx2=10;Theta1=pi/2;
%    Y-axis direction filter:
Ny1=10;Ny2=10;Theta2=0;
% 2. The thresholding parameter alfa:
alfa=0.1;

% Get the initial image lena.bmp 512 512
wr=x(:,:,1);
wg=x(:,:,2);
wb=x(:,:,3);

% X-axis direction edge detection
filterx=d2dgauss(Nx1,Sigmax1,Nx2,Sigmax2,Theta1);
Ixr= conv2(wr,filterx,'same');
Ixg= conv2(wg,filterx,'same');
Ixb= conv2(wb,filterx,'same');
Ix=(Ixr+Ixg+Ixb)/3;


% Y-axis direction edge detection
%subplot(3,2,3)
filtery=d2dgauss(Ny1,Sigmay1,Ny2,Sigmay2,Theta2);
Iyr=conv2(wr,filtery,'same'); 
Iyg=conv2(wg,filtery,'same'); 
Iyb=conv2(wb,filtery,'same'); 
Iy=(Iyr+Iyg+Iyb)/3;

% Norm of the gradient (Combining the X and Y directional derivatives)
%subplot(3,2,4)
NVIr=sqrt(Ixr.*Ixr+Iyr.*Iyr);
NVIg=sqrt(Ixg.*Ixg+Iyg.*Iyg);
NVIb=sqrt(Ixb.*Ixb+Iyb.*Iyb);

% Thresholding
I_maxr=max(max(NVIr));
I_maxg=max(max(NVIg));
I_maxb=max(max(NVIb));

I_minr=min(min(NVIr));
I_ming=min(min(NVIg));
I_minb=min(min(NVIb));

I_max=(I_maxr+I_maxg+I_maxb)/3;
I_min=(I_minr+I_ming+I_minb)/3;

NVI=(NVIr+NVIg+NVIb)/3;

level=alfa*(I_max-I_min)+I_min;
Ibw=max(NVI,level.*ones(size(NVI)));

% Thinning (Using interpolation to find the pixels where the norms of 
% gradient are local maximum.)
[n,m]=size(Ibw);
for i=2:n-1,
   for j=2:m-1,
      if Ibw(i,j) > level,
         X=[-1,0,+1;-1,0,+1;-1,0,+1];
         Y=[-1,-1,-1;0,0,0;+1,+1,+1];
         Z=[Ibw(i-1,j-1),Ibw(i-1,j),Ibw(i-1,j+1);
            Ibw(i,j-1),Ibw(i,j),Ibw(i,j+1);
            Ibw(i+1,j-1),Ibw(i+1,j),Ibw(i+1,j+1)];
         XI=[Ix(i,j)/NVI(i,j), -Ix(i,j)/NVI(i,j)];
         YI=[Iy(i,j)/NVI(i,j), -Iy(i,j)/NVI(i,j)];
         ZI=interp2(X,Y,Z,XI,YI);
         if Ibw(i,j) >= ZI(1) & Ibw(i,j) >= ZI(2)
            I_temp(i,j)=I_max;
         else
            I_temp(i,j)=I_min;
         end
      else
         I_temp(i,j)=I_min;
      end
   end
end
I_temp=uint8(I_temp);
[szx,szy]=size(I_temp);
for i=1:szx
   for j=1:szy
      if I_temp(i,j)>0
         I_temp(i,j)=1;
      end
   end
end



%%%%%%%%%%%%%% End of the main.m file %%%%%%%%%%%%%%%