3、 Canny算法的實現(xiàn)流程
由于本文主要目的在于學(xué)習(xí)和實現(xiàn)算法，而對于圖像讀取、視頻獲取等內(nèi)容不進(jìn)行闡述。因此選用OpenCV算法庫作為其他功能的實現(xiàn)途徑（關(guān)于OpenCV的使用，作者將另文表述）。首先展現(xiàn)本文將要處理的彩色圖片。
圖2 待處理的圖像
3.1 圖像讀取和灰度化
編程時采用上文所描述的第二種方法來實現(xiàn)圖像的灰度化。其中ptr數(shù)組中保存的灰度化后的圖像數(shù)據(jù)。具體的灰度化后的效果如圖3所示。
- IplImage* ColorImage = cvLoadImage( "12.jpg", -1 );
- IplImage* OpenCvGrayImage;
- unsigned char* ptr;
- if (ColorImage == NULL)
- return;
- int i = ColorImage->width * ColorImage->height;
- BYTE data1;
- BYTE data2;
- BYTE data3;
- ptr = new unsigned char[i];
- for(intj=0; j<ColorImage->height; j++)
- {
- for(intx=0; x<ColorImage->width; x++)
- {
- data1 = (BYTE)ColorImage->imageData[j*ColorImage->widthStep + i*3];
- data2 = (BYTE)ColorImage->imageData[j*ColorImage->widthStep + i*3 + 1];
- data3 = (BYTE)ColorImage->imageData[j*ColorImage->widthStep + i*3 + 2];
- ptr[j*ColorImage->width+x]=(BYTE)(0.072169*data1 + 0.715160*data2 + 0.212671*data3);
- }
- }
- OpenCvGrayImage=cvCreateImageHeader(cvGetSize(ColorImage), ColorImage->depth, 1);
- cvSetData(GrayImage,ptr, GrayImage->widthStep);
- cvNamedWindow("GrayImage",CV_WINDOW_AUTOSIZE);
- cvShowImage("GrayImage",OpenCvGrayImage);
- cvWaitKey(0);
- cvDestroyWindow("GrayImage");
圖3 灰度化后的圖像
3.2 圖像的高斯濾波
根據(jù)上面所講的邊緣檢測過程，下一個步驟就是對圖像進(jìn)行高斯濾波?？筛鶕?jù)之前博文描述的方法獲取一維或者二維的高斯濾波核。因此進(jìn)行圖像高斯濾波可有兩種實現(xiàn)方式，以下具體進(jìn)行介紹。 首先定義該部分的通用變量：
- double nSigma = 0.4;
- int nWidowSize = 1+2*ceil(3*nSigma);
- int nCenter = (nWidowSize)/2;
兩種方法都需要用到的變量：
- int nWidth = OpenCvGrayImage->width;
- int nHeight = OpenCvGrayImage->height;
- unsigned char* nImageData = new unsigned char[nWidth*nHeight];
- unsigned char*pCanny = new unsigned char[nWidth*nHeight];
- double* nData = new double[nWidth*nHeight];
- for(int j=0; j<nHeight; j++)
- {
- for(i=0; i<nWidth; i++)
- nImageData[j*nWidth+i] = (unsigned char)OpenCvGrayImage->imageData[j*nWidth+i];
- }
3.2.1 根據(jù)一維高斯核進(jìn)行兩次濾波
1）生成一維高斯濾波系數(shù)
-
- double* pdKernal_1 = new double[nWidowSize];
- double dSum_1 = 0.0;
-
-
-
-
-
-
-
-
- for(int i=0; i<nWidowSize; i++)
- {
- double nDis = (double)(i-nCenter);
- pdKernal_1[i] = exp(-(0.5)*nDis*nDis/(nSigma*nSigma))/(sqrt(2*3.14159)*nSigma);
- dSum_1 += pdKernal_1[i];
- }
- for(i=0; i<nWidowSize; i++)
- {
- pdKernal_1[i] /= dSum_1;
- }
2）分別進(jìn)行x向和y向的一維加權(quán)濾波，濾波后的數(shù)據(jù)保存在矩陣pCanny中
- for(i=0; i<nHeight; i++)
- {
- for(j=0; j<nWidth; j++)
- {
- double dSum = 0;
- double dFilter=0;
- for(int nLimit=(-nCenter); nLimit<=nCenter; nLimit++)
- {
- if((j+nLimit)>=0 && (j+nLimit) < nWidth )
- {
- dFilter += (double)nImageData[i*nWidth+j+nLimit] * pdKernal_1[nCenter+nLimit];
- dSum += pdKernal_1[nCenter+nLimit];
- }
- }
- nData[i*nWidth+j] = dFilter/dSum;
- }
- }
-
- for(i=0; i<nWidth; i++)
- {
- for(j=0; j<nHeight; j++)
- {
- double dSum = 0.0;
- double dFilter=0;
- for(int nLimit=(-nCenter); nLimit<=nCenter; nLimit++)
- {
- if((j+nLimit)>=0 && (j+nLimit) < nHeight)
- {
- dFilter += (double)nData[(j+nLimit)*nWidth+i] * pdKernal_1[nCenter+nLimit];
- dSum += pdKernal_1[nCenter+nLimit];
- }
- }
- pCanny[j*nWidth+i] = (unsigned char)(int)dFilter/dSum;
- }
- }
3.2.2 根據(jù)二維高斯核進(jìn)行濾波
1）生成二維高斯濾波系數(shù)
-
- double* pdKernal_2 = new double[nWidowSize*nWidowSize];
- double dSum_2 = 0.0;
-
-
-
-
-
-
-
- for(i=0; i<nWidowSize; i++)
- {
- for(int j=0; j<nWidowSize; j++)
- {
- int nDis_x = i-nCenter;
- int nDis_y = j-nCenter;
- pdKernal_2[i+j*nWidowSize]=exp(-(1/2)*(nDis_x*nDis_x+nDis_y*nDis_y)
- /(nSigma*nSigma))/(2*3.1415926*nSigma*nSigma);
- dSum_2 += pdKernal_2[i+j*nWidowSize];
- }
- }
- for(i=0; i<nWidowSize; i++)
- {
- for(int j=0; j<nWidowSize; j++)
- {
- pdKernal_2[i+j*nWidowSize] /= dSum_2;
- }
- }
2）采用高斯核進(jìn)行高斯濾波，濾波后的數(shù)據(jù)保存在矩陣pCanny中
- int x;
- int y;
- for(i=0; i<nHeight; i++)
- {
- for(j=0; j<nWidth; j++)
- {
- double dFilter=0.0;
- double dSum = 0.0;
- for(x=(-nCenter); x<=nCenter; x++)
- {
- for(y=(-nCenter); y<=nCenter; y++)
- {
- if( (j+x)>=0 && (j+x)<nWidth && (i+y)>=0 && (i+y)<nHeight)
- {
- dFilter += (double)nImageData [(i+y)*nWidth + (j+x)]
- * pdKernal_2[(y+nCenter)*nWidowSize+(x+nCenter)];
- dSum += pdKernal_2[(y+nCenter)*nWidowSize+(x+nCenter)];
- }
- }
- }
- pCanny[i*nWidth+j] = (unsigned char)dFilter/dSum;
- }
- }
3.3 圖像增強(qiáng)——計算圖像梯度及其方向
根據(jù)上文分析可知，實現(xiàn)代碼如下
-
-
-
-
- double* P = new double[nWidth*nHeight];
- double* Q = new double[nWidth*nHeight];
- int* M = new int[nWidth*nHeight];
- double* Theta = new double[nWidth*nHeight];
-
- for(i=0; i<(nHeight-1); i++)
- {
- for(j=0; j<(nWidth-1); j++)
- {
- P[i*nWidth+j] = (double)(pCanny[i*nWidth + min(j+1, nWidth-1)] - pCanny[i*nWidth+j] + pCanny[min(i+1, nHeight-1)*nWidth+min(j+1, nWidth-1)] - pCanny[min(i+1, nHeight-1)*nWidth+j])/2;
- Q[i*nWidth+j] = (double)(pCanny[i*nWidth+j] - pCanny[min(i+1, nHeight-1)*nWidth+j] + pCanny[i*nWidth+min(j+1, nWidth-1)] - pCanny[min(i+1, nHeight-1)*nWidth+min(j+1, nWidth-1)])/2;
- }
- }
-
- for(i=0; i<nHeight; i++)
- {
- for(j=0; j<nWidth; j++)
- {
- M[i*nWidth+j] = (int)(sqrt(P[i*nWidth+j]*P[i*nWidth+j] + Q[i*nWidth+j]*Q[i*nWidth+j])+0.5);
- Theta[i*nWidth+j] = atan2(Q[i*nWidth+j], P[i*nWidth+j]) * 57.3;
- if(Theta[i*nWidth+j] < 0)
- Theta[i*nWidth+j] += 360;
- }
- }
3.4 非極大值抑制
根據(jù)上文所述的工作原理，這部分首先需要求解每個像素點在其鄰域內(nèi)的梯度方向的兩個灰度值，然后判斷是否為潛在的邊緣，如果不是則將該點灰度值設(shè)置為0.
首先定義相關(guān)的參數(shù)如下：
- unsigned char* N = new unsigned char[nWidth*nHeight];
- int g1=0, g2=0, g3=0, g4=0;
- double dTmp1=0.0, dTmp2=0.0;
- double dWeight=0.0;
其次，對邊界進(jìn)行初始化：
- for(i=0; i<nWidth; i++)
- {
- N[i] = 0;
- N[(nHeight-1)*nWidth+i] = 0;
- }
- for(j=0; j<nHeight; j++)
- {
- N[j*nWidth] = 0;
- N[j*nWidth+(nWidth-1)] = 0;
- }
進(jìn)行局部最大值尋找，根據(jù)上文圖1所述的方案進(jìn)行插值，然后判優(yōu)，實現(xiàn)代碼如下：
- for(i=1; i<(nWidth-1); i++)
- {
- for(j=1; j<(nHeight-1); j++)
- {
- int nPointIdx = i+j*nWidth;
- if(M[nPointIdx] == 0)
- N[nPointIdx] = 0;
- else
- {
-
-
-
-
-
-
- if( ((Theta[nPointIdx]>=90)&&(Theta[nPointIdx]<135)) ||
- ((Theta[nPointIdx]>=270)&&(Theta[nPointIdx]<315)))
- {
-
- g1 = M[nPointIdx-nWidth-1];
- g2 = M[nPointIdx-nWidth];
- g3 = M[nPointIdx+nWidth];
- g4 = M[nPointIdx+nWidth+1];
- dWeight = fabs(P[nPointIdx])/fabs(Q[nPointIdx]);
- dTmp1 = g1*dWeight+g2*(1-dWeight);
- dTmp2 = g4*dWeight+g3*(1-dWeight);
- }
-
-
-
-
-
- else if( ((Theta[nPointIdx]>=135)&&(Theta[nPointIdx]<180)) ||
- ((Theta[nPointIdx]>=315)&&(Theta[nPointIdx]<360)))
- {
- g1 = M[nPointIdx-nWidth-1];
- g2 = M[nPointIdx-1];
- g3 = M[nPointIdx+1];
- g4 = M[nPointIdx+nWidth+1];
- dWeight = fabs(Q[nPointIdx])/fabs(P[nPointIdx]);
- dTmp1 = g2*dWeight+g1*(1-dWeight);
- dTmp2 = g4*dWeight+g3*(1-dWeight);
- }
-
-
-
-
-
- else if( ((Theta[nPointIdx]>=45)&&(Theta[nPointIdx]<90)) ||
- ((Theta[nPointIdx]>=225)&&(Theta[nPointIdx]<270)))
- {
- g1 = M[nPointIdx-nWidth];
- g2 = M[nPointIdx-nWidth+1];
- g3 = M[nPointIdx+nWidth];
- g4 = M[nPointIdx+nWidth-1];
- dWeight = fabs(P[nPointIdx])/fabs(Q[nPointIdx]);
- dTmp1 = g2*dWeight+g1*(1-dWeight);
- dTmp2 = g3*dWeight+g4*(1-dWeight);
- }
-
-
-
-
-
- else if( ((Theta[nPointIdx]>=0)&&(Theta[nPointIdx]<45)) ||
- ((Theta[nPointIdx]>=180)&&(Theta[nPointIdx]<225)))
- {
- g1 = M[nPointIdx-nWidth+1];
- g2 = M[nPointIdx+1];
- g3 = M[nPointIdx+nWidth-1];
- g4 = M[nPointIdx-1];
- dWeight = fabs(Q[nPointIdx])/fabs(P[nPointIdx]);
- dTmp1 = g1*dWeight+g2*(1-dWeight);
- dTmp2 = g3*dWeight+g4*(1-dWeight);
- }
- }
-
- if((M[nPointIdx]>=dTmp1) && (M[nPointIdx]>=dTmp2))
- N[nPointIdx] = 128;
- else
- N[nPointIdx] = 0;
- }
- }
3.5雙閾值檢測實現(xiàn)
1）定義相應(yīng)參數(shù)如下
- int nHist[1024];
- int nEdgeNum;
- int nMaxMag = 0;
- int nHighCount;
2）構(gòu)造灰度圖的統(tǒng)計直方圖，根據(jù)上文梯度幅值的計算公式可知，最大的梯度幅值為：
因此設(shè)置nHist為1024足夠。以下實現(xiàn)統(tǒng)計直方圖：
- for(i=0;i<1024;i++)
- nHist[i] = 0;
- for(i=0; i<nHeight; i++)
- {
- for(j=0; j<nWidth; j++)
- {
- if(N[i*nWidth+j]==128)
- nHist[M[i*nWidth+j]]++;
- }
- }
3）獲取最大梯度幅值及潛在邊緣點個數(shù)
- nEdgeNum = nHist[0];
- nMaxMag = 0;
- for(i=1; i<1024; i++)
- {
- if(nHist[i] != 0)
- {
- nMaxMag = i;
- }
- nEdgeNum += nHist[i];
- }
4）計算兩個閾值
- double dRatHigh = 0.79;
- double dThrHigh;
- double dThrLow;
- double dRatLow = 0.5;
- nHighCount = (int)(dRatHigh * nEdgeNum + 0.5);
- j=1;
- nEdgeNum = nHist[1];
- while((j<(nMaxMag-1)) && (nEdgeNum < nHighCount))
- {
- j++;
- nEdgeNum += nHist[j];
- }
- dThrHigh = j;
- dThrLow = (int)((dThrHigh) * dRatLow + 0.5);
這段代碼的意思是，按照灰度值從低到高的順序，選取前79%個灰度值中的最大的灰度值為高閾值，低閾值大約為高閾值的一半。這是根據(jù)經(jīng)驗數(shù)據(jù)的來的，至于更好地參數(shù)選取方法，作者后面會另文研究。 5）進(jìn)行邊緣檢測
- SIZE sz;
- sz.cx = nWidth;
- sz.cy = nHeight;
- for(i=0; i<nHeight; i++)
- {
- for(j=0; j<nWidth; j++)
- {
- if((N[i*nWidth+j]==128) && (M[i*nWidth+j] >= dThrHigh))
- {
- N[i*nWidth+j] = 255;
- TraceEdge(i, j, dThrLow, N, M, sz);
- }
- }
- }
以上代碼在非極大值抑制產(chǎn)生的二值灰度矩陣的潛在點中按照高閾值尋找邊緣，并以所找到的點為中心尋找鄰域內(nèi)滿足低閾值的點，從而形成一個閉合的輪廓。然后對于不滿足條件的點，可用如下代碼直接刪除掉。
-
- for(i=0; i<nHeight; i++)
- {
- for(j=0; j<nWidth; j++)
- {
- if(N[i*nWidth+j] != 255)
- {
- N[i*nWidth+j] = 0 ;
- }
- }
- }
其中TraceEdge函數(shù)為一個嵌套函數(shù)，用于在每個像素點的鄰域內(nèi)尋找滿足條件的點。其實現(xiàn)代碼如下：
- void TraceEdge(int y, int x, int nThrLow, LPBYTE pResult, int *pMag, SIZE sz)
- {
-
- 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);
- }
- }
- }
以上就從原理上實現(xiàn)了整個Canny算法。其檢測效果如圖4所示。注意：以上代碼僅為作者理解所為，目的是驗證本人對算法的理解，暫時沒有考慮到代碼的執(zhí)行效率的問題。
圖4 邊緣檢測結(jié)果
4、擴(kuò)展
首先看一下OpenCV中cvCanny函數(shù)對該圖像的處理結(jié)果，如圖5所示。
圖5 OpenCV中的Canny邊緣檢測結(jié)果
對比圖4和圖5可以發(fā)現(xiàn)，作者自己實現(xiàn)的邊緣檢測效果沒有OpenCV的好，具體體現(xiàn)在：1）丟失了一些真的邊緣；2）增加了一些假的邊緣。 經(jīng)過對整個算法的來回檢查，初步推斷主要的問題可能在于在進(jìn)行灰度矩陣梯度幅值計算式所采用的模板算子性能不是太好，還有就是關(guān)于兩個閾值的選取方法。關(guān)于這兩個方面的改進(jìn)研究，后文闡述。
5、總結(jié)
本文是過去一段時間，對圖像邊緣檢測方法學(xué)習(xí)的總結(jié)。主要闡述了Canny算法的工作原理，實現(xiàn)過程，在此基礎(chǔ)上基于VC6.0實現(xiàn)了該算法，并給出了效果圖。最后，通過對比發(fā)現(xiàn)本文的實現(xiàn)方法雖然能夠?qū)崿F(xiàn)邊緣檢測，但效果還不是很理想，今后將在閾值選取原則和梯度幅值算子兩個方面進(jìn)行改進(jìn)。
|
