使用索贝尔(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之外的获得。

感谢阅读至此，希望共同进步！