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=∣gradx∣∣grady∣weight=\frac{\lvert gradx \rvert}{\lvert grady \rvert}weight=∣grady∣∣gradx∣

图3、4的 weight=∣grady∣∣gradx∣weight=\frac{\lvert grady \rvert}{\lvert gradx \rvert}weight=∣gradx∣∣grady∣

梯度方向的判断

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

一般情况：
图1 ∣gradx∣<∣grady∣\lvert gradx \rvert \lt \lvert grady \rvert∣gradx∣<∣grady∣ 且两者同号
图2 ∣gradx∣<∣grady∣\lvert gradx \rvert \lt \lvert grady \rvert∣gradx∣<∣grady∣ 且两者异号
图3 ∣gradx∣>∣grady∣\lvert gradx \rvert \gt \lvert grady \rvert∣gradx∣>∣grady∣ 且两者同号
图4 ∣gradx∣>∣grady∣\lvert gradx \rvert \gt \lvert grady \rvert∣gradx∣>∣grady∣ 且两者异号

a、b梯度值的计算

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

图1中，
∣gradx∣<∣grady∣且两者同号grada=(1−weight)∗grad4+weight∗grad3gradb=(1−weight)∗grad8+weight∗grad7\lvert gradx \rvert \lt \lvert grady \rvert 且两者同号\\ grad_a=\left(1-weight\right)*grad_4+weight*grad_3 \\ grad_b=\left(1-weight\right)*grad_8+weight*grad_7 ∣gradx∣<∣grady∣且两者同号grada=(1−weight)∗grad4+weight∗grad3gradb=(1−weight)∗grad8+weight∗grad7
图2中，
∣gradx∣<∣grady∣且两者异号grada=(1−weight)∗grad4+weight∗grad5gradb=(1−weight)∗grad8+weight∗grad9\lvert gradx \rvert \lt \lvert grady \rvert 且两者异号\\ grad_a=\left(1-weight\right)*grad_4+weight*grad_5 \\ grad_b=\left(1-weight\right)*grad_8+weight*grad_9 ∣gradx∣<∣grady∣且两者异号grada=(1−weight)∗grad4+weight∗grad5gradb=(1−weight)∗grad8+weight∗grad9
图3中，
∣gradx∣>∣grady∣且两者同号grada=(1−weight)∗grad2+weight∗grad3gradb=(1−weight)∗grad6+weight∗grad7\lvert gradx \rvert \gt \lvert grady \rvert 且两者同号\\ grad_a=\left(1-weight\right)*grad_2+weight*grad_3\\ grad_b=\left(1-weight\right)*grad_6+weight*grad_7 ∣gradx∣>∣grady∣且两者同号grada=(1−weight)∗grad2+weight∗grad3gradb=(1−weight)∗grad6+weight∗grad7
图4中，
∣gradx∣>∣grady∣且两者异号grada=(1−weight)∗grad2+weight∗grad9gradb=(1−weight)∗grad6+weight∗grad5\lvert gradx \rvert \gt \lvert grady \rvert 且两者异号\\ grad_a=\left(1-weight\right)*grad_2+weight*grad_9\\ grad_b=\left(1-weight\right)*grad_6+weight*grad_5 ∣gradx∣>∣grady∣且两者异号grada=(1−weight)∗grad2+weight∗grad9gradb=(1−weight)∗grad6+weight∗grad5

代码

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;elsedire.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;}}}}
}