susan.txt

susan算法

------susan.m 文件------------

% -----------------------------------------------------------------------
% 
% 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<g
            R(i,j) = g-nmax;
        else
            R(i,j) = 0;
        end
    end
end

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



----------End of susan.m 文件----------

susan_threshold.m 文件

% -------------------------------------------------------------
%
% >> 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);


==========end of susan_threshold.m ==========


========================================susan改进算法===========

function [EDG] = edgemap(IM, TR, KR, NR, OP)
    % =====================================================================
    % File name   : edgemap.m
    % File Type   : m-file (script file for Matlab)
    % Requirements: Matlab Image Processing Toolbox
    % Begin       : 2006-07-07
    % Last Update : 2007-05-31
    % Author      : Nicola Asuni
    % Description : This function detects image edges using an improved 
    %               SUSAN technique.
    % Copyright   : Tecnick.com S.r.l.
    %               Via Ugo Foscolo n.19
    %               09045 Quartu Sant'Elena (CA)
    %               ITALY
    %               http://www.tecnick.com
    %               [e-mail]           
    % License     : GNU GENERAL PUBLIC LICENSE v.2
    %               http://www.gnu.org/copyleft/gpl.html
    % Version     : 1.1.000
    % =====================================================================
    
    % DESCRIPTION
    % --------------------
    % This function detects edges, which are those places in an image that
    % correspond to object boundaries. To find edges, this function looks
    % for places in the image where the intensity changes rapidly, using
    % an improved SUSAN technique.
    %
    % This algorithm is based on the technique described on:
    % S.M. Smith, J.M. Brady, "SUSAN - a new approach to low level image
    % processing", Int Journal of Computer Vision, 23(1):45-78, May 1997.
    % This technicque is subject to a patent:
    % S.M. Smith, "Method for digitally processing images to determine the
    % position of edges and/or corners therein for guidance of unmanned
    % vehicle. UK Patent 2272285. Proprietor: Secretary of State for
    % Defence, UK. 15 January 1997.
    
    % KEYWORDS
    % --------------------
    % SUSAN, image, edge, detection, detector, orientation, direction, matlab, octave.
    
    % WARNING
    % --------------------
    % This function is slow because of high computational complexity.
    
    % USAGE
    % --------------------
    % [EDG] = edgemap(IM)
    % [EDG] = edgemap(IM, TR)
    % [EDG] = edgemap(IM, TR, KR)
    % [EDG] = edgemap(IM, TR, KR, NR)
    % [EDG] = edgemap(IM, TR, KR, NR, OP)
    
    % INPUT
    % --------------------
    % IM : source image (RGB or grayscale)
    % TR : Brightness Threshold (default = 20)
    % KR : USAN Kernel Radius (nucleus excluded) (default = 3)
    % NR : EDG matrix will be normalized to this range of integers
    %      (default = 0 = not normalize)
    % OP : if true removes from USAN the pixels that are not 
    %      directly connected with the nucleus (default = false).
    %      IMPORTANT: This optimization is very slow, so use it carefully
    %      only for small images and when it's really needed. 
    
    % OUTPUT
    % --------------------
    % EDG : edge strength image
    
    % Examples
    % --------------------
    % Please check the edgexample.m and edgexample2.m files on how to use
    % this function.
    
    % NOTES
    % --------------------
    % This implementation is not intended to be used in a production
    % environment. The main purpose of this script is to clearly show how
    % this technique works. Better performaces could be obtained using a
    % compiled version or rewriting this technique using a low-level
    % programming language.
    
    % ---------------------------------------------------------------------
    
    % Some initial tests on the input arguments
    
    if (nargin < 1)
        disp('edgemap function by Nicola Asuni.');
        disp('This function returns an edge map.');
        disp('Usage:');
        disp('[EDG] = edgemap(IM)');
        disp('[EDG] = edgemap(IM, TR)');
        disp('[EDG] = edgemap(IM, TR, KR)');
        disp('[EDG] = edgemap(IM, TR, KR, NR)');
        disp('[EDG] = edgemap(IM, TR, KR, NR, OP)');
        disp('Where:');
        disp('IM : source image (RGB or grayscale)');
        disp('TR : Brightness Threshold (default = 20)');
        disp('KR : USAN Kernel Radius (nucleus excluded) (default = 3)');
        disp('NR : the EDG matrix will be normalized to this range of integers (default = 0 = not normalize)');
        disp('OP : if true removes from USAN the pixels that are not directly connected with the nucleus (default = false)');
        disp('EDG : edge strength image');
        EDG = [];
        return;
    end
    
    % assign default values
    if (nargin > 5)
        error('Too many arguments');
    end
    if (nargin < 5)
        OP = false;
    end
    if (nargin < 4)
        NR = 255;
    end
    if (nargin < 3)
        KR = 3;
    end
    if (nargin < 2)
        TR = 27;
    end
    
    % ---------------------------------------------------------------------
    
    % convert coloured images to grayscale (if needed)
    NCOLORS = ndims(IM);
    if (NCOLORS == 3)
        % RGB image
        IMG = double(rgb2gray(IM));
    elseif (NCOLORS == 2)
        IMG = double(IM);
    else
        error('Unrecognized image type, please use RGB or greyscale images');
    end
    
    % the image leves are traslated to reduce computational errors around
    % zero
    IMG = IMG + 255;
    
    % get image size
    [m,n] = size(IMG);
    
    % kernel width
    KW = (2 * KR) + 1; 
    
    % create a circular kernel mask (KM)
    KM = ones(KW,KW);
    for i = -KR:KR
        for j = -KR:KR
            if (round(sqrt((i.^2) + (j.^2))) > KR)
                KM(i+KR+1,j+KR+1) = 0;
            end
        end
    end
    
    % number of nonzero kernel elements (max kernel area)
    KAREA = nnz(KM);
    
    % calculates geometric threshold
    GT = 3 * KAREA / 4;
    
    % padding: add borders to image to simplify calculations
    IMG = [repmat(IMG(1,:),KR,1);IMG;repmat(IMG(m,:),KR,1)];
    IMG = [repmat(IMG(:,1),1,KR),IMG,repmat(IMG(:,n),1,KR)];
    
    % initialize edge strength image
    EDG = zeros(m,n);
    
    % for each image pixel
    for i = KR+1:m+KR
        for j = KR+1:n+KR
            % USAN region:
            % calculate the number of pixels within the circular mask which
            % have a similar brightness to the nucleus
            USAN = KM .* exp(-((IMG(i-KR:i+KR, j-KR:j+KR) - IMG(i,j)) / TR).^6);
            
            if (OP)
                % remove pixels that are not directly connected with the
                % nucleus (this is an improvement of the original SUSAN
                % technique).
                % select nonzero pixels of USAN region
                USAN_BINARY = ceil(USAN);
                % check if we have minimum condition to have separate regions
                if (nnz(USAN_BINARY) < (KAREA - KR))
                    USAN = bwselect(USAN_BINARY,KR+1,KR+1,8) .* USAN;
                end
            end
            
            % USAN area
            USAN_AREA = sum(sum(USAN));
            
            % subtract the USAN size from the geometric threshold GT to produce
            % an edge strength image
            if (USAN_AREA < GT)
                % this pixel is an edge
                EDG(i-KR,j-KR) = GT - USAN_AREA;    
            end
        end
    end
    
    % normalize edge values to the specified range of integers
    % (this is not included on original SUSAN technique)
    if (NR > 0)
        NSCALE = NR / max(max(EDG));
        EDG = round(NSCALE .* EDG);
        % convert data to best integer type
        if (NR > (2^32))
            EDG = uint64(EDG);
        elseif (NR > (2^16))
            EDG = uint32(EDG);
        elseif (NR > (2^8))
            EDG = uint16(EDG);
        else
            EDG = uint8(EDG);
        end
    end

To test both edge detectors and two dimensional feature detectors a test image which includes two dimensional structures of many different types has been designed. This is shown in Figure 10.