影象特徵提取：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;
}