susancorner sh.m

MATLAB中利用SUSAN实现边缘检测

function SusanCorner(Image);   
%a = imread('d:\配准图片\ren2.bmp'); 
%figure; 
%imshow(a); 
Image=imread('F:\susancorners\zhuan.jpg');
a = double(Image); 
[height width]=size(a); 
b = zeros(7); %b为圆形模板？
susan = zeros(height,width); 
maxsum = 0.0; %响应最大极值？
for i = 1:height 
    for j = 1:width 
        center = a(i,j); %模板的中心像素值
        cgx = 0.0;  
        cgy = 0.0; 
        for k = 1:7 
            for l = 1:7 
                m = i + k - 4; %
                n = j + l - 4; 
                if(m<1) 
                    m = 1; 
                end 
                if(m>height) 
                    m = height; 
                end 
                if(n<1) 
                    n = 1; 
                end 
                if(n>width) 
                    n = width; 
                end 
                 
                if abs(a(m,n)-center)<=20 %t=20,可控制每个像素点所获得的usan区域的大小，由图像的噪声和对比度决定
                    b(k,l) = exp(-((a(m,n)-center)/20)^6); %程序是1
                else b(k,l) = 0; 
                end 
                 
                if ((m==i)&&(n==j)) 
                    b(k,l) = 0; 
                end 
                 
                cgx = cgx + m*b(k,l); %usan区域大小
                cgy = cgy + n*b(k,l); 
            end  
        end 
        susansum = sum(sum(b)); 
        if susansum>0.00000001 
            cgx = cgx/susansum; 
            cgy = cgy/susansum; 
            distance = (cgx-i)^2+(cgy-j)^2; 
        else distance = 10.0; 
        end 
         
        susan(i,j) = susansum; 
         
        if distance<2.0 
            susan(i,j) = 100; 
        elseif susansum>maxsum 
            maxsum = susansum; %寻找角点响应函数极值
        end 
    end 
end 
 
g = 3*maxsum/4; %g
r = zeros(height, width); 
for i = 1:height 
    for j = 1:width 
        if abs(susan(i,j))==0 
            max = r(i,j); 
            indexx = i; 
            indexy = j; 
            for k = -1:1 
                for l = -1:1 
                    m = i + k; 
                    n = j + l; 
                    if(m<1) 
                        m = 1; 
                    end 
                    if(m>height) 
                        m = height; 
                    end 
                    if(n<1) 
                        n = 1; 
                    end 
                    if(n>width) 
                        n = width; 
                    end 
                     
                    if r(m,n)>max 
                        max = r(m,n); 
                        indexx = m; 
                        indexy = n; 
                    end 
                end  
            end 
           % susan(indexx, indexy) = 255; 
           points=[m,n];
        end 
    end  
end 


%save ('E:\work\SUSAN\edge.jpg','susan'); 
%susan = uint8(susan); 
%figure; 
%imshow(susan);
imshow(Image)
hold on;
plot(points,'r+');