Canny Edge Detector Algorithm

Hi I am able to NEARLY implement Canny Edge Detector but I am stuck at the connectivity part in order to make longer edges. I am using a 336*414 image and i am getting result(an image) but this is not the exact result that we get after applying the MATLAB's in built function im=edge(pic,'canny');. Please provide me with some code that connects 'Weak edges' with 'Strong Edges' as given in Gonzalez. I will be very thankful to you.

%Program that finds edges using Canny Edge Detector

clc

s=.17;

for x=-12:12 for y=-12:12 g(x+13,y+13)=exp(-(x^2+y^2)/(2*s^2)); end end

%im=imread('E:\PHD\DIP3E_Original_Images_CH10\Fig1016(a).tif');

im=imread('E:\PHD\MATLAB Images\cannygray2.jpg');

%imshow(im);

%Image smoothing using Gaussian Filter for x=1:336 for y=1:414 sum2=0; for s=12:-1:-12 xps=x-s; if(xps<=0) xps=1; end if(xps>336) xps=336; end for t=12:-1:-12 ypt=y+t; if(ypt<=0) ypt=1; end if(ypt>414) ypt=414; end sum2=sum2+g(s+13,t+13)*double(im(xps,ypt));

            end
        end

         dd(x,y)=sum2; 

    end
end

%Display Gaussian Image % figure; % imshow(uint8(dd));

%Findding the Gradient Image using Sobel Operator
gx=[-1 -2 -1;0 0 0;1 2 1];
gy=[-1 0 1;-2 0 2;-1 0 1];

for x=1:336
  for y=1:414
      sum3=0;
      sum4=0;
      for s=-1:1
          xps=x+s;
          if(xps<=0)
              xps=1;
          end
          if(xps>336)
             xps=336;
          end
          for t=-1:1
              ypt=y+t;
              if(ypt<=0)
                  ypt=1;
              end
              if(ypt>414)
                  ypt=414;
              end
               sum3=sum3+gx(s+2,t+2)*double(dd(xps,ypt));
               sum4=sum4+gy(s+2,t+2)*double(dd(xps,ypt));

            end
        end
        sum5=sqrt(sum3^2+sum4^2);
        sum6=atan(sum4/sum3)*(180/pi);
        ee(x,y)=sum5; 
        aa(x,y)=sum6;

    end
  end

%Display Gradient(Magnitude)Image % figure; % imshow(uint8(ee)); %

 %Display Gradient (Direction)Image
%   figure;
%   imshow(uint8(aa));

for x=1:336 for y=1:414 gn(x,y)=0; end end

 %Non-Maxima Supression
 for x=1:336
    for y=1:414
        for s=-1:1
            xps=x+s;
            if(xps<=0)
                xps=1;
            end
            if(xps>336)
               xps=336;
            end
            for t=-1:1
                ypt=y+t;
                if(ypt<=0)
                    ypt=1;
                end
                if(ypt>414)
                    ypt=414;
                end

               kk(s+2,t+2)=ee(xps,ypt);  

            end
        end

        %As we have already found the direction (edge normal or gradient vector direction) at each and every pixel
        %location (all directions are stored in aa(x,y) image), now we need
        %to check the neighboring pixel values in respective direction dk
        %(like if direction dk is horizontal i.e 90 degree wrt to x axis ( this is actually y axis in our daily sense),
        %then compare the values of (2,1) and (2,3) with sum5.
        %If sum5 (which is sqrt(gx^2 + gy^2)) is greater than both the neighbors in the
        %direction dk, then gn(x,y)=sum5 otherwise gn(x,y)=0.
        if((aa(x,y)>=-22.5 && aa(x,y)<=22.5) || (aa(x,y)>=-157.5 && aa(x,y)<=157.5))
            if(ee(x,y)>kk(1,2) && ee(x,y)>kk(3,2))
                gn(x,y)=ee(x,y);
            else
                gn(x,y)=0;
            end
        end

        if((aa(x,y)>=67.5 && aa(x,y)<=112.5) || (aa(x,y)>=-112.5 && aa(x,y)<=-67.5))
            if(ee(x,y)>kk(2,1) && ee(x,y)>kk(2,3))
                gn(x,y)=ee(x,y);
            else
                gn(x,y)=0;
            end
        end

        if((aa(x,y)>=22.5 && aa(x,y)<67.5) || (aa(x,y)>=-157.5 && aa(x,y)<-112.5))
            if(ee(x,y)>kk(1,1) && ee(x,y)>kk(3,3))
                gn(x,y)=ee(x,y);
            else
                gn(x,y)=0;
            end
        end

        if((aa(x,y)>112.5 && aa(x,y)<=157.5) || (aa(x,y)>-67.5 && aa(x,y)<=-22.5))
            if(ee(x,y)>kk(1,3) && ee(x,y)>kk(3,1))
                gn(x,y)=ee(x,y);
            else
                gn(x,y)=0;
            end

        end
    end
 end

% figure; % imshow(uint8(gn));

TH=150;
TL=60;

for x=1:336
  for y=1:414
      gnh(x,y)=0;
  end
end

for x=1:336
  for y=1:414
      gnl(x,y)=0;
  end
end

 for x=1:336
    for y=1:414
        if(gn(x,y)>=TH)
            gnh(x,y)=gn(x,y);
        end
    end
 end

% figure; % imshow(uint8(gnh)); %

 for x=1:336
    for y=1:414
        if(gn(x,y)>=TL)
            gnl(x,y)=gn(x,y);
        end
    end
 end

% figure; % imshow(uint8(gnl)); % gnl2=gnl-gnh; % % figure; % imshow(uint8(gnl2)); %

    for x=1:336
        for y=1:414
            gnl4(x,y)=0;
        end
    end

for x=1:336
   for y=1:414
       if(gnh(x,y)>0)
          for s=-1:1
               xps=x+s;
               if(xps==0)
                   xps=1;
               end
               if(xps>336)
                    xps=336;
               end
               for t=-1:1
                   ypt=y+t;
                   if(ypt==0)
                       ypt=1;
                   end
                   if(ypt>414)
                       ypt=414;
                   end

                     %Selecting a 3*3 window for comparison as done below
                     kk(s+2,t+2)=gnl2(xps,ypt);  

                 end
             end

         for s=-1:1
             xps2=x+s;
             if(xps2==0)
                xps2=1;
             end
             if(xps2>336)
                xps2=336;
             end
             for t=-1:1
                 ypt2=y+t;
                 if(ypt2==0)
                    ypt2=1;
                 end
                 if(ypt2>414)
                    ypt2=414;
                 end

                 if(kk(s+2,t+2)>0)
                    gnl4(xps2,ypt2)=kk(s+2,t+2);
                 else
                     gnl4(xps2,ypt2)=0;
                 end

             end
         end
         end
     end
  end

   gf=gnh+gnl4;

   for x=1:336
      for y=1:414
          if(gf(x,y)>=150)
              gf(x,y)=255;
          else
              gf(x,y)=0;
          end
      end
   end

     figure;
     imshow(uint8(gf));