使用索貝爾(Sobel)進行梯度運算時的數學意義和代碼實現研究

對於做圖像處理的工程師來說，Sobel非常熟悉且常用。但是當我們需要使用Sobel進行梯度運算，且希望得到“數學結果”（作為下一步運算的基礎）而不是“圖片效果”的時候，就必須深入了解Sobel的知識原理和OpenCV實現的細節（當然我們是OpenCV支持則）。這里對具體內容進行研究。

一、基本原理

一般來說，用來表示微分的最常用的算子是索貝爾(Sobel)算子，它可以實現任意階導數和混合偏導數(例如: ∂²/∂x∂y)。

在x方向上用Sobel算子進行近似一階求導的結果

Sobel算子效果，y方向近似一階導數

OpenCV中給出了函數使用的定義

void cv::Sobel(
  cv::InputArray  src,                 // 源圖像
  cv::OutputArray dst,                 // 目標圖像
  int             ddepth,              // 像素深度 (如CV_8U)
  int             xorder,              // x方向對應的倒數順序
  int             yorder,              // y方向對應的倒數順序
  cv::Size        ksize      = 3,      // 核大小
  double          scale      = 1,      // 閾值
  double          delta      = 0,      // 偏移
  int             borderType = cv::BORDER_DEFAULT  // 邊框外推方法
);

1、其中src和dst是源圖像和目標圖像，可以通過指明參數ddepth來確定目標圖像的深度或類型（如CV_32F）。舉個例子，如果src是一幅8位圖像，那么dst需要至少CV_16S的深度保證不出現溢出；

2、xorder和yorder是求導順序，其取值范圍為0、1和2。0表示在這個方向上不進行求導，那2代表什么？

3、ksize是一個奇數，表示了調用的濾波器的寬和高，目前最大支持到31；

4、閾值和偏移將在把結果存入dst前調用，這有助於你將求導結果可視化.borderType參數的功能與其他卷積操作完全一樣。

上面有一個遺留問題，就是xorder(yorder)取2的時候代表什么？為此翻閱OpenCV源碼

step1

 

 

step 2

step3

 

static void getSobelKernels( OutputArray _kx, OutputArray _ky,
                             int dx, int dy, int _ksize, bool normalize, int ktype )
{
    int i, j, ksizeX = _ksize, ksizeY = _ksize;
    if( ksizeX == 1 && dx > 0 )
        ksizeX = 3;
    if( ksizeY == 1 && dy > 0 )
        ksizeY = 3;
    CV_Assert( ktype == CV_32F || ktype == CV_64F );
    _kx.create(ksizeX, 1, ktype, -1, true);
    _ky.create(ksizeY, 1, ktype, -1, true);
    Mat kx = _kx.getMat();
    Mat ky = _ky.getMat();
    if( _ksize % 2 == 0 || _ksize > 31 )
        CV_Error( CV_StsOutOfRange, "The kernel size must be odd and not larger than 31" );
    std::vector<int> kerI(std::max(ksizeX, ksizeY) + 1);
    CV_Assert( dx >= 0 && dy >= 0 && dx+dy > 0 );
    for( int k = 0; k < 2; k++ )
    {
        Mat* kernel = k == 0 ? &kx : &ky;
        int order = k == 0 ? dx : dy;
        int ksize = k == 0 ? ksizeX : ksizeY;
        CV_Assert( ksize > order );
        if( ksize == 1 )
            kerI[0] = 1;
        else if( ksize == 3 )
        {
            if( order == 0 )
                kerI[0] = 1, kerI[1] = 2, kerI[2] = 1;
            else if( order == 1 )
                kerI[0] = -1, kerI[1] = 0, kerI[2] = 1;
            else
                kerI[0] = 1, kerI[1] = -2, kerI[2] = 1;
        }
        else
        {
            int oldval, newval;
            kerI[0] = 1;
            for( i = 0; i < ksize; i++ )
                kerI[i+1] = 0;
            for( i = 0; i < ksize - order - 1; i++ )
            {
                oldval = kerI[0];
                for( j = 1; j <= ksize; j++ )
                {
                    newval = kerI[j]+kerI[j-1];
                    kerI[j-1] = oldval;
                    oldval = newval;
                }
            }
            for( i = 0; i < order; i++ )
            {
                oldval = -kerI[0];
                for( j = 1; j <= ksize; j++ )
                {
                    newval = kerI[j-1] - kerI[j];
                    kerI[j-1] = oldval;
                    oldval = newval;
                }
            }
        }
        Mat temp(kernel->rows, kernel->cols, CV_32S, &kerI[0]);
        double scale = !normalize ? 1. : 1./(1 << (ksize-order-1));
        temp.convertTo(*kernel, ktype, scale);
    }
}

其中

 

那么，可以看見，當order ==2 時候，生成了[1 -2 1]作為類似的模板，不管是什么，這個不是我想要的。

     除了上面看到的，還可以發現同時設置xorder 和 yorder的時候，最后並沒有看到相加的動作。而如果我們計算的結果是梯度場的時候，就不僅要算xorder，而且要算yorder,並且最后要把這兩個結果求和。如果自己編碼，那么可能如下：

//求x方向偏導
   vx = *(lpSrc + IMGW + 1) - *(lpSrc + IMGW - 1) +
   *(lpSrc + 1)*2 - *(lpSrc - 1)*2 +
   *(lpSrc - IMGW + 1) - *(lpSrc - IMGW - 1);
   //求y方向偏導
   vy = *(lpSrc + IMGW - 1) - *(lpSrc - IMGW - 1) +
   *(lpSrc + IMGW)*2 - *(lpSrc - IMGW)*2 +
   *(lpSrc + IMGW + 1) - *(lpSrc - IMGW + 1);
   gradSum += (abs(vx)+abs(vy));        

 

二、定量研究

    如果直接用自然圖片（比如lena）,sobel計算之后可能啥都看不出來。為了定量研究，必須自己做圖片。

    搞成這樣，看的清楚

//自己生成實驗圖片
    Mat matTst  =  Mat(Size( 11 , 11 ),CV_8UC1,Scalar( 0 ));
    line(matTst,Point( 5 , 0 ),Point( 5 , 11 ),Scalar( 255 ));
    line(matTst,Point( 0 , 5 ),Point( 11 , 5 ),Scalar( 255 ));

求一遍X方向的導數

//x方向求導
    Sobel(matTst,matX,CV_16SC1,1,0);

為什么要用16s?首先，所有的模式為CV_{8U,16S,16U,32S,32F,64F}C{1,2,3}，其中U代表char,S代表short,F代表float，這個和c++里面一樣；8只能是正數，16/32都可以是負數。所有的C1和C都是一樣的，比如cv_8uc1 == cv_8u。那么原始圖像是CV_8UC1，也就是CV_8U了,那么進行卷積計算，也就是原圖像和

-1 0 1

-2 0 2

-1 0 1

或者類似的東西進行卷積計算，那么結果得到最大為 1*255+2*255+1*255 - 0  > 10000，最小的結果為 0 - 1*255+2*255+1*255  < -10000,所以要用CV_16SC1來保存，才合適。那么得到的結果為

​首先是有正有負，集中在中間區域值變化非常的大。取其中一個像素來看

中心位置 使用

-1 0 1

-2 0 2

-1 0 1

進行卷積，那么結果為

0*-1 + 255*0 + 0*1 +255*-2 +255*0 + 255*2 + 0 *-1 +255*0+0*1 = 0+0+0-255+0+255*2+0+0+0 = 0

其他結果也是，那么結果是符合卷積運算的。

計算Y方向卷積

//y方向求導
    Sobel(matTst,matY,CV_16SC1,0,1);

  

計算XY的卷積

   //xy方向
    Sobel(matTst,matXY,CV_16SC1,1,1);

這個結果不是我想要的，我想要的是 (abs(vx)+abs(vy));    

那么這樣實現

//求和

    matXY = abs(matX) + abs(matY);

甚至可以進一步幫助顯示

//方便顯示

 normalize(matXY,matXY,0,255,NORM_MINMAX);

    matXY.convertTo(matXY,CV_8UC1);

無論如何，這個結果和直接 Sobel(matTst,matXY,CV_16SC1, 1 , 1 );都是不一樣的

三、小結反思

應該說，Sobel大概能夠做什么？這個很早之前就已經知道了。但是為什么能夠達到這樣的效果？這個問題，一直到我需要使用Sobel進行相關的數學計算的時候才能夠搞明白。

掌握知識最后要歸結到數學抽象層次，才能夠算是徹底掌握。這是Sobel之外的獲得。

感謝閱讀至此，希望共同進步！