Canny邊緣檢測算法

<strong><span style="font-size:18px;"><span style="color:#6600cc;">clear all;
close all;
clc;

img=rgb2gray(imread('9.jpg'));
subplot(241);imshow(img);title('原圖')
[m n]=size(img);
img=double(img);

%%canny邊緣檢測的前兩步相對不複雜，是以我就直接調用系統函數了
%%高斯濾波
w=fspecial('gaussian',[5 5]);
img=imfilter(img,w,'replicate');
subplot(242);imshow(uint8(img));title('高斯濾波')

%%sobel邊緣檢測
w=fspecial('sobel');
img_w=imfilter(img,w,'replicate');      %求橫邊緣
w=w';%轉置
img_h=imfilter(img,w,'replicate');      %求豎邊緣
img=sqrt(img_w.^2+img_h.^2);        %注意這裡不是簡單的求平均，而是平方和在開方。我曾經好長一段時間都搞錯了
subplot(243);imshow(uint8(img));title('sobel邊緣檢測')


%%下面是非極大抑制
new_edge=zeros(m,n);
for i=2:m-1
    for j=2:n-1
        Mx=img_w(i,j);
        My=img_h(i,j);
        
        if My~=0
            o=atan(Mx/My);      %邊緣的法線弧度
        elseif My==0 && Mx>0
            o=pi/2;
        else
            o=-pi/2;            
        end
        
        %Mx處用My和img進行插值
        adds=get_coords(o);      %邊緣像素法線一側求得的兩點坐标，插值需要       
        M1=My*img(i+adds(2),j+adds(1))+(Mx-My)*img(i+adds(4),j+adds(3));   %插值後得到的像素，用此像素和目前像素比較 
        adds=get_coords(o+pi);%邊緣法線另一側求得的兩點坐标，插值需要
        M2=My*img(i+adds(2),j+adds(1))+(Mx-My)*img(i+adds(4),j+adds(3));   %另一側插值得到的像素，同樣和目前像素比較
        
        isbigger=(Mx*img(i,j)>M1)*(Mx*img(i,j)>=M2)+(Mx*img(i,j)<M1)*(Mx*img(i,j)<=M2); %如果目前點比兩邊點都大置1
        
        if isbigger
           new_edge(i,j)=img(i,j); 
        end        
    end
end
subplot(244);imshow(uint8(new_edge));title('非極大抑制')


%%下面是滞後門檻值處理
up=120;     %上門檻值
low=100;    %下門檻值
set(0,'RecursionLimit',10000);  %設定最大遞歸深度
for i=1:m
    for j=1:n
      if new_edge(i,j)>up &&new_edge(i,j)~=255  %判斷上門檻值
            new_edge(i,j)=255;
            new_edge=connect(new_edge,i,j,low);
      end
    end
end

subplot(245);imshow(new_edge==255);title('滞後門檻值處理')

</span></span></strong><pre name="code" class="plain"><span style="color:#009900;">function nedge=connect(nedge,y,x,low)       %種子定位後的連通分析
    neighbour=[-1 -1;-1 0;-1 1;0 -1;0 1;1 -1;1 0;1 1];  %八連通搜尋
    [m n]=size(nedge);
    for k=1:8
        yy=y+neighbour(k,1);
        xx=x+neighbour(k,2);
        if yy>=1 &&yy<=m &&xx>=1 && xx<=n  
            if nedge(yy,xx)>=low && nedge(yy,xx)~=255   %判斷下門檻值
                nedge(yy,xx)=255;
                nedge=connect(nedge,yy,xx,low);
            end
        end        
    end 

end</span>
           

<span style="color:#3333ff;">function re=get_coords(angle)       %angle是邊緣法線角度，傳回法線前後兩點
    sigma=0.000000001;
    x1=ceil(cos(angle+pi/8)*sqrt(2)-0.5-sigma);
    y1=ceil(-sin(angle-pi/8)*sqrt(2)-0.5-sigma);
    x2=ceil(cos(angle-pi/8)*sqrt(2)-0.5-sigma);
    y2=ceil(-sin(angle-pi/8)*sqrt(2)-0.5-sigma);
    re=[x1 y1 x2 y2];

end</span>