圖像特征提取：Sobel邊緣檢測

前言

點和線是做圖像分析時兩個最重要的特征，而線條往往反映了物體的輪廓，對圖像中邊緣線的檢測是圖像分割與特征提取的基礎。文章主要討論兩個實際工程中常用的邊緣檢測算法：Sobel邊緣檢測和Canny邊緣檢測，Canny邊緣檢測由於算法復雜將在另一篇文章中單獨介紹，文章不涉及太多原理，因為大部分的圖像處理書籍都有相關內容介紹，文章主要通過Matlab代碼，一步一步具體實現兩種經典的邊緣檢測算法。

Sobel邊緣檢測

Soble邊緣檢測算法比較簡，實際應用中效率比canny邊緣檢測效率要高，但是邊緣不如Canny檢測的准確，但是很多實際應用的場合，sobel邊緣卻是首選，尤其是對效率要求較高，而對細紋理不太關心的時候。

Soble邊緣檢測通常帶有方向性，可以只檢測豎直邊緣或垂直邊緣或都檢測。

所以我們先定義兩個梯度方向的系數：

kx=0;ky=1;% horizontal
kx=1;ky=0;% vertical
kx=1;ky=1;% both

然后我們來計算梯度圖像，我們知道邊緣點其實就是圖像中灰度跳變劇烈的點，所以先計算梯度圖像，然后將梯度圖像中較亮的那一部分提取出來就是簡單的邊緣部分。

Sobel算子用了一個3*3的濾波器來對圖像進行濾波從而得到梯度圖像，這里面不再詳細描述怎樣進行濾波及它們的意義等。

豎起方向的濾波器：y_mask=op = [-1 -2 -1;0 0 0;1 2 1]/8;

水平方向的濾波器：op的轉置：x_mask=op’;

定義好濾波器后，我們就開始分別求垂直和豎起方向上的梯度圖像。用濾波器與圖像進行卷積即可：

bx = abs(filter2(x_mask,a));
by = abs(filter2(y_mask,a));

上面bx為水平方向上的梯度圖像，by為垂直方向上的梯度圖像。為了更清楚的說明算法過程，下面給出一張示例圖像的梯度圖像。

原圖：

豎直方向梯度圖像：by

水平方向梯度圖像：bx

而最終的梯度圖像為：b = kx*bx.*bx + ky*by.*by;當然這里如果定義了檢測方向，就會得到對應上面的圖像。

這里值得注意的是為了計算效率並沒有對b開平方。而涉及濾波等操作時對圖像邊界的處理是值得注意的一個地方。我們這里應該將梯度圖像的四周1像素點都設置為0。

得到梯度圖像后，我們需要的是計算閾值，這是Sobel算法很核心的一部分，也是對效果影響較大的地方，同理講到canny邊緣檢測時，用到的雙閾值法也是canny算法的核心。

同樣這里，我並不太多的介紹算法原理，相關文獻可以直接百度。閾值也可以通過自己期待的效果進行自定義閾值，如果沒有，則我們計算默認閾值。

scale = 4;
cutoff = scale*mean2(b);
thresh = sqrt(cutoff);

其中mean2函數是求圖像所有點灰度的平均值。scale是一個系數。

有了閾值后，我們就可以地得到的梯度圖像進行二值化。二值化過程，不做詳細說明，遍歷圖像中的像素點，灰度大於閾值的置為白點，灰度小於閾值的則置為黑點。

下面是示例圖像梯度圖像二值化后的效果：

其實很多介紹Soble算法的文章介紹到這里就結束了，印象中OpenCv的算法也是到此步為止。但是我們注意到上面的邊緣圖像，邊緣較粗，如果我們在做邊界跟蹤或輪廓提取時，上面圖像並不是我們期望的結果。

所以下面要介紹一個很重要的算法，用非極大值抑制算法對邊緣進行1像素化。本人在開始的時候也一直以為這里用一個普通的細化算法就可以了，可是總得到到想要的結果，最后查找matlab早期版本的源碼才找到方法，其實跟canny算法里原理差不多。

我們需要遍歷剛才得到的梯度圖像二值化結果b,對於b內的任意一點(i,j)，如果滿足下面條件，則保持白點，否則置為黑點。條件簡單的說即是：如果是豎直邊緣，則它的梯度值要比左邊和右邊的點都大；如果是水平連續，則該點的梯度值要比上邊和下邊的都大。

bx(i,j)>kx*by(i,j) && b(i,j-1)<b(i,j) && b(i,j+1)<b(i,j)

或者

by(i,j)>ky*bx(i,j) && b(i-1,j)<b(i,j) &&b (i+1,j)<b(i,j)

經過這樣的非極大值抑制我們就可以得到比較理想的邊緣圖像。

下面給出利用OpenCV里的一些濾波函數，從新寫的一個Sobel邊緣檢測的函數：

bool Sobel(const Mat& image,Mat& result,int TYPE)
{
    if(image.channels()!=1)
        return false;
    // 系數設置
    int kx(0);
    int ky(0);
    if( TYPE==SOBEL_HORZ ){
        kx=0;ky=1;
    }
    else if( TYPE==SOBEL_VERT ){
        kx=1;ky=0;
    }
    else if( TYPE==SOBEL_BOTH ){
        kx=1;ky=1;
    }
    else
        return false;

    // 設置mask
    float mask[3][3]={{1,2,1},{0,0,0},{-1,-2,-1}};
    Mat y_mask=Mat(3,3,CV_32F,mask)/8;
    Mat x_mask=y_mask.t(); // 轉置

    // 計算x方向和y方向上的濾波
    Mat sobelX,sobelY;
    filter2D(image,sobelX,CV_32F,x_mask);
    filter2D(image,sobelY,CV_32F,y_mask);
    sobelX=abs(sobelX);
    sobelY=abs(sobelY);
    // 梯度圖
    Mat gradient=kx*sobelX.mul(sobelX)+ky*sobelY.mul(sobelY);

    // 計算閾值
    int scale=4;
    double cutoff=scale*mean(gradient)[0];

    result.create(image.size(),image.type());
    result.setTo(0);
    for(int i=1;i<image.rows-1;i++)
    {
        float* sbxPtr=sobelX.ptr<float>(i);
        float* sbyPtr=sobelY.ptr<float>(i);
        float* prePtr=gradient.ptr<float>(i-1);
        float* curPtr=gradient.ptr<float>(i);
        float* lstPtr=gradient.ptr<float>(i+1);
        uchar* rstPtr=result.ptr<uchar>(i);
        // 閾值化和極大值抑制
        for(int j=1;j<image.cols-1;j++)
        {
            if( curPtr[j]>cutoff && (
                (sbxPtr[j]>kx*sbyPtr[j] && curPtr[j]>curPtr[j-1] && curPtr[j]>curPtr[j+1]) ||
                (sbyPtr[j]>ky*sbxPtr[j] && curPtr[j]>prePtr[j] && curPtr[j]>lstPtr[j]) ))
                rstPtr[j]=255;
        }
    }

    return true;
}