canny算子_OpenCV圖像處理專欄十九 | 手動實作基于Canny算子的邊緣檢測

1. 前言

接着昨天手動構造Sobel算子實作檢測，今天來講講如何手動實作Canny邊緣檢測。由于要實作這個算法的需要的先驗知識比較多，是以在學習這個算法的實作之前我們先來學習一下用于圖像二值化的OSTU大津法。

2. OTSU算法原理

Ostu方法又名最大類間差方法，通過統計整個圖像的直方圖特性來實作全局門檻值

的自動選取。

它将圖像的像素按灰階級分成

個部分，使得

個部分之間的灰階值差異最小，通過方差的計算來尋找一個合适的灰階級别來劃分。

OTSU算法的計算很簡單，不受到圖像亮度和對比度的影響，是以使得類間方差最大的分割意味着錯分機率最小。

設

為灰階級門檻值，從

(一般為

)個灰階級周遊

，使得

為某個值的時候，前景和背景的方差最大，這個

值便是我們要求的灰階級門檻值。

設

表示分開後前景像素點占圖像像素比例，

表示分開後前景像素點的平均灰階，

表示分開後背景像素點數占圖像比例，

表示分開後背景像素點的平均灰階。

圖像的總灰階被定義為：

圖像的方差計算公式為：

并且有：

是以可以把公式(2)展開化簡為：

公式推導完畢，我們來看一下OTSU算法的詳細步驟。

3. OTSU算法流程

• 先計算圖像的直方圖，即将圖像所有的像素點按照

  

  共

  

  個bin，統計落在每個bin的像素點數量。
• 歸一化直方圖，也即将每個bin中像素點數量除以總的像素點。

• 

  表示分類的門檻值，也即一個灰階級，從

  

  開始疊代。
• 通過歸一化的直方圖，統計

  

  灰階級的像素(假設像素值在此範圍的像素叫做前景像素) 所占整幅圖像的比例

  

  ，并統計前景像素的平均灰階

  

  。然後統計

  

  灰階級的像素(假設像素值在此範圍的像素叫做背景像素) 所占整幅圖像的比例

  

  ，并統計背景像素的平均灰階

  

  。
• 計算前景像素和背景像素的方差

  

  （已經推導了）。
• 不斷

  

  ，轉到第

  

  個步驟。
• 将最大

  

  相應的

  

  值作為圖像的全局門檻值。

4. OTSU算法代碼實作

int OTSU(Mat src){
    int row = src.rows;
    int col = src.cols;
    int PixelCount[256]={0};
    float PixelPro[256]={0};
    int total_Pixel = row * col;
    float threshold = 0;
    //統計灰階級中每個像素在整副圖中的個數
    for(int i = 0; i < row; i++){
        for(int j = 0; j < col; j++){
            PixelCount[src.at<int>(i, j)]++;
        }
    }
    //計算每個像素在整副圖中的個數
    for(int i = 0; i < 256; i++){
        PixelPro[i] = (float)(PixelCount[i]) / (float)(total_Pixel);
    }
    //周遊灰階級[0,255]，計算出方差最大的灰階值，為最佳門檻值
    float w0, w1, u0tmp, u1tmp, u0, u1, u, deltaTmp, deltaMax = 0;
    for(int i = 0; i < 256; i++){
        w0 = w1 = u0tmp = u1tmp = u0 = u1 = u = deltaTmp = 0;
        for(int j = 0; j < 256; j++){
            if(j <= i){//以i為門檻值分類，第一類總的機率
                w0 += PixelPro[j];
                u0tmp += j * PixelPro[j];
            }else{//前景部分
                w1 += PixelPro[j];
                u1tmp += j * PixelPro[j];
            }
        }
        u0 = u0tmp / w0; //第一類的平均灰階
        u1 = u1tmp / w1; //第二類的平均灰階
        u = u0 + u1; //整副圖像的平均灰階
        //計算類間方差
        deltaTmp = w0 * (u0 - u) * (u0 - u) + w1 * (u1 - u) * (u1 - u);
        //找出最大類間方差以及對應門檻值
        if(deltaTmp > deltaMax){
            deltaMax = deltaTmp;
            threshold = i;
        }
    }
    return threshold;
}
           

5. 邊緣檢測的一般标準

邊緣檢測有下面幾個标準： (1) 以低的錯誤率檢測邊緣，也即意味着需要盡可能準确的捕獲圖像中盡可能多的邊緣。 (2) 檢測到的邊緣應精确定位在真實邊緣的中心。 (3) 圖像中給定的邊緣應隻被标記一次，并且在可能的情況下，圖像的噪聲不應産生假的邊緣。

6. Canny算子邊緣檢測步驟

有了上面的鋪墊，我們來到今天的主題，我們直接來看一下基于Canny算子進行邊緣檢測的步驟，我會在第6節詳細講解每一個步驟。基于Canny算子邊緣檢測的步驟如下：

• 使用高斯濾波算法，以平滑圖像，濾除噪聲
• 計算圖像中每個像素點的梯度強度和方向
• 應用非極大值抑制，以消除邊緣檢測帶來的雜散響應
• 應用雙門檻值檢測來确定真正的邊緣和潛在的邊緣
• 通過抑制孤立的弱邊緣最終完成邊緣檢測

7. 對應算法步驟的詳細解釋

7.1 高斯濾波

首先高斯函數的定義為

，其中

是圖像中像素點的坐标，

為标準差，高斯模闆就是利用這個函數來計算的。

接下來看高斯模闆是如何構造的，假設高斯模闆的尺寸為

。

為什麼長寬都為奇數，這主要是保證整個模闆有唯一中心元素，便于計算。

高斯模闆中的元素值為：

然後在實作生成高斯模闆時，又有兩種形式，即整數形式和小數形式。對于小數形式的就是按照上面的公式直接構造，不需要其他處理，而整數形式的需要歸一化，即将模闆左上角的值歸一化為

。

使用整數高斯模闆時，需要在模闆前面加一個系數，這個系數為

，就是模闆系數和的導數。

生成小數高斯模闆代碼如下：

#define PI 3.1415926
//生成小數形式的高斯模闆
void generateGaussianTemplate(double window[][11], int ksize, double sigma){
    int center = ksize / 2; //模闆的中心位置，也就是坐标原點
    double x2, y2;
    for(int i = 0; i < ksize; i++){
        x2 = pow(i - center, 2);
        for(int j = 0; j < ksize; j++){
            y2 = pow(j - center, 2);
            double g = exp(-(x2 + y2)) / (2 * sigma * sigma);
            g /= 2 * PI * sigma;
            window[i][j] = g;
        }
    }
    //歸一化左上角的數為1
    double k = 1 / window[0][0];
    for(int i = 0; i < ksize; i++) {
        for (int j = 0; j < ksize; j++) {
            window[i][j] *= k;
        }
    }
}
           

整數部分基本一樣，就不說了。對于

，

的高斯模闆結果為：

1.00000 2.71828 1.00000 
2.71828 7.38906 2.71828 
1.00000 2.71828 1.00000 
           

這裡想說一下

的作用，當

比較小的時候，生成的高斯模闆中心的系數比較大，而周圍的系數比較小，這樣對圖像的平滑效果不明顯。而當

比較大時，生成的模闆的各個系數相差就不是很大，比較類似于均值模闆，對圖像的平滑效果比較明顯。

7.2 計算梯度強度和方向

利用Sobel算子傳回水準

和垂直

的一階導數值。以此來計算梯度強度

和梯度方向。

關于Sobel算子的構造請看昨天的推文：

OpenCV圖像處理專欄十八 | 手動構造Sobel算子完成邊緣檢測

7.3 應用非極大值抑制，以消除邊緣檢測帶來的雜散響應

非極大值抑制是一種邊緣稀疏技術，作用在于

瘦邊

。在對圖檔進行梯度計算之後，僅僅基于梯度值提取的邊緣仍然比較模糊。根據

标準(3)

，對邊緣有且應當隻有一個準确的響應。而非極大值抑制可以幫助将将局部最大值之外的所有梯度值抑制為

，對梯度圖像中每個像素進行非極大值抑制的算法步驟為：

(1) 将目前像素的梯度強度與沿正負梯度方向上的兩個像素進行比較。

(2) 如果目前像素的梯度值與另外兩個像素值相比最大，則該像素點保留為邊緣點，否則該像素點被抑制。

通常為了更加精确的計算，

在跨梯度方向的兩個相鄰像素之間使用線性插值來得到要比較的像素梯度

。

如Fig1所示，将梯度分為

個方向，分别為E、NE、N、NW、W、SW、S、SE。其中

代表

,

代表

，

代表

，

代表

。像素點

的梯度方向為

，則像素點

和

的梯度線性插值為：

是以非極大值抑制的僞代碼如下：

需要注意的是，如何标志方向并不重要，重要的是梯度方向的計算要和梯度算子的選取保持一緻。

7.4 雙門檻值檢測

在施加非極大值抑制之後，剩餘的像素可以更準确地表示圖像中的實際邊緣。然而，仍然存在由于噪聲和顔色變化引起的一些邊緣像素。為了解決這些雜散響應，必須用弱梯度值過濾邊緣像素，并保留具有高梯度值的邊緣像素，可以通過選擇高低門檻值來實作。如果邊緣像素的梯度值高于高門檻值，則将其标記為強邊緣像素；如果邊緣像素的梯度值小于高門檻值并且大于低門檻值，則将其标記為弱邊緣像素；如果邊緣像素的梯度值小于低門檻值，則會被抑制。門檻值的選擇取決于給定輸入圖像的内容。

雙門檻值檢測的僞代碼如下：

7.5 抑制孤立弱邊緣完成邊緣檢測

到目前為止，被劃分為強邊緣的像素點已經被确定為邊緣，因為它們是從圖像中的真實 邊緣中提取出來的。然而，對于弱邊緣像素，将會有一些争論，因為這些像素可以從真實邊緣提取也可以是因噪聲或顔色變化引起的。為了獲得準确的結果，應該抑制由後者引起的弱邊緣。通常，由真實邊緣引起的弱邊緣像素将連接配接到強邊緣像素，而噪聲響應未連接配接。為了跟蹤邊緣連接配接，通過檢視弱邊緣像素及其

個鄰域像素，隻要其中一個為強邊緣像素，則該弱邊緣點就可以保留為真實的邊緣。

這部分的僞代碼如下：

8. C++代碼實作

按照上面的步驟一步步實作基于Canny的邊緣檢測算法，代碼如下：

const double PI = 3.1415926;

double getSum(Mat src){
    double sum = 0;
    for(int i = 0; i < src.rows; i++){
        for(int j = 0; j < src.cols; j++){
            sum += (double)src.at<double>(i, j);
        }
    }
    return sum;
}

Mat CannyEdgeDetection(cv::Mat src, int kSize, double hightThres, double lowThres){
//    if(src.channels() == 3){
//        cvtColor(src, src, CV_BGR2GRAY);
//    }
    cv::Rect rect;
    Mat gaussResult;
    int row = src.rows;
    int col = src.cols;
    printf("%d %dn", row, col);
    Mat filterImg = Mat::zeros(row, col, CV_64FC1);
    src.convertTo(src, CV_64FC1);
    Mat dst = Mat::zeros(row, col, CV_64FC1);
    int gaussCenter = kSize / 2;
    double  sigma = 1;
    Mat guassKernel = Mat::zeros(kSize, kSize, CV_64FC1);
    for(int i = 0; i < kSize; i++){
        for(int j = 0; j < kSize; j++){
            guassKernel.at<double>(i, j) = (1.0 / (2.0 * PI * sigma * sigma)) * (double)exp(-(((double)pow((i - (gaussCenter+1)), 2) + (double)pow((j-(gaussCenter+1)), 2)) / (2.0*sigma*sigma)));
        }
    }
    Scalar sumValueScalar = cv::sum(guassKernel);
    double sum = sumValueScalar.val[0];
    guassKernel = guassKernel / sum;
    for(int i = gaussCenter; i < row - gaussCenter; i++){
        for(int j = gaussCenter; j < col - gaussCenter; j++){
            rect.x = j - gaussCenter;
            rect.y = i - gaussCenter;
            rect.width = kSize;
            rect.height = kSize;
            //printf("%d %dn", i, j);
            //printf("%d %d %.5fn", i, j, cv::sum(guassKernel.mul(src(rect))).val[0]);
            filterImg.at<double>(i, j) = cv::sum(guassKernel.mul(src(rect))).val[0];
        }
    }
    Mat gradX = Mat::zeros(row, col, CV_64FC1); //水準梯度
    Mat gradY = Mat::zeros(row, col, CV_64FC1); //垂直梯度
    Mat grad = Mat::zeros(row, col, CV_64FC1); //梯度幅值
    Mat thead = Mat::zeros(row, col, CV_64FC1); //梯度角度
    Mat whereGrad = Mat::zeros(row, col, CV_64FC1);//區域
    //x方向的Sobel算子
    Mat Sx = (cv::Mat_<double>(3, 3) << -1, 0, 1, -2, 0, 2, -1, 0, 1);
    //y方向的Sobel算子
    Mat Sy = (cv::Mat_<double >(3, 3) << 1, 2, 1, 0, 0, 0, -1, -2, -1);
    //計算梯度指派和角度
    for(int i=1; i < row-1; i++){
        for(int j=1; j < col-1; j++){
            rect.x = j-1;
            rect.y = i-1;
            rect.width = 3;
            rect.height = 3;
            Mat rectImg = Mat::zeros(3, 3, CV_64FC1);
            filterImg(rect).copyTo(rectImg);
            //梯度和角度
            gradX.at<double>(i, j) += cv::sum(rectImg.mul(Sx)).val[0];
            gradY.at<double>(i, j) += cv::sum(rectImg.mul(Sy)).val[0];
            grad.at<double>(i, j) = sqrt(pow(gradX.at<double>(i, j), 2) + pow(gradY.at<double>(i, j), 2));
            thead.at<double>(i, j) = atan(gradY.at<double>(i, j)/gradX.at<double>(i, j));
            if(0 <= thead.at<double>(i, j) <= (PI/4.0)){
                whereGrad.at<double>(i, j) = 0;
            }else if(PI/4.0 < thead.at<double>(i, j) <= (PI/2.0)){
                whereGrad.at<double>(i, j) = 1;
            }else if(-PI/2.0 <= thead.at<double>(i, j) <= (-PI/4.0)){
                whereGrad.at<double>(i, j) = 2;
            }else if(-PI/4.0 < thead.at<double>(i, j) < 0){
                whereGrad.at<double>(i, j) = 3;
            }
        }
    }
    //printf("successn");
    //梯度歸一化
    double gradMax;
    cv::minMaxLoc(grad, &gradMax);
    if(gradMax != 0){
        grad = grad / gradMax;
    }
    //雙門檻值确定
    cv::Mat caculateValue = cv::Mat::zeros(row, col, CV_64FC1); //grad變成一維
    resize(grad, caculateValue, Size(1, grad.rows * grad.cols));
    cv::sort(caculateValue, caculateValue, CV_SORT_EVERY_COLUMN + CV_SORT_ASCENDING);//升序
    long long highIndex= row * col * hightThres;
    double highValue = caculateValue.at<double>(highIndex, 0); //最大門檻值
    double lowValue = highValue * lowThres;
    //NMS
    for(int i = 1 ; i < row-1; i++ ){
        for( int j = 1; j<col-1; j++){
            // 八個方位
            double N = grad.at<double>(i-1, j);
            double NE = grad.at<double>(i-1, j+1);
            double E = grad.at<double>(i, j+1);
            double SE = grad.at<double>(i+1, j+1);
            double S = grad.at<double>(i+1, j);
            double SW = grad.at<double>(i-1, j-1);
            double W = grad.at<double>(i, j-1);
            double NW = grad.at<double>(i -1, j -1); // 區域判斷，線性插值處理
            double tanThead; // tan角度
            double Gp1; // 兩個方向的梯度強度
            double Gp2; // 求角度，絕對值
            tanThead = abs(tan(thead.at<double>(i,j)));
            switch ((int)whereGrad.at<double>(i,j)) {
                case 0: Gp1 = (1- tanThead) * E + tanThead * NE; Gp2 = (1- tanThead) * W + tanThead * SW; break;
                case 1: Gp1 = (1- tanThead) * N + tanThead * NE; Gp2 = (1- tanThead) * S + tanThead * SW; break;
                case 2: Gp1 = (1- tanThead) * N + tanThead * NW; Gp2 = (1- tanThead) * S + tanThead * SE; break;
                case 3: Gp1 = (1- tanThead) * W + tanThead *NW; Gp2 = (1- tanThead) * E + tanThead *SE; break;
                default: break;
            }
            // NMS -非極大值抑制和雙門檻值檢測
            if(grad.at<double>(i, j) >= Gp1 && grad.at<double>(i, j) >= Gp2){
                //雙門檻值檢測
                if(grad.at<double>(i, j) >= highValue){
                    grad.at<double>(i, j) = highValue;
                    dst.at<double>(i, j) = 255;
                } else if(grad.at<double>(i, j) < lowValue){
                    grad.at<double>(i, j) = 0;
                } else{
                    grad.at<double>(i, j) = lowValue;
                }
            } else{
                grad.at<double>(i, j) = 0;
            }
        }
    }
    //抑制孤立低門檻值點 3*3. 找到高門檻值就255
    for(int i = 1; i < row - 1; i++){
        for(int j = 1; j < col - 1; j++){
            if(grad.at<double>(i, j) == lowValue){
                //3*3 區域強度
                rect.x = i-1;
                rect.y = j-1;
                rect.width = 3;
                rect.height = 3;
                for(int x = 0; x < 3; x++){
                    for(int y = 0; y < 3; y++){
                        if(grad(rect).at<double>(x, y)==highValue){
                            dst.at<double>(i, j) = 255;
                            break;
                        }
                    }
                }
            }
        }
    }
    return dst;
}

           

9. 效果

邊緣檢測的效果不如OpenCV自帶的，需要再研究研究。

10. 結論

本文介紹了圖像處理算法中最常用的邊緣檢測算法的原理以及一個C++複現，然而可惜的是效果沒有OpenCV自帶算法的效果好，并且速度也會稍慢，一定是某個細節沒有處理完美，但作為科普來講大概已經夠用了，如果想實作完美的Sobel邊緣檢測請檢視OpenCV源碼或者和我一起讨論。

11. 參考

• https://blog.csdn.net/jmu201521121021/article/details/80622011
• https://www.cnblogs.com/techyan1990/p/7291771.html

歡迎關注GiantPandaCV, 在這裡你将看到獨家的深度學習分享，堅持原創，每天分享我們學習到的新鮮知識。( • ̀ω•́ )✧

有對文章相關的問題，或者想要加入交流群，歡迎添加BBuf微信：

https://u.wechat.com/MPWFDnmCPu6zgf5YUtdpT_U (二維碼自動識别)