Canny邊緣檢測演算法原理及C語言實現詳解
Canny運算元是John Canny在1986年提出的，那年老大爺才28歲，該文章發表在PAMI頂級期刊上的(1986. A computational approach to edge detection. IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 8, 1986, pp. 679-698)。老大爺目前在加州伯克利做machine learning,80-90年代視覺都是影象處理，現在做視覺都是機器學習的天下，大爺的主頁(http://www.cs.berkeley.edu/~jfc/)。
Canny運算元與Marr（LoG）邊緣檢測方法類似（Marr大爺號稱計算機視覺之父），也屬於是先平滑後求導數的方法。John Canny研究了最優邊緣檢測方法所需的特性，給出了評價邊緣檢測效能優劣的三個指標：
1 好的信噪比，即將非邊緣點判定為邊緣點的概率要低，將邊緣點判為非邊緣點的概率要低；
2 高的定位效能，即檢測出的邊緣點要儘可能在實際邊緣的中心；
3 對單一邊緣僅有唯一響應，即單個邊緣產生多個響應的概率要低，並且虛假響應邊緣應該得到最大抑制。
用一句話說，就是希望在提高對景物邊緣的敏感性的同時，可以抑制噪聲的方法才是好的邊緣提取方法。
Canny運算元求邊緣點具體演算法步驟如下：
1. 用高斯濾波器平滑影象．
2. 用一階偏導有限差分計算梯度幅值和方向
3. 對梯度幅值進行非極大值抑制
4. 用雙閾值演算法檢測和連線邊緣．
具體的步驟是能容易理解，現在就是用C語言怎麼實現了，在參考了網上諸多教程的基礎下，寫了個程式碼給大家參考，肯定有不少問題，希望能得到大家的指點。
首先用白話敘述下Canny運算元的原理：看作者寫的paper題目就是邊緣檢測，何為邊緣？圖象區域性區域亮度變化顯著的部分，對於灰度影象來說，也就是灰度值有一個明顯變化，既從一個灰度值在很小的緩衝區域內急劇變化到另一個灰度相差較大的灰度值。依據我僅有的數學知識，怎麼表徵這種灰度值的變化呢？導數就是表徵變化率的，但是數字影象都是離散的，也就是導數肯定會用差分來代替。也就是具體演算法中的步驟2。但是在真實的影象中，一般會有噪聲，噪聲會影響梯度（換不嚴格說法 偏導數）的計算，所以步驟1上來先濾波。理論上將影象梯度幅值的元素值越大，說明影象中該點的梯度值越大，但這不能說明該點就是邊緣。在Canny演算法中，非極大值抑制（步驟3）是進行邊緣檢測的重要步驟，通俗意義上是指尋找畫素點區域性最大值，沿著梯度方向，比較它前面和後面的梯度值進行了
基本原理簡單說完，上程式碼，程式碼按照下面6步驟敘述。
第一步：灰度化
第二步：高斯濾波
第三步：計算梯度值和方向
第四步：非極大值抑制
第五步：雙閾值的選取
第六步：邊緣檢測
1 把彩色影象變成灰度影象。該部分主要是為像我這樣的小菜鳥準備的。
該部分是按照Canny演算法通常處理的影象為灰度圖，如果獲取的彩色影象，那首先就得進行灰度化。以RGB格式的彩圖為例，通常灰度化採用的公式是：
Gray=0.299R+0.587G+0.114B;
說個我經常出問題的程式碼：OpenCvGrayImage->imageData[i*OpenCvGrayImage->widthStep+j] 這是opencv iplimage格式通過直接訪問記憶體讀取畫素值的方式，我一直搞不清楚，i*widthStep還是j*widthStep。
記住一點，是高*widthStep就行。而且是*widthStep，而不是乘以width.如果影象的寬度不是4的倍數，opencv貌似還有補齊這一說法，所以mark一下。
程式碼如下：
1 ////////第一步：灰度化
2 IplImage * ColorImage=cvLoadImage("c:\\photo.bmp",1);
3 if (ColorImage==NULL)
4 {
5 printf("image read error");
6 return 0;
7 }
8 cvNamedWindow("Sourceimg",0);
9 cvShowImage("Sourceimg",ColorImage); //
10 IplImage * OpenCvGrayImage;
11 OpenCvGrayImage=cvCreateImage(cvGetSize(ColorImage),ColorImage->depth,1);
12 float data1,data2,data3;
13 for (int i=0;i<ColorImage->height;i++)
14 {
15 for (int j=0;j<ColorImage->width;j++)
16 {
17 data1=(uchar)(ColorImage->imageData[i*ColorImage->widthStep+j*3+0]);
18 data2=(uchar)(ColorImage->imageData[i*ColorImage->widthStep+j*3+1]);
19 data3=(uchar)(ColorImage->imageData[i*ColorImage->widthStep+j*3+2]);
20 OpenCvGrayImage->imageData[i*OpenCvGrayImage->widthStep+j]=(uchar)(0.07*data1 + 0.71*data2 + 0.21*data3);
21 }
22 }
23 cvNamedWindow("GrayImage",0);
24 cvShowImage("GrayImage",OpenCvGrayImage); //顯示灰度圖
25 cvWaitKey(0);
26 cvDestroyWindow("GrayImage");
2 對影象高斯濾波，影象高斯濾波的實現可以用兩個一維高斯核分別兩次加權實現，也就是先一維X方向卷積，得到的結果再一維Y方向卷積。當然也可以直接通過一個二維高斯核一次卷積實現。也就是二維卷積模板，由於水平有限，只說二維卷積模板怎麼算。
首先，一維高斯函式：
二維高斯函式：
是不是和一維很像，其實就是兩個一維乘積就是二維，有的書和文章中將r2=x2+y2 來表示距離的平方，也就是當前要卷積的點離核心的距離的平方，如果是3*3的卷積模板，核心就是中間的那個元素，那左上角的點到核心的距離是多少呢，就是sqrt(1+1)=sqrt(2),距離的平方 r2=2。基於這個理論，那麼模板中每一個點的高斯係數可以由上面的公式計算，這樣得到的是不是最終的模板呢？答案不是，需要歸一化，也即是每一個點的係數要除以所有係數之和，這樣才是最終的二維高斯模板。
這個裡面有個小知識點，要想計算上面的係數，需要知道高斯函式的標準差σ （sigma），還需要知道選3*3還是5*5的模板，也就是模板要多大，實際應用的時候，這兩者是有關係的，根據數理統計的知識，高斯分佈的特點就是數值分佈在（μ—3σ,μ+3σ)中的概率為0.9974，也就是模板的大小其實就是6σ這麼大就OK了，但是6σ可能不是奇數，因為我們一定要保證有核心。所以模板視窗的大小一般採用1+2*ceil(3*nSigma) ceil是向上取整函式，例如ceil(0.6)=1。
計算得到模板，那就是直接卷積就OK，卷積的意思就是影象中的點附近的模板大小區域乘以高斯模板區域，得到的結果就是該點卷積後的結果。卷積的核心意義就是獲取原始影象中像模板特徵的性質。插一句話，目前深度學習中很火的一個卷積神經網路CNN，就是利用卷積的原理，通過學習出這些卷積模板來識別檢測。
程式碼如下：
////////第二步：高斯濾波
///////
double nSigma=0.2;
int nWindowSize=1+2*ceil(3*nSigma);//通過sigma得到視窗大小
int nCenter=nWindowSize/2;
int nWidth=OpenCvGrayImage->width;
int nHeight=OpenCvGrayImage->height;
IplImage * pCanny;
pCanny=cvCreateImage(cvGetSize(ColorImage),ColorImage->depth,1);
//生成二維濾波核
double *pKernel_2=new double[nWindowSize*nWindowSize];
double d_sum=0.0;
for(int i=0;i<nWindowSize;i++)
{
for (int j=0;j<nWindowSize;j++)
{
double n_Disx=i-nCenter;//水平方向距離中心畫素距離
double n_Disy=j-nCenter;
pKernel_2[j*nWindowSize+i]=exp(-0.5*(n_Disx*n_Disx+n_Disy*n_Disy)/(nSigma*nSigma))/(2.0*3.1415926)*nSigma*nSigma;
d_sum=d_sum+pKernel_2[j*nWindowSize+i];
}
}
for(int i=0;i<nWindowSize;i++)//歸一化處理
{
for (int j=0;j<nWindowSize;j++)
{
pKernel_2[j*nWindowSize+i]=pKernel_2[j*nWindowSize+i]/d_sum;
}
}
//輸出模板
for (int i=0;i<nWindowSize*nWindowSize;i++)
{
if (i%(nWindowSize)==0)
{
printf("\n");
}
printf("%.10f ",pKernel_2[i]);
}
//濾波處理
for (int s=0;s<nWidth;s++)
{
for (int t=0;t<nHeight;t++)
{
double dFilter=0;
double dSum=0;
//當前座標（s,t）
//獲取8鄰域
for (int x=-nCenter;x<=nCenter;x++)
{
for (int y=-nCenter;y<=nCenter;y++)
{
if ((x+s>=0)&&(x+s<nWidth)&&(y+t>=0)&&(y+t<nHeight))//判斷是否越界
{
double currentvalue=(double)OpenCvGrayImage->imageData[(y+t)*OpenCvGrayImage->widthStep+x+s];
dFilter+=currentvalue*pKernel_2[x+nCenter+(y+nCenter)*nCenter];
dSum+=pKernel_2[x+nCenter+(y+nCenter)*nCenter];
}
}
}
pCanny->imageData[t*pCanny->widthStep+s]=(uchar)(dFilter/dSum);
}
}
cvNamedWindow("GaussImage",0);
cvShowImage("GaussImage",pCanny); //顯示高斯圖
cvWaitKey(0);
cvDestroyWindow("GaussImage");
3 計算梯度值和方向，影象灰度值的梯度一般使用一階有限差分來進行近似，這樣就可以得影象在x和y方向上偏導數的兩個矩陣。本文采用非常簡單的運算元，當然你也可以選擇高大上的運算元或者自己寫的：
其中f為影象灰度值，P代表X方向梯度幅值，Q代表Y方向 梯度幅值，M是該點幅值，Θ是梯度方向，也就是角度。
程式碼如下：
1 //////////////////////////////////////////////////////////
2 ////////第三步：計算梯度值和方向
3 //////////////////同樣可以用不同的檢測器////////////////加上把影象放到0-255之間//////
4 ///// P[i,j]=(S[i+1,j]-S[i,j]+S[i+1,j+1]-S[i,j+1])/2 /////
5 ///// Q[i,j]=(S[i,j]-S[i,j+1]+S[i+1,j]-S[i+1,j+1])/2 /////
6 /////////////////////////////////////////////////////////////////
7 double *P=new double[nWidth*nHeight];
8 double *Q=new double[nWidth*nHeight];
9 int *M=new int[nWidth*nHeight];
10 //IplImage * M;//梯度結果
11 //M=cvCreateImage(cvGetSize(ColorImage),ColorImage->depth,1);
12 double *Theta=new double[nWidth*nHeight];
13 int nwidthstep=pCanny->widthStep;
14 for(int iw=0;iw<nWidth-1;iw++)
15 {
16 for (int jh=0;jh<nHeight-1;jh++)
17 {
18 P[jh*nWidth+iw]=(double)(pCanny->imageData[min(iw+1,nWidth-1)+jh*nwidthstep]-pCanny->imageData[iw+jh*nwidthstep]+
19 pCanny->imageData[min(iw+1,nWidth-1)+min(jh+1,nHeight-1)*nwidthstep]-pCanny->imageData[iw+min(jh+1,nHeight-1)*nwidthstep])/2;
20 Q[jh*nWidth+iw]=(double)(pCanny->imageData[iw+jh*nwidthstep]-pCanny->imageData[iw+min(jh+1,nHeight-1)*nwidthstep]+
21 pCanny->imageData[min(iw+1,nWidth-1)+jh*nwidthstep]-pCanny->imageData[min(iw+1,nWidth-1)+min(jh+1,nHeight-1)*nwidthstep])/2;
22 }
23 }
24 //計算幅值和方向
25 for(int iw=0;iw<nWidth-1;iw++)
26 {
27 for (int jh=0;jh<nHeight-1;jh++)
28 {
29 M[jh*nWidth+iw]=(int)sqrt(P[jh*nWidth+iw]*P[jh*nWidth+iw]+Q[jh*nWidth+iw]*Q[jh*nWidth+iw]+0.5);
30 Theta[jh*nWidth+iw]=atan2(Q[jh*nWidth+iw],P[jh*nWidth+iw])*180/3.1415;
31 if (Theta[jh*nWidth+iw]<0)
32 {
33 Theta[jh*nWidth+iw]=360+Theta[jh*nWidth+iw];
34 }
35 }
36 }
4 非極大值抑制是進行邊緣檢測的一個重要步驟，通俗意義上是指尋找畫素點區域性最大值。沿著梯度方向，比較它前面和後面的梯度值進行了。見下圖。
上圖中左右圖：g1、g2、g3、g4都代表畫素點，很明顯它們是c的八領域中的4個，左圖中c點是我們需要判斷的點，藍色的直線是它的梯度方向，也就是說c要是區域性極大值，它的梯度幅值M需要大於直線與g1g2和g2g3的交點，dtmp1和dtmp2處的梯度幅值。但是dtmp1和dtmp2不是整畫素，而是亞畫素，也就是座標是浮點的，那怎麼求它們的梯度幅值呢？線性插值，例如dtmp1在g1、g2之間，g1、g2的幅值都知道，我們只要知道dtmp1在g1、g2之間的比例，就能得到它的梯度幅值，而比例是可以靠夾角計算出來的，夾角又是梯度的方向。
寫個線性插值的公式：設g1的幅值M(g1)，g2的幅值M(g2),則dtmp1可以很得到：
M(dtmp1)=w*M(g2)+(1-w)*M(g1)
其中w=distance(dtmp1,g2)/distance(g1,g2)
distance(g1,g2) 表示兩點之間的距離。實際上w是一個比例係數，這個比例係數可以通過梯度方向(幅角的正切和餘切)得到。
右邊圖中的4個直線就是4個不同的情況，情況不同，g1、g2、g3、g4代表c的八領域中的4個座標會有所差異，但是線性插值的原理都是一致的。
程式碼如下：
1 ////////第四步：非極大值抑制
2 //注意事項 權重的選取，也就是離得近權重大
3 /////////////////////////////////////////////////////////////////
4 IplImage * N;//非極大值抑制結果
5 N=cvCreateImage(cvGetSize(ColorImage),ColorImage->depth,1);
6 IplImage * OpencvCannyimg;//非極大值抑制結果
7 OpencvCannyimg=cvCreateImage(cvGetSize(ColorImage),ColorImage->depth,1);
8 int g1=0, g2=0, g3=0, g4=0; //用於進行插值，得到亞畫素點座標值
9 double dTmp1=0.0, dTmp2=0.0; //儲存兩個亞畫素點插值得到的灰度資料
10 double dWeight=0.0; //插值的權重
11
12 for(int i=1;i<nWidth-1;i++)
13 {
14 for(int j=1;j<nHeight-1;j++)
15 {
16 //如果當前點梯度為0，該點就不是邊緣點
17 if (M[i+j*nWidth]==0)
18 {
19 N->imageData[i+j*nwidthstep]=0;
20 }else
21 {
22 ////////首先判斷屬於那種情況，然後根據情況插值///////
23 ////////////////////第一種情況///////////////////////
24 ///////// g1 g2 /////////////
25 ///////// C /////////////
26 ///////// g3 g4 /////////////
27 /////////////////////////////////////////////////////
28 if((Theta[i+j*nWidth]>=90&&Theta[i+j*nWidth]<135)||(Theta[i+j*nWidth]>=270&&Theta[i+j*nWidth]<315))
29 {
30 g1=M[i-1+(j-1)*nWidth];
31 g2=M[i+(j-1)*nWidth];
32 g3=M[i+(j+1)*nWidth];
33 g4=M[i+1+(j+1)*nWidth];
34 dWeight=fabs(P[i+j*nWidth])/fabs(Q[i+j*nWidth]);
35 dTmp1=g1*dWeight+(1-dWeight)*g2;
36 dTmp2=g4*dWeight+(1-dWeight)*g3;
37 ////////////////////第二種情況///////////////////////
38 ///////// g1 /////////////
39 ///////// g2 C g3 /////////////
40 ///////// g4 /////////////
41 /////////////////////////////////////////////////////
42 }else if((Theta[i+j*nWidth]>=135&&Theta[i+j*nWidth]<180)||(Theta[i+j*nWidth]>=315&&Theta[i+j*nWidth]<360))
43 {
44 g1=M[i-1+(j-1)*nWidth];
45 g2=M[i-1+(j)*nWidth];
46 g3=M[i+1+(j)*nWidth];
47 g4=M[i+1+(j+1)*nWidth];
48 dWeight=fabs(Q[i+j*nWidth])/fabs(P[i+j*nWidth]);
49 dTmp1=g1*dWeight+(1-dWeight)*g2;
50 dTmp2=g4*dWeight+(1-dWeight)*g3;
51 ////////////////////第三種情況///////////////////////
52 ///////// g1 g2 /////////////
53 ///////// C /////////////
54 ///////// g4 g3 /////////////
55 /////////////////////////////////////////////////////
56 }else if((Theta[i+j*nWidth]>=45&&Theta[i+j*nWidth]<90)||(Theta[i+j*nWidth]>=225&&Theta[i+j*nWidth]<270))
57 {
58 g1=M[i+1+(j-1)*nWidth];
59 g2=M[i+(j-1)*nWidth];
60 g3=M[i+(j+1)*nWidth];
61 g4=M[i-1+(j+1)*nWidth];
62 dWeight=fabs(P[i+j*nWidth])/fabs(Q[i+j*nWidth]);
63 dTmp1=g1*dWeight+(1-dWeight)*g2;
64 dTmp2=g4*dWeight+(1-dWeight)*g3;
65 ////////////////////第四種情況///////////////////////
66 ///////// g1 /////////////
67 ///////// g4 C g2 /////////////
68 ///////// g3 /////////////
69 /////////////////////////////////////////////////////
70 }else if((Theta[i+j*nWidth]>=0&&Theta[i+j*nWidth]<45)||(Theta[i+j*nWidth]>=180&&Theta[i+j*nWidth]<225))
71 {
72 g1=M[i+1+(j-1)*nWidth];
73 g2=M[i+1+(j)*nWidth];
74 g3=M[i-1+(j)*nWidth];
75 g4=M[i-1+(j+1)*nWidth];
76 dWeight=fabs(Q[i+j*nWidth])/fabs(P[i+j*nWidth]);
77 dTmp1=g1*dWeight+(1-dWeight)*g2;
78 dTmp2=g4*dWeight+(1-dWeight)*g3;
79
80 }
81
82 }
83
84 if ((M[i+j*nWidth]>=dTmp1)&&(M[i+j*nWidth]>=dTmp2))
85 {
86 N->imageData[i+j*nwidthstep]=200;
87
88 }else N->imageData[i+j*nwidthstep]=0;
89
90 }
91 }
92
93
94 //cvNamedWindow("Limteimg",0);
95 //cvShowImage("Limteimg",N); //顯示非抑制
96 //cvWaitKey(0);
97 //cvDestroyWindow("Limteimg");
5 雙閾值的選取。
雙閾值的選取是按照直方圖來選擇的，首先把梯度幅值的直方圖（扯點題外話：梯度的幅值直方圖和角度直方圖也是SIFT運算元中的一個環節）求出來，選取佔直方圖總數%多少（自己定，程式碼中定義70%）所對應的梯度幅值為高閾值，高閾值的一半為低閾值，這只是一種簡單策略。也可以採用其他的。
程式碼如下：
1 ///////第五步：雙閾值的選取
2 //注意事項 注意是梯度幅值的閾值
3 /////////////////////////////////////////////////////////////////
4 int nHist[1024];//直方圖
5 int nEdgeNum;//所有邊緣點的數目
6 int nMaxMag=0;//最大梯度的幅值
7 for(int k=0;k<1024;k++)
8 {
9 nHist[k]=0;
10 }
11 for (int wx=0;wx<nWidth;wx++)
12 {
13 for (int hy=0;hy<nHeight;hy++)
14 {
15 if((uchar)N->imageData[wx+hy*N->widthStep]==200)
16 {
17 int Mindex=M[wx+hy*nWidth];
18 nHist[M[wx+hy*nWidth]]++;//獲取了梯度直方圖
19
20 }
21 }
22 }
23 nEdgeNum=0;
24 for (int index=1;index<1024;index++)
25 {
26 if (nHist[index]!=0)
27 {
28 nMaxMag=index;
29 }
30 nEdgeNum+=nHist[index];//經過non-maximum suppression後有多少邊緣點畫素
31 }
32 //計算兩個閾值 注意是梯度的閾值
33 int nThrHigh;
34 int nThrLow;
35 double dRateHigh=0.7;
36 double dRateLow=0.5;
37 int nHightcount=(int)(dRateHigh*nEdgeNum+0.5);
38 int count=1;
39 nEdgeNum=nHist[1];
40 while((nEdgeNum<=nHightcount)&&(count<nMaxMag-1))
41 {
42 count++;
43 nEdgeNum+=nHist[count];
44 }
45 nThrHigh=count;
46 nThrLow= (int)(nThrHigh*dRateLow+0.5);
47 printf("\n直方圖的長度 %d \n ",nMaxMag);
48 printf("\n梯度的閾值幅值大 %d 小閾值 %d \n ",nThrHigh,nThrLow);
6 邊緣檢測,也即是根據步驟5得到的雙閾值進行連線。
首先判斷該點是否超過高閾值，然後判斷該點的8鄰域點中尋找滿足超過低閾值的點，再根據此點收集新的邊緣，直到整個影象邊緣閉合。整個影象找完後，將非邊緣點剔除，即灰度值置0.
程式碼如下：
1 //////////////////////////////////////////////////////////
2 ////////第六步：邊緣檢測
3 /////////////////////////////////////////////////////////////////
4
5 for(int is=1;is<nWidth-1;is++)
6 {
7 for (int jt=1;jt<nHeight-1;jt++)
8 {
9 //CvScalar s=cvGet2D(N,jt,is);
10 //int currentvalue=s.val[0];
11 int currentvalue=(uchar)(N->imageData[is+jt*N->widthStep]);
12 if ((currentvalue==200)&&(M[is+jt*nWidth]>=nThrHigh))
13 //是非最大抑制後的點且 梯度幅值要大於高閾值
14 {
15 N->imageData[is+jt*nwidthstep]=255;
16 //鄰域點判斷
17 TraceEdge(is, jt, nThrLow, N, M);
18 }
19 }
20 }
21 for (int si=1;si<nWidth-1;si++)
22 {
23 for (int tj=1;tj<nHeight-1;tj++)
24 {
25 if ((uchar)N->imageData[si+tj*nwidthstep]!=255)
26 {
27 N->imageData[si+tj*nwidthstep]=0;
28 }
29 }
30
31 }
32
33 cvNamedWindow("Cannyimg",0);
34 cvShowImage("Cannyimg",N);
其中，鄰域跟蹤演算法給出了兩個，一種是遞迴，一種是非遞迴。
遞迴演算法。解決了堆疊溢位問題。之前找了很多canny程式碼參考，其中有一個版本流傳很廣，但是對於大影象堆疊溢位。
1 int TraceEdge(int w, int h, int nThrLow, IplImage* pResult, int *pMag)
2 {
3 int n,m;
4 char flag = 0;
5 int currentvalue=(uchar)pResult->imageData[w+h*pResult->widthStep];
6 if ( currentvalue== 0)
7 {
8 pResult->imageData[w+h*pResult->widthStep]= 255;
9 flag=0;
10 for (n= -1; n<=1; n++)
11 {
12 for(m= -1; m<=1; m++)
13 {
14 if (n==0 && m==0) continue;
15 int curgrayvalue=(uchar)pResult->imageData[w+n+(h+m)*pResult->widthStep];
16 int curgrdvalue=pMag[w+n+(h+m)*pResult->width];
17 if (curgrayvalue==200&&curgrdvalue>nThrLow)
18 if (TraceEdge(w+n, h+m, nThrLow, pResult, pMag))
19 {
20 flag=1;
21 break;
22 }
23 }
24 if (flag) break;
25 }
26 return(1);
27 }
28 return(0);
29 }
非遞迴演算法。如下：
1 void TraceEdge(int w, int h, int nThrLow, IplImage* pResult, int *pMag)
2 {
3 //對8鄰域畫素進行查詢
4 int xNum[8] = {1,1,0,-1,-1,-1,0,1};
5 int yNum[8] = {0,1,1,1,0,-1,-1,-1};
6 int xx=0;
7 int yy=0;
8 bool change=true;
9 while(change)
10 {
11 change=false;
12 for(int k=0;k<8;k++)
13 {
14 xx=w+xNum[k];
15 yy=h+yNum[k];
16 // 如果該象素為可能的邊界點，又沒有處理過
17 // 並且梯度大於閾值
18 int curgrayvalue=(uchar)pResult->imageData[xx+yy*pResult->widthStep];
19 int curgrdvalue=pMag[xx+yy*pResult->width];
20 if(curgrayvalue==200&&curgrdvalue>nThrLow)
21 {
22 change=true;
23 // 把該點設定成為邊界點
24 pResult->imageData[xx+yy*pResult->widthStep]=255;
25 h=yy;
26 w=xx;
27 break;
28 }
29 }
30 }
31 }
到此，整個演算法寫完了。打擊下信心，整個演算法跑起來沒問題，但是沒有opencv 的cvCanny 一個函式效果好。分析了下原因，一個是梯度運算元選的太簡單，opencv一般選用的是3*3 sobel。二是邊緣連線性還是不夠好，出現了很多斷的，也就是鄰域跟蹤演算法不夠好。希望有高手能改進。
附上程式碼：http://www.pudn.com/downloads726/sourcecode/graph/detail2903773.html