1 下午的报告，YY做了有关Susan算法的报告，介绍了该算法的原理。
2 从zbg168的Baidu Blog中恰好看到了Susan算法以及其改进的算法。
现引如下，后面给出测试和一些实验的实例。
关于Susan算法：
Abhishek Ivaturi
E-Mail: ilv_shek@yahoo.com
Company/University: Kalmar University
算法来源：
http://www.mathworks.com/matlabcentral/fileexchange/loadAuthor.do?objectType=author&objectId=1094521
关于Susan改进算法：
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
% info@tecnick.com
=====================Susan算法的Matlab源程序代码===========================================
| Author: |
Abhishek Ivaturi |
| |
Summary: |
SUSAN Edge detection in gray scale images. |
| |
MATLAB Release: |
R13 |
| |
Required Products: |
Image Processing Toolbox |
| |
Description: |
Edge detection in gray scale images using the SUSAN algorithm. (takes some time to compute, but i hope to fix it...code is rather crude right now)Does not yet include non maximal suppresion. |
------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
% info@tecnick.com
% 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.
Figure 10: An image designed to test corner and edge detectors.
Note the corners of widely varying angles, various junctions with more than two regions meeting at a point, and the region created by taking a brightness ramp and raising a rectangle in its centre by a uniform amount. An exact description of this test image can be found in [63].
Figure 11 shows the output of the SUSAN edge detector. Note that where a reported edge appears to jitter between one side of an image edge and the other, the effect is due simply to quantization of the marked edge position. In reality, giving the figure at a higher resolution (i.e. with sub-pixel accuracy) would show the edges to lie on the image edges. It can be seen from this diagram that the SUSAN edge detector possesses the following attributes:
- Edge connectivity at junctions is complete.
- The reported edges lie exactly on the image edges. (Because of image quantization, they lie on pixels on one side of the edge or the other.)
- The edges around and inside the brightness ramp are correctly found and no false edges are reported. Thus it is shown that the edge detector copes with finding edges of regions with low or smooth variation of brightness.
The output of a standard implementation of the Canny edge detector using a very small Gaussian smoothing value (to give the best possible output) and with the final threshold optimized to just pick up the weakest correct edges is shown in Figure 12. Note the following points of
Figure 11: Output of the SUSAN edge finder (t=10) given the test image.
Figure 12: Output of the Canny edge finder (
=0.5) given the test image.
interest:
- Very few junctions involving more than two edges are correctly connected.
- Some sharper corners have broken edges.
- Incorrect edges within the brightness ramp were nearly as strong as the weakest correct edges. They have been eliminated by selecting exactly the right final threshold for binarization.
- For normal real images a larger value of
is needed, usually at least 
. This results in even greater breaks at junctions, and incorrectly rounded corners.
- The Canny edge detection is approximately 10 times slower than SUSAN edge detection.
Figure 13 shows the output of the SUSAN edge finder when Gaussian noise is added to the test image. The noise added gives a minimum signal to noise ratio of 5. This amount of noise is greater than would be expected in a normal real image by a factor of about 3. The results are good. A lot of the noise in the reported edges is only apparent, as the noise simply forces the edge to toggle between lying on one side of the real edge or the other, with no loss of accuracy. Continuity is still preserved, in even the sharpest angles (note, of course, that pixels touching diagonally constitute continuity). The output of the Canny edge finder degrades slightly more than the output of the SUSAN algorithm; see Figure 14.
Figure 13: Output of the SUSAN edge finder (t=10) given a test image with Gaussian noise of minimum SNR=5 added.
Figure 14: Output of the Canny edge finder (
=0.5) given the test image with Gaussian noise of minimum SNR=5 added.
Figures 15 and 16 show the result of running the SUSAN edge finder on a digitized real image of some printed text. Note that edges (here and in some of the other examples) have been marked in white with a black border so that all the marked edges are visible. The small (three by three) mask has been used for the SUSAN edge detector, and Canny has been carefully optimized to give the best possible results, by choosing a small value for 
and the best response threshold. Note that not only does the Canny detector have broken edges at several places, but the boundaries of the dots comprising the colon are not very clean circles. The original dots' boundaries are fairly good circles, and are better represented by the output of the SUSAN edge finder. This suggests that the SUSAN edge detector is more isotropic than a derivative based detector, as one would expect.
Figure 15: Output of the SUSAN edge finder (t=20, 3 by 3 mask) given the ``letters'' image.
Figure 16: Small sections of the outputs of the SUSAN (t=20, 3 by 3 mask) and Canny (
=0.5) edge finders when given the ``letters'' image.
Two more real examples of testing the SUSAN edge detector are given in Figures 17 and 18. Figure 17 is taken from the testing of the motion segmentation and tracking system ASSET (see [65]), where the SUSAN edge detector is used to fine-tune motion segmentation boundaries, enhancing the final accuracy. In Figure 18 note that junctions are correctly handled. In both pictures very slight breaks in the reported edges appear due to the dithering process used by the printer.
Figure 17: Output of the SUSAN edge finder (t=25) given the road image.
Figure 18: Output of the SUSAN edge finder (t=20) given a real test image. The image was created by laser printing a synthetic test image and then digitizing it using a video camera.
The results of the tests described show the SUSAN edge detector to be accurate, stable and very fast on both test and real images. It has given very good results in all the tests.
Finally, the idempotence of the SUSAN principle is investigated. It has been suggested (e.g., see [70]) that feature enhancement should be idempotent. This means that after one application of the enhancement process, further applications should make no change to the processed image. The supposed desirability of this property is based on three notions:
- Human vision interprets a scene as a series of edges, and will interpret a line drawing of the same scene in the same way. Thus it may be argued that the human visual process when fed by its interpretation of a scene will give the same interpretation. See [70].
- If a process can be applied to its output with no change then its non-linear parts (e.g. thresholding) are functioning correctly.
- If a process is idempotent then it is clear that only one application of the process is necessary to achieve the maximum ``enhancement'' possible. This could be seen as signifying efficient design.
