OpenCV学习笔记（九）——Sobel边缘检测

前言：

      Sobel算子是像素图像边缘检测中最重要的算子之一，在机器学习、数字媒体、计算机视觉等信息科技领域起着举足轻重的作用。在技术上，它是一个离散的一阶差分算子，用来计算图像亮度函数的一阶梯度之近似值。在图像的任何一点使用此算子，将会产生该点对应的梯度矢量或是其法矢量。

在边缘检测中，常用的一种模板是Sobel 算子。Sobel 算子有两个，一个是检测水平边缘的 ；另一个是检测垂直边缘的 。与Prewitt算子相比，Sobel算子对于像素的位置的影响做了加权，可以降低边缘模糊程度，因此效果更好。 由于Sobel算子是滤波算子的形式，用于提取边缘，可以利用快速 卷积 函数， 简单有效，因此应用广泛。美中不足的是，Sobel算子并没有将图像的主体与背景严格地区分开来，换言之就是Sobel算子没有基于 图像灰度 进行处理，由于Sobel算子没有严格地模拟人的视觉生理特征，所以提取的图像轮廓有时并不能令人满意。

一、sobel算子的性质及构建方法

1.1 sobel算子的分离性

        我们先来看一个3阶的Sobel边缘检测算子，如下所示：

    

      显然，3*3的Sobel算子是可分离的，它是Sobel算子的标准形式，可以利用二项式展开式的系数构建窗口更大的Sobel算子，如5*5、7*7等，但是有一点必须要注意，窗口大小要为奇数。

      那么怎么来构建任意大小的sobel算子呢？

 1.2  sobel算子的构建过程

       Sobel算子是在一个坐标轴的方向上进行非归一化的高斯平滑，在另外一个坐标轴方向上进行差分处理，

的Sobel算子是由平滑算子和差分算子full卷积而得到的，其中为奇数。对于窗口大小为的非归一化的

Sobel平滑算子等于阶的二项式展开式的系数，那么问题只剩下怎么构建窗口大小为的Sobel差分算子？其

实，窗口大小为的Sobel差分算子是在阶的二项式展开式的系数两侧补零，然后向后差分得到的。

      接下来我们来构建一个4阶的非均一化的Sobel平滑算子（其实就是高斯平滑算子）和Sobel差分算子，来理解整个构建过程。

      第一步：取二项式的指数n=3，然后计算展开式的系数，如下所示：

也就是：

这就是4阶的非均一化的Sobel平滑算子（其实就是高斯平滑算子）。

      第二步： 取二项式的指数，然后计算展开式的系数，即：

也就是：

然后在两侧补零，得到：

接着后向差分（后面的数值减去前面相邻的数值），即得到差分后的结果为：

这就是4阶的Sobel差分算子。

        第三步：将4阶的Sobel平滑算子（其实就是高斯平滑算子）和Sobel差分算子进行full卷积，即可得到的Sobel算子，即： 

Sobel平滑算子和差分算子总结如下：

二、代码实现

    

#include <opencv2/opencv.hpp>
#include <iostream>
#include <math.h>

using namespace cv;
using namespace std;




/*********************************************************************************************************/
/*     用sobel算子对灰度图像进行滤波           */

int factorial(int n)
{
	// factorial()函数实现阶乘
	int fac = 1;
	if (n == 0)
	{
		return fac;
	}
	for (int i = 1; i <= n; ++i)
	{
		fac *= i;
	}
	return fac;
}

Mat getPascalSmooth(int n)
{
	//  getPascalSmooth()函数用来创建sobel平滑算子
	Mat pascalSmooth = Mat::zeros(Size(n, 1), CV_32FC1);
	for (int i = 0; i < n; ++i)
	{
		pascalSmooth.at<float>(0, i) =float( factorial(n - 1) / (factorial(i)*factorial(n - 1 - i)));

	}
	return pascalSmooth;
}

Mat getPascalDiff(int n)
{
	//  getPascalDiff()函数用来创建sobel差分算子
	Mat pascalDiff = Mat::zeros(Size(n, 1), CV_32FC1);
	Mat pascalSmooth_previous = getPascalSmooth(n - 1);
	for (int i = 0; i < n; ++i)
	{
		if (i == 0)
		{
			pascalDiff.at<float>(0, i) = 1;
		}
		else if (i == n - 1)
		{
			pascalDiff.at<float>(0, i) = -1;
		}
		else
		{
			pascalDiff.at<float>(0, i) = pascalSmooth_previous.at<float>(0, i) - pascalSmooth_previous.at<float>(0, i - 1);
		}
	}
	return  pascalDiff;
}


void conv2D(InputArray src, InputArray kernel, OutputArray dst, int ddepth, Point anchor = Point(-1, -1), int borderType = BORDER_DEFAULT)
{
	//  conv2D()函数用来完成same卷积运算
	Mat kernelFlip;
	flip(kernel, kernelFlip, -1);    // 卷积运算的第一步，将卷积核逆时针翻转180°
	filter2D(src, dst, ddepth, kernelFlip, anchor, 0.0, borderType);   //卷积运算的第二步
}

void sepConv2D_Y_X(InputArray src, OutputArray  src_kerY_kerX, int ddepth, InputArray kernelY, InputArray kernelX, Point anchor = Point(-1, -1), int borderType = BORDER_DEFAULT)
{
	// 对于可分离的离散二维卷积，先进行垂直方向上的卷积，再进行水平方向上的卷积
	Mat src_kerY;
	conv2D(src, kernelY, src_kerY, ddepth, anchor, borderType);   // 输入矩阵与垂直方向上的卷积核的卷积
	conv2D(src_kerY, kernelX, src_kerY_kerX, ddepth, anchor, borderType);  // 把从上面得到的卷积结果和水平方向上的卷积核卷积
}

void sepConv2D_X_Y(InputArray src, OutputArray  src_kerX_kerY, int ddepth, InputArray kernelX, InputArray kernelY, Point anchor = Point(-1, -1), int borderType = BORDER_DEFAULT)
{
	// 对于可分离的离散二维卷积，先进行水平方向上的卷积，再进行垂直方向上的卷积
	Mat src_kerX;
	conv2D(src, kernelX, src_kerX, ddepth, anchor, borderType);   // 输入矩阵与水平方向上的卷积核的卷积
	conv2D(src_kerX, kernelY, src_kerX_kerY, ddepth, anchor, borderType);  // 把从上面得到的卷积结果和垂直方向上的卷积核卷积
}



Mat sobel(Mat& image, int x_flag,int y_flag,int winSize, int borderType)
{
	// sobel函数用来完成图像灰度矩阵与sobel核的卷积 
	CV_Assert(winSize >= 3 && winSize % 2 == 1);   // sobel卷积核的窗口大小为大于3的奇数
	Mat pascalSmooth = getPascalSmooth(winSize);   // sobel平滑算子
	Mat pascalDiff = getPascalSmooth(winSize);     // sobel差分算子
	Mat image_con_sobel;                           // 输出矩阵
	if (x_flag != 0)
	{

		sepConv2D_Y_X(image, image_con_sobel, CV_32FC1, pascalSmooth.t(), pascalDiff, Point(-1, -1), borderType);     //先进行一维垂直方向上的平滑，再进行一维水平方向上的差分，即图像与sobel_x进行卷积运算
	}
	if (x_flag == 0 && y_flag != 0)
	{
		sepConv2D_X_Y(image, image_con_sobel, CV_32FC1, pascalSmooth, pascalDiff.t(), Point(-1, -1), borderType);    //先进行一维水平方向上的平滑，再进行一维垂直方向上的差分,即图像与sobel_y进行卷积运算
	}
	return  image_con_sobel;
}
/********************************************************************************************************************/




//  主函数
int main( )
{
	const Mat src_img = imread("test11.png");
	if (src_img.empty())
	{
		printf("could not load image...\n");
		return -1;
	}
	namedWindow("原图：", CV_WINDOW_AUTOSIZE);
	imshow("原图：", src_img);
	
	// 将彩色图转化为灰度图，调用OpenCV提供的cvtColor接口
	Mat gray_img;
	cvtColor(src_img, gray_img,CV_BGR2GRAY);
	namedWindow("灰度图", CV_WINDOW_AUTOSIZE);
	imshow("灰度图", gray_img);

	// 用sobel算子计算度图像的水平和垂直方向上的差分
	Mat image_con_sobel_Ix, image_con_sobel_Iy;
	image_con_sobel_Ix = sobel(gray_img, 1, 0, 3, BORDER_DEFAULT);
	image_con_sobel_Iy = sobel(gray_img, 0, 1, 3, BORDER_DEFAULT);
	
	//  水平方向和垂直方向上的边缘强度
	// 数据类型转换，边缘强度的灰度级显示
	Mat scale_sobel_Ix, scale_sobel_Iy;
	convertScaleAbs(image_con_sobel_Ix, scale_sobel_Ix);              //  转化为8位灰度级显示
	convertScaleAbs(image_con_sobel_Iy, scale_sobel_Iy);
	imshow("垂直方向的边缘", scale_sobel_Ix);
	imshow("水平方向的边缘", scale_sobel_Iy);
	waitKey(0);
	return 0;
}