Matlab Implementation for SIFT "Scale-invariant feature transform "

Scale-invariant feature transform (or SIFT) is an algorithm in computer vision to detect and describe local features in images.

function [outImage,keypoints] = my_sift(A,maxS,t)
if (~exist('t', 'var'))
        t=0;
end
 if(size(A, 3)>1)
A=rgb2gray(A);
end
m=1;
sigmas=(1:1:maxS);
%A=padarray(A,[1,1]);
A=double(A);
keypoints=zeros(size(A,1)*size(A,2),4);
logImages=zeros(size(A,1),size(A,2),size(sigmas,2));
t=t*max(max(A));
for k=1 :size(sigmas,2)
    h=sigmas(1,k)*sigmas(1,k)* fspecial('log', 2*ceil(3*sigmas(1,k))+1,sigmas(1,k));
     logImages(:,:,k)= imfilter(A,h);
end

for i=2 :size(A,1)-1,
    for j=2:size(A,2)-1,
       isKeypoint=1;
        [maxVal,maxInd]=max(abs(logImages(i,j,:)));
        if(maxInd~=1 && maxInd~=maxS)
         
           for q=(i-1):i+1
               for w=(j-1):j+1
                 
                 if(abs(logImages(i,j,maxInd))<(t*maxInd*maxInd)|| abs(logImages(i,j,maxInd))< abs(logImages(q,w,maxInd-1))||abs(logImages(i,j,maxInd))< abs(logImages(q,w,maxInd))||abs(logImages(i,j,maxInd))< abs(logImages(q,w,maxInd+1)))
                     isKeypoint=0;
                     break;
                 end
               end
                if(isKeypoint==0)
                    break;
                end
           end
        else
            isKeypoint=0;
        end
   
        if(isKeypoint==1)
            keypoints(m,1)=j;
             keypoints(m,2)=i;
              keypoints(m,3)=sqrt(2)*sigmas(1,maxInd);
              keypoints(m,4)=sigmas(1,maxInd);
             m=m+1;
                 
           
        end
   end
     
   
end
%yellow = uint8([255 255 0]);
%shapeInserter = vision.ShapeInserter('Shape','Circles','BorderColor','Custom','CustomBorderColor',yellow,'LineWidth',0.1);
A=uint8(A);
outImage=A;
keypoints= keypoints(1:min(find(keypoints==0))-1,1:4);
%RGB = repmat(A(2:size(A,1)-1,2:size(A,1)-1),[1,1,3]); % convert I to an RGB image
for b=1:size(keypoints,1)  %min(find(keypoints==min(keypoints,1)))
outImage = insertShape(outImage, 'circle',keypoints(b,1:3) );


end
%yellow = uint8([255 255 0]);
%shapeInserter = vision.ShapeInserter('Shape','Circles','BorderColor','white');
%RGB = repmat(A,[1,1,3]);
%outImage = step(shapeInserter, A(2:size(A,1)-1,2:size(A,1)-1), keypoints);

end