图像特征提取：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边缘检测的函数：
1 bool Sobel(const Mat& image,Mat& result,int TYPE)
2 {
3 if(image.channels()!=1)
4 return false;
5 // 系数设置
6 int kx(0);
7 int ky(0);
8 if( TYPE==SOBEL_HORZ ){
9 kx=0;ky=1;
10 }
11 else if( TYPE==SOBEL_VERT ){
12 kx=1;ky=0;
13 }
14 else if( TYPE==SOBEL_BOTH ){
15 kx=1;ky=1;
16 }
17 else
18 return false;
19
20 // 设置mask
21 float mask[3][3]={{1,2,1},{0,0,0},{-1,-2,-1}};
22 Mat y_mask=Mat(3,3,CV_32F,mask)/8;
23 Mat x_mask=y_mask.t(); // 转置
24
25 // 计算x方向和y方向上的滤波
26 Mat sobelX,sobelY;
27 filter2D(image,sobelX,CV_32F,x_mask);
28 filter2D(image,sobelY,CV_32F,y_mask);
29 sobelX=abs(sobelX);
30 sobelY=abs(sobelY);
31 // 梯度图
32 Mat gradient=kx*sobelX.mul(sobelX)+ky*sobelY.mul(sobelY);
33
34 // 计算阈值
35 int scale=4;
36 double cutoff=scale*mean(gradient)[0];
37
38 result.create(image.size(),image.type());
39 result.setTo(0);
40 for(int i=1;i<image.rows-1;i++)
41 {
42 float* sbxPtr=sobelX.ptr<float>(i);
43 float* sbyPtr=sobelY.ptr<float>(i);
44 float* prePtr=gradient.ptr<float>(i-1);
45 float* curPtr=gradient.ptr<float>(i);
46 float* lstPtr=gradient.ptr<float>(i+1);
47 uchar* rstPtr=result.ptr<uchar>(i);
48 // 阈值化和极大值抑制
49 for(int j=1;j<image.cols-1;j++)
50 {
51 if( curPtr[j]>cutoff && (
52 (sbxPtr[j]>kx*sbyPtr[j] && curPtr[j]>curPtr[j-1] && curPtr[j]>curPtr[j+1]) ||
53 (sbyPtr[j]>ky*sbxPtr[j] && curPtr[j]>prePtr[j] && curPtr[j]>lstPtr[j]) ))
54 rstPtr[j]=255;
55 }
56 }
57
58 return true;
59 }