Susan code from internet

% -----------------------------------------------------------------------
%
% This function uses the SUSAN algorithm to find edges within an image
%
%
% >>image_out = susan(image_in,threshold)
%
%
% Input parameters ... The gray scale image, and the threshold
% image_out .. (class: double) image indicating found edges
% typical threshold values may be from 10 to 30
%
%
%The following steps are performed at each image pixel:
% ( from the SUSAN webpage, http://www.fmrib.ox.ac.uk/~steve/susan/susan/node4.html )
%
% Place a circular mask around the pixel in question.
% Calculate the number of pixels within the circular mask which have similar brightness to
% the nucleus. These define the USAN.
% Subtract USAN size from geometric threshold to produce edge strength image.
%
% Estimating moments to find the edge direction has not been implemented .
% Non-maximal suppresion to remove weak edges has not been implemented yet.
%
% example:
%
% >> image_in=imread('test_pattern.tif');
% >> image = susan(image_in,27);
% >> imshow(image,[])
%
%
% Abhishek Ivaturi
%
% -------------------------------------------------------------------------

function image_out = susan(im,threshold)

%close all

%clc

% check to see if the image is a color image...
d = length(size(im));
if d==3
    image=double(rgb2gray(im));
elseif d==2
    image=double(im);
end

% mask for selecting the pixels within the circular region (37 pixels, as
% used in the SUSAN algorithm

mask = ([   0 0 1 1 1 0 0 ;
            0 1 1 1 1 1 0;
            1 1 1 1 1 1 1;
            1 1 1 1 1 1 1;
            1 1 1 1 1 1 1;
            0 1 1 1 1 1 0;
            0 0 1 1 1 0 0]); 

% the output image indicating found edges
R=zeros(size(image));

% define the USAN area
nmax = 3*37/4;

% padding the image
[a b]=size(image);
new=zeros(a+7,b+7);
[c d]=size(new);
new(4:c-4,4:d-4)=image;
for i=4:c-4
    for j=4:d-4
        current_image = new(i-3:i+3,j-3:j+3);
        current_masked_image = mask.*current_image;

%   Uncomment here to implement binary thresholding

%         current_masked_image(find(abs(current_masked_image-current_masked_image(4,4))>threshold))=0;        
%         current_masked_image(find(abs(current_masked_image-current_masked_image(4,4))<=threshold))=1;
%   This thresholding is more stable
        current_thresholded = susan_threshold(current_masked_image,threshold);
        g=sum(current_thresholded(:));
        if nmax            R(i,j) = g-nmax;
        else
            R(i,j) = 0;
        end
    end
end

image_out=R(4:c-4,4:d-4);

 

 

 

 

 

 

 

% -------------------------------------------------------------
%
%  >> thresholded_image = susan_threshold(image,threshold);
%
%    applies the thresholding scheme,
%
%              -{(I(r) - I(r0))/t}^(5/6)
%        c = e
%    to the current neighborhood of the center pixel
%    within the circle defined by the mask
%
%  Abhishek Ivaturi
% 
% ------------------------------------------------------------------

 

function thresholded = susan_threshold(image,threshold)

[a b]=size(image);
intensity_center = image((a+1)/2,(b+1)/2);

temp1 = (image-intensity_center)/threshold;
temp2 = temp1.^6;
thresholded = exp(-1*temp2);