Sobel算子的理解

大家好，又见面了，我是你们的朋友全栈君。

sobel算子主要用于获得数字图像的一阶梯度，常见的应用和物理意义是边缘检测。

原理

算子使用两个33的矩阵（图1）算子使用两个33的矩阵（图1）去和原始图片作卷积，分别得到横向Ｇ（ｘ）和纵向Ｇ（ｙ）的梯度值，如果梯度值大于某一个阈值，则认为该点为边缘点

Gx方向的相关模板：

Gy方向的相关模板：

看了网上的很多资料，都把卷积和相关的概念给弄糊了，书上给的Sobel的模板不是卷积模板，而是相关模板，因为卷积的话要先将模板旋转180度以后再与图像做相关的操作。

所以Sobel的卷积模板是：

Gx=

具体计算如下：

图像的每一个像素的横向及纵向灰度值通过以下公式结合，来计算该点灰度的大小：

通常，为了提高效率使用不开平方的近似值：

然后可用以下公式计算梯度方向：

opencv还提供了一个scharr函数，比Sobel算子更为精准，也是3×3的模板。

下面放上源代码：

#include "opencv2/imgproc/imgproc.hpp"  #include "opencv2/highgui/highgui.hpp"  #include <iostream>  #include <cmath>using namespace cv;using namespace std;bool sobelEdge(Mat&  srcImage, Mat& resultImageX, Mat& resultImageY, uchar threshold){	CV_Assert(srcImage.channels() == 1);	// 初始化水平核因子	Mat sobelx = (Mat_<double>(3, 3) << -1, 0,		1, -2, 0, 2, -1, 0, 1);	// 初始化垂直核因子	Mat sobely = (Mat_<double>(3, 3) << -1, -2, -1,		0, 0, 0, 1, 2, 1);	resultImageX = Mat::zeros(srcImage.rows - 2,		srcImage.cols - 2, srcImage.type());	resultImageY = Mat::zeros(srcImage.rows - 2,		srcImage.cols - 2, srcImage.type());	double edgeX = 0;	double edgeY = 0;	double graMagX = 0;// 垂直方向上的梯度模长	double graMagY = 0;// 水平方向上的梯度模长	for (int k = 1; k < srcImage.rows - 1; ++k) 	{		for (int n = 1; n < srcImage.cols - 1; ++n) 		{			edgeX = 0;			edgeY = 0;			// 遍历计算水平与垂直梯度			for (int i = -1; i <= 1; ++i) 			{				for (int j = -1; j <= 1; ++j)				{					edgeX += srcImage.at<uchar>(k + i, n + j) *						sobelx.at<double>(1 + i, 1 + j);					edgeY += srcImage.at<uchar>(k + i, n + j) *						sobely.at<double>(1 + i, 1 + j);				}			}			// 计算垂直方向上的梯度模长			graMagX = sqrt(pow(edgeX, 2));			// 计算水平方向上的梯度模长			graMagY = sqrt(pow(edgeY, 2));			// 二值化			resultImageX.at<uchar>(k - 1, n - 1) =				((graMagX > threshold) ? 255 : 0);			// 二值化			resultImageY.at<uchar>(k - 1, n - 1) =				((graMagY > threshold) ? 255 : 0);		}	}	return true;}int OTSU(Mat &srcImage){	int nRows = srcImage.rows;	int nCols = srcImage.cols;	int threshold = 0;	double max = 0.0;	double AvePix[256];	int nSumPix[256];	double nProDis[256];	double nSumProDis[256];						 	for (int i = 0; i < 256; i++)	{		AvePix[i] = 0.0;		nSumPix[i] = 0;		nProDis[i] = 0.0;		nSumProDis[i] = 0.0;	}		for (int i = 0; i < nRows; i++)	{		for (int j = 0; j < nCols; j++)		{			nSumPix[(int)srcImage.at<uchar>(i, j)]++;		}	}		for (int i = 0; i < 256; i++)	{		nProDis[i] = (double)nSumPix[i] / (nRows*nCols);	}		AvePix[0] = 0;	nSumProDis[0] = nProDis[0];		for (int i = 1; i < 256; i++)	{		nSumProDis[i] = nSumProDis[i - 1] + nProDis[i];		AvePix[i] = AvePix[i - 1] + i*nProDis[i];	}	double mean = AvePix[255];	for (int k = 1; k < 256; k++)	{		double PA = nSumProDis[k];		double PB = 1 - nSumProDis[k];		double value = 0.0;		if (fabs(PA) > 0.001 && fabs(PB) > 0.001)		{			double MA = AvePix[k];//前一半的平均			double MB = (mean - PA*MA) / PB;//后一半的平均			value = value = (double)(PA * PB * pow((MA - MB), 2));//类间方差  		         //或者这样value = (double)(PA * PB * pow((MA-MB),2));//类间方差			 //pow(PA,1)* pow((MA - mean),2) + pow(PB,1)* pow((MB - mean),2)			if (value > max)			{				max = value;				threshold = k;			}		}	}	return threshold;}int main(){	Mat srcImage = cv::imread("building.jpg");	if (!srcImage.data)		return -1;	Mat srcGray;	cvtColor(srcImage, srcGray, CV_BGR2GRAY);	imshow("srcGray", srcGray);	//调用二值化函数得到最佳阈值	int otsuThreshold = OTSU(srcGray);	cout << otsuThreshold << endl;	Mat XresultImage;	Mat YresultImage;	sobelEdge(srcGray, XresultImage, YresultImage, otsuThreshold);	Mat resultImage;	//水平垂直边缘叠加	addWeighted(XresultImage, 0.5, YresultImage, 0.5, 0.0, resultImage);	imshow("resx", XresultImage);	imshow("resy", YresultImage);	imshow("res", resultImage);	waitKey(0);	return 0;}

原图灰度图：

垂直边缘图像：

水平边缘图像：

水平垂直合一起：

从上图可以看出 边缘图还有很多不是边缘点的也被取出来了，下一步就是做非极大值抑制