哈夫變換和Canny邊緣檢測算法

http://blog.chinaunix.net/uid-24517893-id-2846453.html  

摘   要 在圖象邊緣檢測中往往要求所檢測到的邊緣具有封閉特性， 本文詳細地分析了目前常用的兩種算法： 哈夫變換和Canny邊緣檢測算法，最後， 探討邊緣算子應滿足的准則。
關鍵詞 邊緣檢測；閉合性；哈夫變換；Canny算子

1引言

       圖象的邊緣是指圖象局部區域亮度變化顯著的部分， 該區域的灰度剖面一般可以看作是一個階躍， 既從一個灰度值在很小的緩沖區域內急劇變化到另一個灰度相差較 大的灰度值。圖象的邊緣部分集中了圖象的大部分信息， 圖象邊緣的確定與提取對於整個圖象場景的識別與理解是非常重要的 ，同時也是圖象分割所依賴的重要特 征，邊緣檢測主要是圖象的灰度變化的度量、檢測和定位， 自從1959提出邊緣檢測以來，經過五十多年的發展， 已有許多中不同的邊緣檢測方法。在我們常用的 幾種用於邊緣檢測的算子中Laplace算子常常會產生雙邊界； 而其他一些算子如Sobel算子又往往會形成不閉合區域。 本文主要討論了在邊緣檢測中，獲 取封閉邊界區域的算法。

2 圖象邊緣檢測的基本步驟

     （1）濾波。邊緣檢測主要基於導數計算，但受噪聲影響。 但濾波器在降低噪聲的同時也導致邊緣強度的損失。

     （2）增強。增強算法將鄰域中灰度有顯著變化的點突出顯示。 一般通過計算梯度幅值完成。

     （3）檢測。但在有些圖象中梯度幅值較大的並不是邊緣點。 最簡單的邊緣檢測是梯度幅值閾值判定。

     （4）定位。精確確定邊緣的位置。

圖1 邊緣檢測酸法的基本步驟

3 邊界閉合的算法

3.1 哈夫變換^[3]

      由於噪聲的存在，用各種算子得到的邊緣象素不連續， 但是由於邊緣象素之間有一定的連續性， 我們就可以根據邊緣象素在梯度幅度或梯度方向上的連續性把他們連 接起來。具體說來，如果象素(s,t)在象素(x，y) 的領域且它們的梯度幅度與梯度方向在給定的閾值下滿足：

T是幅度閾值；A是角度閾值；

     那麼， 如對所有的邊緣象素都進行上述的判斷和連接就可以得到一個閉合的 邊界。

哈夫變換方法是利用圖像得全局特性而對目標輪廓進行直接檢測的方 法，在已知區域形狀的條件下， 哈夫變換可以准確地捕獲到目標的邊界（連續的獲不連續的）， 並最終以連續曲線的形式輸出變換結果， 該變換可以從強噪聲環境中將已知形狀的目標准確得分割提取出來。

     哈夫變換的核心思想是: 點—線的對偶性(duality)。 通過變換將圖象從圖像控件轉換到參數空間, 在圖像空間中一條過點(x,y)的直線方程為y=px+q, 通過代數變換可 以轉換為另一種形式p=-px+y,即參數空間中過點(p,q) 的一條直線,如果在圖像空間中保持直線的斜率和截距的不變, 其在參數空間必定過點 (p,q),這也就說明, 在圖像空間中共線的點對應參數空間共點的線. 哈夫變換就是根據上述點— 線的對偶性把在圖象空間中存在的直線檢測問題轉換為參數空間中存 在的點檢測問題,後者的處理要比前者簡單易行得多,只需簡單地累 加統計即可實現對邊緣的檢測.

哈夫變換不僅能檢測直線等一階曲線的目標,對於園、 橢圓等高階的曲線都可以檢測出來。如圓的方程為：

    其參數空間是一個3D空間A(a，b，r)， 原理與檢測直線上的點相同，只是復雜性增加了。 如果圓的半徑r己知，則問題又回到了2D空間A(a，b)

     哈夫變換對已知目標的檢測過程受隨機噪聲和曲線中斷等不利因素的 影響很小，而且分割出的目標是直接放到另一個「干淨」 的緩存中的，因此可以做到零噪聲，是 相當有優勢的。 常規的哈夫變換在理論上能對所有可以寫出具體解析表達式的曲線進 行目標檢測，但是在實際處理時， 經常待檢測的目標不規則或是很難獲取甚至根 本沒有解析式，此時就要采取廣義上的哈夫變換來檢測目標，

3．2最優的階梯型邊緣檢測算法(canny邊緣檢測)

     1.Canny邊緣檢測基本原理

     (1)圖象邊緣檢測必須滿足兩個條件：一能有效地抑制噪聲； 二必須盡量精確確定邊緣的位置。

     (2)根據對信噪比與定位乘積進行測度，得到最優化逼近算子。 這就是Canny邊緣檢測算子。

     (3)類似與Marr（LoG）邊緣檢測方法， 也屬於先平滑後求導數的方法。

     2.Canny邊緣檢測算法：

     step1:用高斯濾波器平滑圖象；

     step2:用一階偏導的有限差分來計算梯度的幅值和方向；

     step3:對梯度幅值進行非極大值抑制；

     step4:用雙閾值算法檢測和連接邊緣。

     step1:高斯平滑函數

       

      step3:非極大值抑制

      僅僅得到全局的梯度並不足以確定邊緣，因此為確定邊緣， 必須保留局部梯度最大的點，而抑制非極大值。（non- maxima suppression,NMS）

解決方法：利用梯度的方向。

圖2非極大值抑制

四個扇區的標號為0到3，對應3*3鄰域的四種可能組合。

      在每一點上，鄰域的中心象素M與沿著梯度線的兩個象素相比。 如果M的梯度值不比沿梯度線的兩個相鄰象素梯度值大，則令M= 0。

即： 

       step4:閾值化

       減少假邊緣段數量的典型方法是對N[i，j]使用一個閾值。 將低於閾值的所有值賦零值。但問題是如何選取閾值？

       解決方法：雙閾值算法。 雙閾值算法對非極大值抑制圖象作用兩個閾值τ1和τ2， 且2τ1≒τ2，從而可以得到兩個閾值邊緣圖象N1［i,j］和 N2［i，j］。由於N2［i，j］使用高閾值得到， 因而含有很少的假邊緣，但有間斷(不閉合)。雙閾值法要在N2［ i，j］中把邊緣連接成輪廓，當到達 輪廓的端點時，該算法就在N1［i,j］ 的8鄰點位置尋找可以連接到輪廓上的邊緣，這樣， 算法不斷地在N1［i,j］中收集邊緣，直到將N2［i,j］連 接起來為止。

4 邊緣算子應滿足的准則

    若滿足此准則，就能保證單邊緣只有一個響應。

對一個算法的性能評價可分為兩個階段進行： 計算假邊緣與丟失邊緣的數目；測量用於估計位置和方向的誤差( 或誤差分布)。邊緣檢測算法的優劣也可 用品質因數( Figure of Merit)來描述。Pratt品質因數是其中一種， 它著重考慮了丟失了有效的邊緣、 邊緣定位誤差和將噪聲判斷為邊緣等三種誤差。

5 結束語

       邊緣檢測在圖象分割、模式識別、機器視覺等中都有重要作用， 人們已研究出很多種邊緣檢測算法， 而哈夫變換和canny邊緣算子等是最經典的算法， 人們已在這些經典算法基礎上提出一些新的改進算法。

參考文獻

[1]賈雲德．機器視覺［Ｍ］．北京：科學出版社，2000
[2]章毓晉．圖象處理和分析［Ｍ］．北京：清華大學出版社， 1999

[3]郎銳．數字圖象處理學．北京：希望電子出版社，2002
[4]王娜，李霞．一種新的改進Canny邊緣檢測算法． 深圳大學學報，2005，4（2），149-152

邊緣提取以及邊緣增強是不少圖像處理軟件都具有的基本功能，它的增強效果很明顯，在用於
識別的應用中，圖像邊緣也是非常重要的特征之一。圖像邊緣保留了原始圖像中相當重要的部分信息，
而又使得總的數據量減小了很多，這正符合特征提取的要求。在以後要談到的霍夫變換(檢測圖像中的幾
何形狀)中，邊緣提取就是前提步驟。

這裡我們只考慮灰度圖像，用於圖像識別的邊緣提取比起僅僅用於視覺效果增強的邊緣提取要
復雜一些。要給圖像的邊緣下一個定義還挺困難的，從人的直觀感受來說，邊緣對應於物體的邊界。圖
像上灰度變化劇烈的區域比較符合這個要求，我們一般會以這個特征來提取圖像的邊緣。但在遇到包含
紋理的圖像上，這有點問題，比如說，圖像中的人穿了黑白格子的衣服，我們往往不希望提取出來的邊
緣包括衣服上的方格。但這個比較困難，涉及到紋理圖像的處理等方法。

好了，既然邊緣提取是要保留圖像的灰度變化劇烈的區域，從數學上，最直觀的方法就是微分
(對於數字圖像來說就是差分)，在信號處理的角度來看，也可以說是用高通濾波器，即保留高頻信號。
這是最關鍵的一步，在此之前有時需要對輸入圖像進行消除噪聲的處理。

用於圖像識別的邊緣提取往往需要輸出的邊緣是二值圖像，即只有黑白兩個灰度的圖像，其中
一個灰度代表邊緣，另一個代表背景。此外，還需要把邊緣細化成只有一個像素的寬度。總的說來邊緣
提取的步驟如下：

1，去噪聲
2，微分運算
3，2值化處理
4，細化

第二步是關鍵，有不少書把第二步就直接稱為邊緣提取。實現它的算法也有很多，一般的圖像
處理教科書上都會介紹好幾種，如拉普拉茲算子，索貝爾算子，羅伯特算子等等。這些都是模板運算，
首先定義一個模板，模板的大小以3*3的較常見，也有2*2，5*5或更大尺寸的。運算時，把模板中心對
應到圖像的每一個像素位置，然後按照模板對應的公式對中心像素和它周圍的像素進行數學運算，算出
的結果作為輸出圖像對應像素點的值。

需要說明的是，模板運算是圖像的一種處理手段--鄰域處理，有許多圖像增強效果都可以采用
模板運算實現，如平滑效果，中值濾波(一種消除噪聲的方法)，油畫效果，圖像的凹凸效果等等。這些
算法都比較簡單，為人們常用。

關於前面提到的幾種邊緣提取算子(拉普拉茲算子，索貝爾算子，羅伯特算子)，教科書上都有
較為詳細的介紹，我這裡不多說了，(手頭上沒有教科書，也懶得翻譯英文資料)，如果你們有時間，可
以把這些方法的具體情況仔細介紹一下。這裡對拉普拉茲算子和索貝爾算子補充兩句。拉普拉茲算子是
2階微分算子，也就是說，相當於求取2次微分，它的精度還算比較高，但對噪聲過於敏感(有噪聲的情
況下效果很差)是它的重大缺點，所以這種算子並不是特別常用。索貝爾算子是最常用的算子之一(它是
一種一階算子)，方法簡單效果也不錯，但提取出的邊緣比較粗，要進行細化處理。另外，索貝爾算子
也可提取出圖像邊緣的方向信息來，有文章論證過，在不考慮噪聲的情況下，它取得的邊緣信息誤差不
超過7度。

順便說一句，往往我們在進行邊緣提取時只注意到位置信息，而忽略了邊緣的方向。事實上，
圖像的邊緣總有一定的走向，我們可以用邊緣曲線的法線方向(和切線垂直的直線)來代表邊緣點的方向
。在圖像識別的應用中，這個方向是非常重要的信息。


上面的幾種算子是屬於比較簡單的方法，邊緣提取的精度都不算特別高，下面介紹幾種高級算
法。首先是馬爾(Marr)算子，馬爾是計算機視覺這門學問的奠基人，很了不起，但這些理論很難懂。他
提出的邊緣提取方法可以看成兩個步驟，一個是平滑作用來消除噪聲，另一個是微分提取邊緣，也可以
說是由兩個濾波器組成，低通濾波去除噪聲，高通濾波提取邊緣。人們也稱這種方法為LOG濾波器，這也
是根據它數學表達式和濾波器形狀起的名字。也可以采用模板運算來實現這種算法，但模板的大小一般
要在7*7以上，所以運算復雜程度比索貝爾算子等要大不少，運算時間當然也長許多。

另外一種非常重要的算法是坎尼(Canny)算子,這是坎尼在1986年寫的一篇論文裡仔細論述的。
他給出了判斷邊緣提取方法性能的指標。而坎尼算子也是圖像處理領域裡的標准方法，也可以說是默認
的方法。比較奇怪的是，國內的圖像處理教科書中，介紹坎尼算子的很少。本人見過的書中，鄭南寧的
『計算機視覺與模式識別』(1998年)，算是介紹的比較詳細的。坎尼算子在使用時要提供給一些參數，
用於控制算法的性能，實際上，對於不同的圖像或不同的邊緣提取目的，應該提供不同的參數，以達到
最佳效果。它也有模板運算方法，模板的大小也比較大,和提供的參數有關，標准的大小差不多是17*17
，可以根據算子的可分離性用快速算法(否則就會慢的一塌糊涂)，坎尼算子的2值化也很有特色，具有
一定的智能性。

還有一種算法：Shen-Castan算子，大概可稱為沈峻算子，總之是中國人的成果，效果和坎尼
算子不相上下，這種算法在對邊緣提取好壞的判別標准上有些不同。(這種方法我沒用過，好象編起程
序來，要比坎尼算子還復雜)

在實際的圖像處理與識別應用中，有時需要根據被處理圖像的種類以及實際目的，量身定做算
法，邊緣提取也是一樣，但是基本原理都是一樣的。

canny算子代碼

void CreatGauss(double sigma, double **pdKernel, int *pnWidowSize);

void GaussianSmooth(SIZE sz, LPBYTE pGray, LPBYTE pResult, double sigma);

void Grad(SIZE sz, LPBYTE pGray, int *pGradX, int *pGradY, int *pMag);

void NonmaxSuppress(int *pMag, int *pGradX, int *pGradY, SIZE sz, LPBYTE pNSRst);

void EstimateThreshold(int *pMag, SIZE sz, int *pThrHigh, int *pThrLow, LPBYTE pGray, 
          double dRatHigh, double dRatLow);

void Hysteresis(int *pMag, SIZE sz, double dRatLow, double dRatHigh, LPBYTE pResult);

void TraceEdge(int y, int x, int nThrLow, LPBYTE pResult, int *pMag, SIZE sz);

void Canny(LPBYTE pGray, SIZE sz, double sigma, double dRatLow,
       double dRatHigh, LPBYTE pResult);

#include "afx.h"
#include "math.h"
#include "canny.h"

// 一維高斯分布函數，用於平滑函數中生成的高斯濾波系數
void CreatGauss(double sigma, double **pdKernel, int *pnWidowSize)
{

LONG i;

//數組中心點
int nCenter;

//數組中一點到中心點距離
double dDis;

//中間變量
double dValue;
double dSum;
dSum = 0;

// [-3*sigma,3*sigma] 以內數據，會覆蓋絕大部分濾波系數
*pnWidowSize = 1+ 2*ceil(3*sigma);

nCenter = (*pnWidowSize)/2;

*pdKernel = new double[*pnWidowSize];

//生成高斯數據
for(i=0;i<(*pnWidowSize);i++)
{
   dDis = double(i - nCenter);
   dValue = exp(-(1/2)*dDis*dDis/(sigma*sigma))/(sqrt(2*3.1415926)*sigma);
   (*pdKernel)[i] = dValue;
   dSum+=dValue;

}
//歸一化
for(i=0;i<(*pnWidowSize);i++)
{
   (*pdKernel)[i]/=dSum;
}

}

//用高斯濾波器平滑原圖像
void GaussianSmooth(SIZE sz, LPBYTE pGray, LPBYTE pResult, double sigma)
{
LONG x, y;
LONG i;

//高斯濾波器長度
int nWindowSize;

//窗口長度
int nLen;

//一維高斯濾波器
double *pdKernel;

//高斯系數與圖像數據的點乘
double dDotMul;

//濾波系數總和 
double dWeightSum;

double *pdTemp;
pdTemp = new double[sz.cx*sz.cy];

//產生一維高斯數據
CreatGauss(sigma, &pdKernel, &nWindowSize);

nLen = nWindowSize/2;

//x方向濾波
for(y=0;y<sz.cy;y++)
{
   for(x=0;x<sz.cx;x++)
   {
    dDotMul = 0;
    dWeightSum = 0;
    for(i=(-nLen);i<=nLen;i++)
    {
     //判斷是否在圖像內部
     if((i+x)>=0 && (i+x)<sz.cx)
     {
      dDotMul+=(double)pGray[y*sz.cx+(i+x)] * pdKernel[nLen+i];
      dWeightSum += pdKernel[nLen+i];
     }
    }
    pdTemp[y*sz.cx+x] = dDotMul/dWeightSum;
   }
}

//y方向濾波
for(x=0; x<sz.cx;x++)
{
   for(y=0; y<sz.cy; y++)
   {
    dDotMul = 0;
    dWeightSum = 0;
    for(i=(-nLen);i<=nLen;i++)
    {
     if((i+y)>=0 && (i+y)< sz.cy)
     {
      dDotMul += (double)pdTemp[(y+i)*sz.cx+x]*pdKernel[nLen+i];
      dWeightSum += pdKernel[nLen+i];
     }
    }
    pResult[y*sz.cx+x] = (unsigned char)dDotMul/dWeightSum;
   }
}

delete []pdKernel;
pdKernel = NULL;

delete []pdTemp;
pdTemp = NULL;

}

// 方向導數,求梯度
void Grad(SIZE sz, LPBYTE pGray,int *pGradX, int *pGradY, int *pMag)
{
LONG y,x;

//x方向的方向導數
for(y=1;y<sz.cy-1;y++)
{
   for(x=1;x<sz.cx-1;x++)
   {
    pGradX[y*sz.cx +x] = (int)( pGray[y*sz.cx+x+1]-pGray[y*sz.cx+ x-1] );
   }
}

//y方向方向導數
for(x=1;x<sz.cx-1;x++)
{
   for(y=1;y<sz.cy-1;y++)
   {
    pGradY[y*sz.cx +x] = (int)(pGray[(y+1)*sz.cx +x] - pGray[(y-1)*sz.cx +x]);
   }
}

//求梯度

//中間變量
double dSqt1;
double dSqt2;

for(y=0; y<sz.cy; y++)
{
   for(x=0; x<sz.cx; x++)
   {
    //二階范數求梯度
    dSqt1 = pGradX[y*sz.cx + x]*pGradX[y*sz.cx + x];
    dSqt2 = pGradY[y*sz.cx + x]*pGradY[y*sz.cx + x];
    pMag[y*sz.cx+x] = (int)(sqrt(dSqt1+dSqt2)+0.5);
   }
}
}

//非最大抑制
void NonmaxSuppress(int *pMag, int *pGradX, int *pGradY, SIZE sz, LPBYTE pNSRst)
{
LONG y,x;
int nPos;

//梯度分量
int gx;
int gy;

//中間變量
int g1,g2,g3,g4;
double weight;
double dTmp,dTmp1,dTmp2;

//設置圖像邊緣為不可能的分界點
for(x=0;x<sz.cx;x++)
{
   pNSRst[x] = 0;
   pNSRst[(sz.cy-1)*sz.cx+x] = 0;

}
for(y=0;y<sz.cy;y++)
{
   pNSRst[y*sz.cx] = 0;
   pNSRst[y*sz.cx + sz.cx-1] = 0;
}

for(y=1;y<sz.cy-1;y++)
{
   for(x=1;x<sz.cx-1;x++)
   {
    //當前點
    nPos = y*sz.cx + x;

    //如果當前像素梯度幅度為0，則不是邊界點
    if(pMag[nPos] == 0)
    {
     pNSRst[nPos] = 0;
    }
    else
    {
     //當前點的梯度幅度
     dTmp = pMag[nPos];

     //x,y方向導數
     gx = pGradX[nPos];
     gy = pGradY[nPos];

     //如果方向導數y分量比x分量大，說明導數方向趨向於y分量
     if(abs(gy) > abs(gx))
     {
      //計算插值比例
      weight = fabs(gx)/fabs(gy);

      g2 = pMag[nPos-sz.cx];
      g4 = pMag[nPos+sz.cx];

      //如果x,y兩個方向導數的符號相同
      //C 為當前像素，與g1-g4 的位置關系為：
      //g1 g2
      //      C
      //       g4 g3
      if(gx*gy>0)
      {
       g1 = pMag[nPos-sz.cx-1];
       g3 = pMag[nPos+sz.cx+1];
      }

      //如果x,y兩個方向的方向導數方向相反
      //C是當前像素，與g1-g4的關系為：
      //      g2 g1
      //        C
      //    g3 g4
      else
      {
       g1 = pMag[nPos-sz.cx+1];
       g3 = pMag[nPos+sz.cx-1];
      }
     }

     //如果方向導數x分量比y分量大，說明導數的方向趨向於x分量
     else
     {
      //插值比例
      weight = fabs(gy)/fabs(gx);

      g2 = pMag[nPos+1];
      g4 = pMag[nPos-1];

      //如果x,y兩個方向的方向導數符號相同
      //當前像素C與 g1-g4的關系為
      //   g3
      //   g4 C g2
      //       g1
      if(gx * gy > 0)
      {
       g1 = pMag[nPos+sz.cx+1];
       g3 = pMag[nPos-sz.cx-1];
      }
     
      //如果x,y兩個方向導數的方向相反
      // C與g1-g4的關系為
      //    g1
      //    g4 C g2
      //     g3
      else
      {
       g1 = pMag[nPos-sz.cx+1];
       g3 = pMag[nPos+sz.cx-1];
      }
     }

     //利用 g1-g4 對梯度進行插值
     {
      dTmp1 = weight*g1 + (1-weight)*g2;
      dTmp2 = weight*g3 + (1-weight)*g4;

      //當前像素的梯度是局部的最大值
      //該點可能是邊界點
      if(dTmp>=dTmp1 && dTmp>=dTmp2)
      {
       pNSRst[nPos] = 128;
      }
      else
      {
       //不可能是邊界點
       pNSRst[nPos] = 0;
      }
     }
    }
   }
}
}

// 統計pMag的直方圖，判定閾值
void EstimateThreshold(int *pMag, SIZE sz, int *pThrHigh, int *pThrLow, LPBYTE pGray, 
          double dRatHigh, double dRatLow)
{
LONG y,x,k;

//該數組的大小和梯度值的范圍有關，如果采用本程序的算法
//那麼梯度的范圍不會超過pow(2,10)
int nHist[256];

//可能邊界數
int nEdgeNum;

//最大梯度數
int nMaxMag;

int nHighCount;

nMaxMag = 0;

//初始化
for(k=0;k<256;k++)
{
   nHist[k] = 0;
}
//統計直方圖,利用直方圖計算閾值
for(y=0;y<sz.cy;y++)
{
   for(x=0;x<sz.cx;x++)
   {
    if(pGray[y*sz.cx+x]==128)
    {
     nHist[pMag[y*sz.cx+x]]++;
    }
   }
}

nEdgeNum = nHist[0];
nMaxMag = 0;

//統計經過「非最大值抑制」後有多少像素
for(k=1;k<256;k++)
{
   if(nHist[k] != 0)
   {
    nMaxMag = k;
   }

   //梯度為0的點是不可能為邊界點的
   //經過non-maximum suppression後有多少像素
   nEdgeNum += nHist[k];

}

//梯度比高閾值*pThrHigh 小的像素點總書目
nHighCount = (int)(dRatHigh * nEdgeNum + 0.5);

k=1;
nEdgeNum = nHist[1];

//計算高閾值
while((k<(nMaxMag-1)) && (nEdgeNum < nHighCount))
{
   k++;
   nEdgeNum += nHist[k];
}

*pThrHigh = k;

//低閾值
*pThrLow = (int)((*pThrHigh) * dRatLow + 0.5);

}

//利用函數尋找邊界起點
void Hysteresis(int *pMag, SIZE sz, double dRatLow, double dRatHigh, LPBYTE pResult)
{
LONG y,x;

int nThrHigh,nThrLow;

int nPos;
//估計TraceEdge 函數需要的低閾值，以及Hysteresis函數使用的高閾值
EstimateThreshold(pMag, sz,&nThrHigh,&nThrLow,pResult,dRatHigh,dRatLow);

//尋找大於dThrHigh的點，這些點用來當作邊界點，
//然後用TraceEdge函數跟蹤該點對應的邊界
for(y=0;y<sz.cy;y++)
{
   for(x=0;x<sz.cx;x++)
   {
    nPos = y*sz.cx + x;

    //如果該像素是可能的邊界點，並且梯度大於高閾值，
    //該像素作為一個邊界的起點
    if((pResult[nPos]==128) && (pMag[nPos] >= nThrHigh))
    {
     //設置該點為邊界點
     pResult[nPos] = 255;
     TraceEdge(y,x,nThrLow,pResult,pMag,sz);
    }

   }
}

//其他點已經不可能為邊界點
for(y=0;y<sz.cy;y++)
{
   for(x=0;x<sz.cx;x++)
   {
    nPos = y*sz.cx + x;

    if(pResult[nPos] != 255)
    {
     pResult[nPos] = 0;
    }
   }
}
}

//根據Hysteresis 執行的結果，從一個像素點開始搜索，搜索以該像素點為邊界起點的一條邊界的
//一條邊界的所有邊界點，函數采用了遞歸算法
//       從（x,y)坐標出發，進行邊界點的跟蹤，跟蹤只考慮pResult中沒有處理並且可能是邊界
//   點的像素（=128），像素值為0表明該點不可能是邊界點，像素值為255表明該點已經是邊界點

void TraceEdge(int y, int x, int nThrLow, LPBYTE pResult, int *pMag, SIZE sz)
{
//對8鄰域像素進行查詢
int xNum[8] = {1,1,0,-1,-1,-1,0,1};
int yNum[8] = {0,1,1,1,0,-1,-1,-1};

LONG yy,xx,k;

for(k=0;k<8;k++)
{
   yy = y+yNum[k];
   xx = x+xNum[k];

   if(pResult[yy*sz.cx+xx]==128 && pMag[yy*sz.cx+xx]>=nThrLow )
   {
    //該點設為邊界點
    pResult[yy*sz.cx+xx] = 255;

    //以該點為中心再進行跟蹤
    TraceEdge(yy,xx,nThrLow,pResult,pMag,sz);
   }
}
}

// Canny算子
void Canny(LPBYTE pGray, SIZE sz, double sigma, double dRatLow,
       double dRatHigh, LPBYTE pResult)
{
//經過高斯濾波後的圖像
LPBYTE pGaussSmooth;

pGaussSmooth = new unsigned char[sz.cx*sz.cy];

//x方向導數的指針
int *pGradX;
pGradX = new int[sz.cx*sz.cy];

//y方向
int *pGradY;
pGradY = new int[sz.cx*sz.cy];

//梯度的幅度
int *pGradMag;
pGradMag = new int[sz.cx*sz.cy];

//對原圖高斯濾波
GaussianSmooth(sz,pGray,pGaussSmooth,sigma);

//計算方向導數和梯度的幅度
Grad(sz,pGaussSmooth,pGradX,pGradY,pGradMag);

//應用非最大抑制
NonmaxSuppress(pGradMag,pGradX,pGradY,sz,pResult);

//應用Hysteresis，找到所有邊界
Hysteresis(pGradMag,sz,dRatLow,dRatHigh,pResult);

delete[] pGradX;
pGradX = NULL;
delete[] pGradY;
pGradY = NULL;
delete[] pGradMag;
pGradMag = NULL;
delete[] pGaussSmooth;
pGaussSmooth = NULL;

}

/*
void CChildWnd::OnCanny() 
{
if (! m_fOpenFile) 
{
   return;
}
m_fDone = TRUE; 
RGBToGray(szImg, aRGB, aGray, BPP);
Canny(aGray,szImg,0.1,0.9,0.76,aBinImg);
  
ShowGrayImage("l",szImg,aBinImg);
}
//*/