Canny边缘检测Step3——梯度幅值的非极大值抑制（NMS）详解

一、背景知识

  经过Canny边缘检测的前两步：

• Canny边缘检测Step1——二维高斯滤波
• Canny边缘检测Step2——Sobel算子计算梯度信息

  我们到了第三步——梯度幅值的非极大值抑制（NMS）
  Canny算子中的非极大值抑制是在对图像进行梯度求取之后，在梯度方向进行的运算。在John Canny提出Canny算子的论文中，非最大值抑制只是在0、90、45、135四个梯度方向上进行，每个像素点梯度方向按照相近程度用这四个方向来代替。
  这样做的好处是简单， 但是这种简化的方法无法达到最好的抑制效果。 因为自然图像中的边缘梯度方向不一定只可能沿着这四个方向。实际图像中的像素点是离散的二维矩阵，处于中心位置C处的梯度方向的两侧的点不一定存在，为了方便表示，我们可以称这两个点为亚像素点，而这两个亚像素点的梯度值就必须通过使用其两侧像素点的梯度值进行插值计算得到。
  设C为中心像素点，其梯度值为grad(i,j)，水平方向梯度为gradx(i,j)，竖直方向梯度值为grady(i,j)，设插值使用的权重为weight，设亚像素点梯度值分别为gradTemp1(i,j)和gradTemp2(i,j)，设插值用到的四个C的邻域点的梯度值分别为grad1(i,j)，grad2(i,j)，grad3(i,j)，grad4(i,j)。
  过C点的梯度方向线共有下述4种情况：
  ① |grady|>|gradx|且gradx与grady同号；
  ② |grady|>|gradx|且gradx与grady异号；
  ③ |grady|<|gradx|且gradx与grady同号；
  ④ |grady|<|gradx|且gradx与grady异号。
  下面这幅图可以辅助我们理解文字描述：

  NMS算法步骤如下（看起来复杂，但编程很简单）：
  若|grady|>|gradx|，则weight(i,j) = |gradx|/|grady|，grad2(i,j) = grad(i-1,j)，grad4(i,j) = grad(i+1,j)。在这种条件下，若gradx与grady同号，则grad1(i,j) = grad(i-1,j-1)，grad3(i,j) = grad(i+1,j+1)，否则grad1(i,j) = grad(i-1,j+1)，grad3(i,j) = grad(i+1,j-1)。
  若|grady|≤|gradx|，则weight(i,j) = |grady|/|gradx|，grad2(i,j) = grad(i,j-1)，grad4(i,j) = grad(i,j+1)。在这种条件下，若gradx与grady同号，则grad1(i,j) = grad(i+1,j-1)，grad3(i,j) = grad(i-1,j+1)，否则grad1(i,j) = grad(i-1,j-1)，grad3(i,j) = grad(i+1,j+1)。
  得到每种情况下的weight(i,j)值以及四个用于插值计算的点的梯度值grad1(i,j)，grad2(i,j)，grad3(i,j)，grad4(i,j)后，便可以计算出亚像素点的梯度值：

  得到两个亚像素点的梯度值后，将中心像素点梯度值与这两个亚像素点梯度值进行比较，若中心点梯度值是这三个点中的最大值，则在结果矩阵中保留中心点梯度值，否则在结果中置0。

二、Matlab实现代码（函数形式）

% 非极大值抑制函数
function output = NMS_algo(grad,gradx,grady)
% 非极大值抑制
grad = double(grad);
[h,w] = size(grad);
% 初始化结果
result = zeros(h,w);
% 抑制算法实现
% 用线性插值算出梯度方向上亚像素点的梯度值
% 权重初始化
weight = zeros(h,w);
% 各像素点梯度值初始化
grad1 = zeros(h,w);
grad2 = zeros(h,w);
grad3 = zeros(h,w);
grad4 = zeros(h,w);
gradTemp1 = zeros(h,w);
gradTemp2 = zeros(h,w);
for i = 2:(h-1)
    for j = 2:(w-1)
        % 梯度为0的点直接在结果中置0
        if grad(i,j) == 0
            result(i,j) = 0;
        % 对于梯度不为0的点
        else
            % y方向梯度大于x方向梯度
            if abs(grady(i,j)) > abs(gradx(i,j))
                % 计算权重
                weight(i,j) = abs(gradx(i,j))/abs(grady(i,j));
                % 中心点上下两点
                grad2(i,j) = grad(i-1,j);
                grad4(i,j) = grad(i+1,j);
                % gradx和grady同号
                if gradx(i,j) * grady(i,j) > 0
                    % 插值用到的另外两点
                    grad1(i,j) = grad(i-1,j-1);
                    grad3(i,j) = grad(i+1,j+1);
                % gradx和grady异号
                else
                    grad1(i,j) = grad(i-1,j+1);
                    grad3(i,j) = grad(i+1,j-1);
                end
            % y方向梯度小于x方向梯度
            else
                weight(i,j) = abs(grady(i,j))/abs(gradx(i,j));
                grad2(i,j) = grad(i,j-1);
                grad4(i,j) = grad(i,j+1);
                % gradx和grady同号
                if gradx(i,j) * grady(i,j) > 0
                    grad1(i,j) = grad(i+1,j-1);
                    grad3(i,j) = grad(i-1,j+1);
                % gradx和grady异号    
                else
                    grad1(i,j) = grad(i-1,j-1);
                    grad3(i,j) = grad(i+1,j+1);
                end
            end
        end
        % 插值计算计算两个亚像素点梯度值
        gradTemp1(i,j) = weight(i,j) * grad1(i,j) + (1 - weight(i,j)) * grad2(i,j);
        gradTemp2(i,j) = weight(i,j) * grad3(i,j) + (1 - weight(i,j)) * grad4(i,j);
        % 比较中心像素点和两个亚像素点的梯度值
        if (grad(i,j) >= gradTemp1(i,j) && grad(i,j) >= gradTemp2(i,j))
            % 中心点在其邻域内为极大值，在结果中保留其梯度值
            result(i,j) = grad(i,j);
        else
            % 否则的话，在结果中置0
            result(i,j) = 0;
        end     
    end
end
output = uint8(result);
end

三、实验结果

  非极大值抑制后的图像如下图所示：