一、概要：

使用Canny邊緣檢測演算法作為例子，介紹影象的平滑方法和邊緣檢測。

Canny邊緣檢測演算法分為四步：

step1:用高斯濾波器平滑影象； step2:用一階偏導的有限差分來計算梯度的幅值和方向；(在橫豎兩個方向上計算邊緣，再求平方和的開方) step3:對梯度幅值進行非極大值抑制； step4:用雙閾值演算法檢測和連線邊緣。 

demo&效果：

原圖(lenna.bmp):

高斯濾波後的影象:

初步求邊緣後的影象：

非極大值抑制後的影象：

雙閾值檢測後的影象：

以下對這四步進行詳細介紹：

二、影象平滑與高斯濾波

影象平滑的目的是消除或儘量減少噪聲的影響，改善影象的質量。在假定加性噪聲是隨機獨立分佈的條件下，利用鄰域的平均或加權平均可以有效的抑制噪聲干擾。

方法是做鄰域計算，比如：

高斯濾波：

採用高斯函式作為加權函式。 原因一：二維高斯函式具有旋轉對稱性，保證濾波時各方向平滑程度相同； 原因二：離中心點越遠權值越小。確保邊緣細節不被模糊。

設定σ2和n，確定高斯模板權值。如σ2 =2和n=5，整數化和歸一化後得：

程式碼：

1.直接用系統函式的方法

img=imread('lenna.bmp');
imshow(img);
[m n]=size(img);
img=double(img);
%高斯濾波
w=fspecial('gaussian',[5 5]);    %[ 5 5 ]是模板尺寸，預設是[ 3 3 ]  模板即上文所提的模板
 
img=imfilter(img,w,'replicate');
figure;
imshow(uint8(img))

2.自己實現的高斯濾波（使用上面的模板）

img=imread('lenna.bmp');
img=im2double(img);
[m n]=size(img);
h=zeros(m,n);
for i=1:m
    for j=1:n
       if i<3 || i>(m-3) || j<3 || j>(n-3)
          h(i,j)=img(i,j);
       else     %程式碼長度原因，分行相加。
          h(i,j)=h(i,j)+img(i-2
 
,j-2)+img(i-2,j+2)+img(i+2,j-2)+img(i+2,j+2);
          h(i,j)=h(i,j)+(img(i-1,j-2)+img(i+1,j-2)+img(i-2,j-1)+img(i+2,j-1)+img(i-2,j+1)+img(i+2,j+1)+img(i-1,j+2)+img(i+1,j+2))*2;
          h(i,j)=h(i,j)+3*(img(i,j-2)+img(i,j+2)+img(i+2,j)+img(i-2,j));
          h(i,j)=h(i,j)+4*(img(i-1,j-1)+img(i-1,j+1)+img(i+1,j-1)+img(i+1,j+1));
          h(i,j)=h(i,j)+6*(img(i,j-1)+img(i,j+1)+img(i+1,j)+img(i-1,j));
          h(i,j)=h(i,j)+7*img(i,j);
          h(i,j)=h(i,j)/79;
       end
    end
end
imshow(h,[]);

三、初步邊緣檢測

對高斯平滑後的影象進行sobel邊緣檢測。這裡需要求橫的和豎的還有聯合的，所以一共三個需要sobel邊緣檢測影象。

程式碼：

img=imread('lenna.bmp');
[m n]=size(img);
img=double(img);
w=fspecial('gaussian',[5 5]);
img=imfilter(img,w,'replicate');
figure;
%%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);        %平方和再開方   .^表示每一位都和自己乘，不清楚的自己百度
figure;
imshow(uint8(img))

四、非極大值抑制以及雙閾值檢測邊緣

什麼是非極大值抑制？

非極大值抑制是在梯度方向上的極大值。

程式碼：（覺得看看註釋就差不多明白思路了）

img=imread('lenna.bmp');
imshow(img);
[m n]=size(img);
img=double(img);
%%canny邊緣檢測的前兩步相對不復雜，所以我就直接呼叫系統函數了
%%高斯濾波
w=fspecial('gaussian',[5 5]);
img=imfilter(img,w,'replicate');
figure;
imshow(uint8(img))
%%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);        %注意這裡不是簡單的求平均，而是平方和在開方。我曾經好長一段時間都搞錯了
figure;
imshow(uint8(img))
%%下面是非極大抑制
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
figure;
imshow(uint8(new_edge))
%%下面是滯後閾值處理
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
figure;
imshow(new_edge==255)

get_coords.m:

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

connect.m:

functionnedge=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

參考資料： 1.http://www.cnblogs.com/tiandsp/archive/2012/12/13/2817240.html 2.南京大學軟體學院數字影象處理課程PPT 一些更詳細的原理性的東西可以參考： http://blog.csdn.net/likezhaobin/article/details/6892176