图像梯度之Sobel算子，及其Matlab和OpenCV实现

1、Sobel算子实现步骤
（1）需要用到计算$x,y$方向导数值的两个滤波核

（2）分别与图像进行卷积

（3）把得到的值卸载一起得到图像梯度向量

梯度方向为：$\theta ={\tan }^{-1}\left(\frac{\partial f}{\partial x}/\frac{\partial f}{\partial y}\right)$
梯度幅度为：$\Vert \nabla f\Vert =\sqrt{{\left(\frac{\partial f}{\partial x}\right)}^2+{\left(\frac{\partial f}{\partial y}\right)}^2}$
通常为了提高效率，使用不开平方的近似值：

如果梯度大于某一阈值，则认为该点$(x,y)$为边缘点。
2、Matlab实现
（1）for循环遍历

img=imread("lena.png");
[x,y,z]=size(img);
img=im2double(img); % 数据类型转化uint8-->double，uint8类型参与运算容易溢出
if(z==3)    % rgb-->gray
   img_gray=(img(:,:,1).*229+img(:,:,2).*578+img(:,:,3).*114)./1000;
end
%% Sobel算子
img_sobel=zeros(x,y);
img_sobelx=img_sobel;
img_sobely=img_sobel;
for i=2:x-1
    for j=2:y-1
       img_sobelx(i,j)=abs(img_gray(i-1,j+1)-img_gray(i-1,j-1)+2*img_gray(i,j+1)-2*img_gray(i,j-1)+img_gray(i+1,j+1)-img_gray(i+1,j-1));
       img_sobely(i,j)=abs(img_gray(i-1,j-1)-img_gray(i+1,j-1)+2*img_gray(i-1,j)-2*img_gray(i+1,j)+img_gray(i-1,j+1)-img_gray(i+1,j+1));
    end    
end
img_sobel=img_sobelx+img_sobely;
for i=1:x
    for j=1:y
        if img_sobel(i,j)>0.7
           img_sobel(i,j)=1;
        else
           img_sobel(i,j)=0;
        end
    end    
end
figure(1)
imshow(img)
figure(2)
imshow(img_sobel)

结果如图所示：

注意：上述程序的实现结果与Matlab函数edge()实现Sobel算子的结果还是有较大的出入，Matlab函数做了进一步的处理。结果如图所示：

（2）利用卷积函数conv2()

sobelx=[-1,0,1;-2,0,2;-1,0,1];	% 卷积核
sobely=[1,2,1;0,0,0;-1,-2,-1];
img_sobelx=conv2(img_gray,sobelx);
img_sobely=conv2(img_gray,sobely);
img_sobel=abs(img_sobelx)+abs(img_sobely);
img_sobel=img_sobel(2:end-1,2:end-1);
for i=1:x
    for j=1:y
        if img_sobel(i,j)>0.7
           img_sobel(i,j)=1;
        else
           img_sobel(i,j)=0;
        end
    end    
end
figure(4)
imshow(img_sobel)

利用卷积函数conv2()和for循环求解出的结果是一样的。如图所示：

3、OpenCV实现
（1）for循环遍历

#include <iostream>
#include <opencv2/opencv.hpp>
using namespace std;
using namespace cv;
int main()
{
	Mat img = imread("lena.png");
	imshow("lena",img);
	Size s = img.size();
	int c = img.channels();
	Mat img_gray = Mat::zeros(s, CV_32FC1);
	Mat img_sobel, img_sobelx, img_sobely;
	img_gray.copyTo(img_sobel);
	img_gray.copyTo(img_sobelx);
	img_gray.copyTo(img_sobely);
	// 灰度处理
	for (int i = 0; i < s.height; i++)
	{
		for (int j = 0; j < s.width; j++)
		{
			img_gray.at<float>(i, j) = (img.at<Vec3b>(i, j)[0] * 114 + img.at<Vec3b>(i, j)[1] * 578 + img.at<Vec3b>(i, j)[2] * 229) / 1000;
		}
	}
	for (int i = 1; i < s.height - 1; i++)
	{
		for (int j = 1; j < s.width - 1; j++)
		{
			img_sobelx.at<float>(i, j) = fabsf(img_gray.at<float>(i - 1, j + 1) - img_gray.at<float>(i - 1, j - 1) + 2 * img_gray.at<float>(i, j + 1) - 2 * img_gray.at<float>(i, j - 1) + img_gray.at<float>(i + 1, j + 1) - img_gray.at<float>(i + 1, j - 1));
			img_sobely.at<float>(i, j) = fabsf(img_gray.at<float>(i - 1, j - 1) - img_gray.at<float>(i + 1, j - 1) + 2 * img_gray.at<float>(i - 1, j) - 2 * img_gray.at<float>(i + 1, j) + img_gray.at<float>(i - 1, j + 1) - img_gray.at<float>(i + 1, j + 1));
		}
	}
	img_sobel = img_sobelx + img_sobely;
	for (int i = 1; i < s.height - 1; i++)
	{
		for (int j = 1; j < s.width - 1; j++)
		{
			if (img_sobel.at<float>(i, j) > 180)
				img_sobel.at<float>(i, j) = 1;
			else
				img_sobel.at<float>(i, j) = 0;
		}
	}
	imshow("img_sobel", img_sobel);
	waitKey(0);
}

结果如图所示：

（2）Sobel()函数

	// 求grad_x
	Sobel(img_gray, grad_x, img_gray.depth(), 0, 1);
	// 求grad_y
	Sobel(img_gray, grad_y, img_gray.depth(), 1, 0);
	addWeighted(grad_x, 0.5, grad_y, 0.5, 0, grad);
	for (int i = 0; i < s.height; i++)
	{
		for (int j = 0; j < s.width; j++)
		{
			if (grad.at<float>(i, j) > 80)
				grad.at<float>(i, j) = 1;
			else
				grad.at<float>(i, j) = 0;
		}
	}
	imshow("grad", grad);
	waitKey(0);

结果如图所示：