canny 算子實現圖像邊緣檢測（詳細過程附源碼）
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('修正後的分割圖像');
至此，程序完成。
效果如下：
至於程序中閾值的求取，大家自己定義就好了。
希望前輩們能給點學習圖像處理的經驗，謝謝了。