Canny算子边缘检测——非极大值抑制Non-Maximum Suppression

Canny算子边缘检测——非极大值抑制Non-Maximum Suppression

Canny算子计算前需要先通过Sobel算子得到图像的梯度值图和方向图。非极大值抑制指的是在一个像素点上A，将它的梯度值和它梯度方向上的前后两个点B和C的进行比较，如果A的梯度值大于B和C的梯度值，则保留其值，否则将A的值置为0。

简单情况

像素点的梯度方向在0°，45°，90°和135°上，此时像素点梯度方向上的前后两点就是它上、下、左、右、左上、左下、右上和右下八个点的其中两个，直接比较大小即可。例如：如果Grad(1)>=Grad(3) && Grad(1)>=Grad(7)，那么就保留1的梯度值，否则将其梯度值置为0。

一般情况

梯度方向不是0°，45°，90°或135°，则1点前后的点不能落在像素点上，因此需要添加如图的a、b两点(称为子像素或亚像素)，通过计算一定的权重和周边点的梯度值来计算这两个点的梯度值，从而与1进行比较。

具体分类及计算过程如下：

分成4类：

图1

图2

图3

图4

梯度方向向量经过了哪四个点的连线，我们就将用哪四个点来计算，例如图1的向量经过了点3和4、点7和8的连线。

权重weight

weight定义为一点处min(abs(水平梯度)，abs(竖直梯度))/max(abs(水平梯度)，abs(竖直梯度))。

图1、2的 $weight=\frac{|gradx|}{|grady|}$

图3、4的 $weight=\frac{|grady|}{|gradx|}$

梯度方向的判断

0°方向时grady为0
45°方向时梯度值为1
90°方向时gradx为0
135°方向时梯度值为-1

一般情况：
图1 $|gradx|<|grady|$ 且两者同号
图2 $|gradx|<|grady|$ 且两者异号
图3 $|gradx|>|grady|$ 且两者同号
图4 $|gradx|>|grady|$ 且两者异号

a、b梯度值的计算

该点水平(或竖直)线上的点的的梯度值乘上1-weight加上与之连线的点的梯度值乘上weight即可。

图1中，
$$\begin{gathered} |gradx|<|grady|且两者同号 \\ gra{d}_a=\left(1-weight\right)\ast gra{d}_4+weight\ast gra{d}_3 \\ gra{d}_b=\left(1-weight\right)\ast gra{d}_8+weight\ast gra{d}_7 \end{gathered}$$
图2中，
$$\begin{gathered} |gradx|<|grady|且两者异号 \\ gra{d}_a=\left(1-weight\right)\ast gra{d}_4+weight\ast gra{d}_5 \\ gra{d}_b=\left(1-weight\right)\ast gra{d}_8+weight\ast gra{d}_9 \end{gathered}$$
图3中，
$$\begin{gathered} |gradx|>|grady|且两者同号 \\ gra{d}_a=\left(1-weight\right)\ast gra{d}_2+weight\ast gra{d}_3 \\ gra{d}_b=\left(1-weight\right)\ast gra{d}_6+weight\ast gra{d}_7 \end{gathered}$$
图4中，
$$\begin{gathered} |gradx|>|grady|且两者异号 \\ gra{d}_a=\left(1-weight\right)\ast gra{d}_2+weight\ast gra{d}_9 \\ gra{d}_b=\left(1-weight\right)\ast gra{d}_6+weight\ast gra{d}_5 \end{gathered}$$

代码

void SobelGradientDirection(const Mat& img, Mat& dire)
{
/*
img是灰度图，dire是每一个像素点的梯度方向
sobelEdgeX和sobelEdgY是通过Sobel算子得到X和Y方向上的梯度值(自己写的函数)
*/
	Mat EdgeX, EdgeY;
	sobelEdgeX(img, EdgeX);
	sobelEdgeY(img, EdgeY);
	int row = EdgeX.size().height;
	int col = EdgeX.size().width;
	dire = Mat(row, col, CV_64FC1);
	for (int i = 0; i < row; i++)
		for (int j = 0; j < col; j++)
			//由于atan的值为-pi/2~pi/2，所以给加90使角度均为正值 范围在0~180之间
			if (int(EdgeX.ptr<uchar>(i)[j] == 0))
					dire.ptr<double>(i)[j] = 90;
			else
				dire.ptr<double>(i)[j] = int(EdgeY.ptr<double>(i)[j]) / int(EdgeX.ptr<uchar>(i)[j]) * 180 / 3.14159265 + 90;
}

/*
Grad是梯度大小图像  
GradX是X梯度幅值图像  
GradY是Y梯度幅值图像  
dire是梯度方向图像
refined是非极大值抑制后的图像
*/
void NonMaxSup(const Mat& Grad, const Mat& GradX, const Mat& GradY, const Mat& dire, Mat& refined)
{
	int row = Grad.size().height;
	int col = Grad.size().width;
	refined = Mat(row, col, CV_8UC1);
	for (int i = 0; i < row; i++)
	{
		for (int j = 0; j < col; j++)
		{
			uchar grad, grad1, grad2, grad3, grad4, grada, gradb;
			double weight;
			if (dire.ptr<double>(i)[j] == 0)
			{
				grad = Grad.ptr<uchar>(i)[j];
				grad1 = Grad.ptr<uchar>(i)[j - 1];
				grad2 = Grad.ptr<uchar>(i)[j + 1];
				refined.ptr<uchar>(i)[j] = max(grad, max(grad1, grad2)) == grad ? grad : 0;
			}
			else if (dire.ptr<double>(i)[j] == 45)
			{
				grad = Grad.ptr<uchar>(i)[j];
				grad1 = Grad.ptr<uchar>(i - 1)[j + 1];
				grad2 = Grad.ptr<uchar>(i + 1)[j - 1];
				refined.ptr<uchar>(i)[j] = max(grad, max(grad1, grad2)) == grad ? grad : 0;
			}
			else if (dire.ptr<double>(i)[j] == 90)
			{
				grad = Grad.ptr<uchar>(i)[j];
				grad1 = Grad.ptr<uchar>(i - 1)[j];
				grad2 = Grad.ptr<uchar>(i + 1)[j];
				refined.ptr<uchar>(i)[j] = max(grad, max(grad1, grad2)) == grad ? grad : 0;
			}
			else if (dire.ptr<double>(i)[j] == 135)
			{
				grad = Grad.ptr<uchar>(i)[j];
				grad1 = Grad.ptr<uchar>(i - 1)[j - 1];
				grad2 = Grad.ptr<uchar>(i + 1)[j + 1];
				refined.ptr<uchar>(i)[j] = max(grad, max(grad1, grad2)) == grad ? grad : 0;
			}
			else if (abs(int(GradX.ptr<uchar>(i)[j])) < abs(int(GradY.ptr<double>(i)[j])))
				if (int(GradX.ptr<uchar>(i)[j]) * int(GradY.ptr<uchar>(i)[j]) < 0)
				{
					grad = Grad.ptr<uchar>(i)[j];
					grad1 = Grad.ptr<uchar>(i - 1)[j];
					grad2 = Grad.ptr<uchar>(i + 1)[j];
					grad3 = Grad.ptr<uchar>(i - 1)[j + 1];
					grad4 = Grad.ptr<uchar>(i + 1)[j - 1];
					weight = double(abs(int(GradX.ptr<uchar>(i)[j]))) / double(abs(int(GradY.ptr<double>(i)[j])));
					grada = weight * grad3 + (1 - weight) * grad1;
					gradb = weight * grad4 + (1 - weight) * grad2;
					refined.ptr<uchar>(i)[j] = max(grad, max(grada, gradb)) == grad ? grad : 0;
				}
				else
				{
					grad = Grad.ptr<uchar>(i)[j];
					grad1 = Grad.ptr<uchar>(i - 1)[j];
					grad2 = Grad.ptr<uchar>(i + 1)[j];
					grad3 = Grad.ptr<uchar>(i - 1)[j - 1];
					grad4 = Grad.ptr<uchar>(i + 1)[j + 1];
					weight = double(abs(int(GradX.ptr<uchar>(i)[j]))) / double(abs(int(GradY.ptr<double>(i)[j])));
					grada = weight * grad3 + (1 - weight) * grad1;
					gradb = weight * grad4 + (1 - weight) * grad2;
					refined.ptr<uchar>(i)[j] = max(grad, max(grada, gradb)) == grad ? grad : 0;
				}
			else
			{
				if (int(GradX.ptr<uchar>(i)[j]) * int(GradY.ptr<uchar>(i)[j]) < 0)
				{
					grad = Grad.ptr<uchar>(i)[j];
					grad1 = Grad.ptr<uchar>(i)[j + 1];
					grad2 = Grad.ptr<uchar>(i)[j - 1];
					grad3 = Grad.ptr<uchar>(i - 1)[j + 1];
					grad4 = Grad.ptr<uchar>(i + 1)[j - 1];
					weight = double(abs(int(GradY.ptr<uchar>(i)[j]))) / double(abs(int(GradX.ptr<double>(i)[j])));
					grada = weight * grad3 + (1 - weight) * grad1;
					gradb = weight * grad4 + (1 - weight) * grad2;
					refined.ptr<uchar>(i)[j] = max(grad, max(grada, gradb)) == grad ? grad : 0;
				}
				else
				{
					grad = Grad.ptr<uchar>(i)[j];
					grad1 = Grad.ptr<uchar>(i)[j - 1];
					grad2 = Grad.ptr<uchar>(i)[j + 1];
					grad3 = Grad.ptr<uchar>(i - 1)[j - 1];
					grad4 = Grad.ptr<uchar>(i + 1)[j + 1];
					weight = double(abs(int(GradY.ptr<uchar>(i)[j]))) / double(abs(int(GradX.ptr<double>(i)[j])));
					grada = weight * grad3 + (1 - weight) * grad1;
					gradb = weight * grad4 + (1 - weight) * grad2;
					refined.ptr<uchar>(i)[j] = max(grad, max(grada, gradb)) == grad ? grad : 0;
				}
			}
		}
	}
}