canny 算子實現圖像邊緣檢測（詳細過程附源碼）

canny 算子實現圖像邊緣檢測（詳細過程附源碼）

http://www.cnblogs.com/blue-lg/archive/2011/12/25/2300899.html

canny邊緣檢測法是高斯函數的一階微分,它能在噪聲抑制和邊緣檢測之間取得較好的平衡.

環境：windows xp+matlab 2010b

時間：2011/12/25

canny算法檢測邊緣主要步驟：

1）用3x3高斯濾波器進行濾波，消除噪聲；

2）針對每一個像素，計算橫向與縱向兩方向的微分近似，以得到像素的梯度大小和方向；

3）對梯度進行"非極大抑制"（非局部最大值置0）；

4）對梯度取兩次閾值；

5）對邊緣進行連接；

下面詳細說明各個步驟：

0）讀入圖像：

clear;
clc;
i=imread('light.jpg');
k=rgb2y(i);%獲取h分量，即亮度分量

 根據邊緣的定義，邊緣檢測的目的是標識數字圖像中亮度變化明顯的點。（參考維基百科http://zh.wikipedia.org/wiki/%E8%BE%B9%E7%BC%98%E6%A3%80%E6%B5%8B）；

 

又根據公式Brightness = 0.3 * R + 0.6 * G + 0.1 * B；計算出亮度分量y;

function k=rgb2y(z)
%i必須為rgb三維矩陣
[m,n,p]=size(z);
k=zeros(m,n);
z=double(z);
for i=1:m
    for j=1:n
       k(i,j)=0.3*z(i,j,1)+0.6*z(i,j,2)+0.1*z(i,j,3);
    end
end

 

1）用3x3高斯濾波器進行濾波，消除噪聲；

function j=gaosi(i);
%i必須為二維double矩陣
j=i;
[h,w]=size(i);
for m=2;h-1
    for n=2:w-1
        j(m,n)=(i(m,n-1)+2*i(m,n)+i(m,n+1))/4;%橫向高斯濾波
    end
end

利用上面自定義的gaosi函數對圖像進行二維3x3濾波，

k1=gaosi(k);%橫向濾波
k1=k1';%對圖像進行轉置，為下一步縱向濾波作准備（縱向濾波==轉置後橫向濾波）
k1=gaosi(k1);
k1=k1';%還原

 

2）針對每一個像素，計算橫向與縱向兩方向的微分近似，以得到像素的梯度大小和方向；

利用上式，易知p,q分別為計算出的橫向、縱向的微分近似，由此再計算出梯度的大小和方向。

%計算梯度的大小和方向
[h,w]=size(k);
for m=2:h-1
    for n=2:w-1
        zz1=k1(m,n-1)+k1(m+1,n-1);
        zz2=k1(m,n)+k1(m+1,n);
        zz3=k1(m,n-1)+k1(m,n);
        zz4=k1(m+1,n-1)+k1(m+1,n);
        kp(m,n)=0.5*(zz2-zz1);
        kq(m,n)=0.5*(zz3-zz4);
        kfu(m,n)=sqrt((kp(m,n)^2)+(kq(m,n)^2));%梯度大小
        kjiao(m,n)=atan(kq(m,n)/(kp(m,n)+0.001));%梯度方向,0.001防止分母為0
    end
end

3）對梯度進行"非極大抑制"（非局部最大值置0）；

1.先將梯度方向歸類為四個主要方向，左右、上下、左斜、右斜。

%非極大值抑制
%首先將梯度方向劃分為4個方向0,45,90,135（以及他們的反向延長線）
for m=2:h-1
    for n=2:w-1
        if kjiao(m,n)>=3/8*pi
            kjiao(m,n)=2;
        else if kjiao(m,n)>=1/8*pi
            kjiao(m,n)=1;
        else if  kjiao(m,n)>=-1/8*pi
             kjiao(m,n)=0;
        else if  kjiao(m,n)>=-3/8*pi
             kjiao(m,n)=3;
            else
                kjiao(m,n)=2;
            end
            end
            end
        end
    end
end

根據劃分後的4個方向，判斷該點是否是8鄰域的局部最大值（梯度方向），比如，梯度方向為左右方向的點，判斷其是否比左右兩點的值來的大，如果不是，使該點的值為0.

%按照各個方向分別判斷
k2=k1;
for m=2:h-1
    for n=2:w-1
       if kjiao(m,n)==0
           if k1(m,n)>k1(m,n-1)&&k1(m,n)>k1(m,n+1);
           else k2(m,n)=0;
           end
       end
       if kjiao(m,n)==1
           if k1(m,n)>k1(m+1,n-1)&&k1(m,n)>k1(m-1,n+1);
           else k2(m,n)=0;
           end
       end
       if kjiao(m,n)==2
           if k1(m,n)>k1(m-1,n)&&k1(m,n)>k1(m+1,n);
           else k2(m,n)=0;
           end
       end
       if kjiao(m,n)==3
           if k1(m,n)>k1(m-1,n-1)&&k1(m,n)>k1(m+1,n+1);
           else k2(m,n)=0;
           end
       end
    end
end

4）對梯度取兩次閾值；

用兩個閾值t1和t2(t2>t1,一般取t2=2*t1),我們把梯度值小於t1的像素的灰度設為0,得到圖像1,然後我們把梯度值小於t2的像素的灰度設為0,得到圖像2。由於圖像2的閾值較高,噪音較少(但同時也損失了有用的邊緣信息,而圖像1的閾值較低,保留了較多信息,因此我們可以以圖像2為基礎,以圖像1為補充來連接圖像的邊緣。

%兩次閾值分割
k3=k2;%以t1為閾值分割後的矩陣
k4=k2;%以t2為閾值分割後的矩陣
t1=50;
t2=2*t1;
for m=2:h-1
    for n=2:w-1
       if kfu(m,n)<t1
           k3(m,n)=0;
       end
       if kfu(m,n)<t2
           k4(m,n)=0;
       end
    end
end

5）對邊緣進行連接；

a.掃描圖像2,當我們遇到一個非零值的像素p時(跟蹤以p為開始點的輪廓線直到該輪廓線的終點q;
b.在圖像1中,考察與圖像2中p點位置對應的點p'的8鄰域,如果在p'點的8鄰域中有非零像素q'存在,將其包括到圖像2中,作為點r,從r開始(重復第a步,直到我們在圖像1和圖像2中都無法繼續為止;
c.我們已經結束了對包含p的輪廓線的連接,將這條輪廓標記為已訪問過,回到第a步，尋找下一條輪廓線，重復第(a)(b)(c)步直到圖像2中再也找不到新輪廓線為止.

findline.m:

function [ff,flag1]=findline(k3,k4,flag,m,n)

flag1=flag;
  m1=m+1;n1=n+1;
  while(m~=m1||n~=n1)%若m和n都不發生變化，表明line已到終點
  flagg=0;
  for i=1:3
      if(flagg==1) break;          
      end
      for j=1:3
          if k3(m-2+i,n-2+j)~=0
              k4(m-2+i,n-2+j)=255;
              m1=m-2+i;n1=n-2+j;%新的[m,n]點
              flag1(m,n)=1;%標記已檢測過
              flagg=1;break;
      end
      end
  end
    m=m1;n=n1;
  end
  ff=k4;

主函數裡寫上：

figure;
subplot(221);imshow(i);title('原圖像');
subplot(222);imshow(k3,[]);title('閾值為50的分割圖像');
subplot(223);imshow(k4,[]);title('閾值為100的分割圖像');
flag=zeros(h,w);%標記該點是否以檢測過，1表示檢測過
for m=2:h-1
    for n=2:w-1
        if k4(m,n)~=0&&flag(m,n)==0
            [k4,flag]=findline(k3,k4,flag,m,n);
        end
    end
end   
subplot(224);imshow(k4,[]);title('修正後的分割圖像');

至此，程序完成。

效果如下：

 

至於程序中閾值的求取，大家自己定義就好了。

希望前輩們能給點學習圖像處理的經驗，謝謝了。