使用索贝尔(Sobel)进行梯度运算时的数学意义和代码实现研究

对于做图像处理的工程师来说，Sobel非常熟悉且常用。但是当我们需要使用Sobel进行梯度运算，且希望得到“数学结果”（作为下一步运算的基础）而不是“图片效果”的时候，就必须深入了解Sobel的知识原理和OpenCV实现的细节（当然我们是OpenCV支持则）。这里对具体内容进行研究。

一、基本原理
一般来说，用来表示微分的最常用的算子是索贝尔(Sobel)算子，它可以实现任意阶导数和混合偏导数(例如: ∂2/∂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之外的获得。
感谢阅读至此，希望共同进步！