图像边缘Canny算子提取

图像边缘Canny算子提取

---

Blog: https://blog.youkuaiyun.com/AnimateX
Email: animatex.deng@gmail.com

本次项目中我觉得最有意思的部分就是梯度计算和边缘提取，如何提取更精细更准确的边缘则是最终的目的！

---

1 边缘检测

1.1 边缘定义

边缘检测是图像处理和计算机视觉中的基本问题，边缘检测的目的是标识数字图像中亮度变化明显的点。图像属性中的显著变化通常反映了属性的重要事件和变化。这些包括：
(1) 深度上的不连续；
(2) 表面方向不连续；
(3) 物质属性变化；
(4) 场景照明变化。

1.2 边缘检测方法分类

边缘检测是图像处理和计算机视觉中，尤其是特征检测中的一个研究领域。图像边缘检测大幅度地减少了数据量，并且剔除了可以认为不相关的信息，保留了图像重要的结构属性。有许多方法用于边缘检测，它们的绝大部分可以划分为两类：
(1) 基于查找一类；
(2) 基于零穿越的一类。
基于查找的方法通过寻找图像一阶导数中的最大和最小值来检测边界，通常是将边界定位在梯度最大的方向。基于零穿越的方法通过寻找图像二阶导数零穿越来寻找边界，通常是Laplacian过零点或者非线性差分表示的过零点。

1.3 影响边缘的因素

自然界图像的边缘并不总是理想的阶梯边缘。相反，它们通常受到一个或多个下面所列因素的影响：
(1) 有限场景深度带来的聚焦模糊.
(2) 非零半径光源产生的阴影带来的半影模糊.
(3) 光滑物体边缘的阴影.
(4) 物体边缘附近的局部镜面反射或者漫反射.

如果将边缘认为是一定数量点亮度发生变化的地方，那么边缘检测大体上就是计算这个亮度变化的导数。为简化起见，我们可以先在一维空间分析边缘检测。在这个例子中，我们的数据是一行不同点亮度的数据。例如，在下面的1维数据中我们可以直观地说在第4与第5个点之间有一个边界：

$$\hspace{2cm} 5 \hspace{1cm} 7 \hspace{1cm} 6 \hspace{1cm} 4 \hspace{1cm} 152 \hspace{1cm} 148 \hspace{1cm} 149 \hspace{2cm}$$
除非场景中的物体非常简单并且照明条件得到了很好的控制，否则确定一个用来判断两个相邻点之间有多大的亮度变化才算是有边界的阈值，并不是一件容易的事。实际上，这也是为什么边缘检测不是一个简单问题的原因之一。

2 边缘检测常见方法

有许多用于边缘检测的方法，他们大致可分为两类：基于搜索和基于零交叉。

基于搜索的边缘检测方法首先计算边缘强度，通常用一阶导数表示，例如梯度模;然后，用计算估计边缘的局部方向，通常采用梯度的方向，并利用此方向找到局部梯度模的最大值。

基于零交叉的方法找到由图像得到的二阶导数的零交叉点来定位边缘.通常用拉普拉斯算子或非线性微分方程的零交叉点，我们将在后面的小节中描述。

滤波做为边缘检测的预处理通常是必要的，通常采用高斯滤波。
已发表的边缘检测方法应用计算边界强度的度量，这与平滑滤波有本质的不同。正如许多边缘检测方法依赖于图像梯度的计算，他们用不同种类的滤波器来估计x-方向和y-方向的梯度。

2.1 一阶微分梯度算子

许多边缘检测操作都是基于亮度的一阶导数，方便得到原始数据的亮度梯度。一阶微分边缘算子也称为梯度边缘算子，它是利用图像在边缘处的阶跃性，即图像梯度在边缘取得极大值的特性进行边缘检测。梯度是一个矢量，它具有方向 $\theta$ 和模 $|\Delta I|$ ：

$$\Delta I = \begin{pmatrix} \frac{\partial I}{\partial x} \\ \frac{\partial I}{\partial y} \end{pmatrix}$$

$$|\Delta I| = \sqrt{(\frac{\partial I}{\partial x} )^2+(\frac{\partial I}{\partial y} )^2} = \sqrt{I_x^2+I_y^2}$$

$$\theta = arctan(I_y / I_x)$$

梯度的方向提供了边缘的趋势信息，因为梯度方向始终是垂直于边缘方向，梯度的模值大小提供了边缘的强度信息。

在实际使用中，通常利用有限差分进行梯度近似。对于上面的公式，我们有如下的近似：

$$\frac{\partial I}{\partial x} = \lim_{h\to 0}\frac{I(x+\Delta x,y)-I(x,y)}{\Delta x} \approx I(x+1,y)-I(x,y),(\Delta x = 1)$$
$$\frac{\partial I}{\partial y} = \lim_{h\to 0}\frac{I(x,y+\Delta xy)-I(x,y)}{\Delta y} \approx I(x,y+1)-I(x,y),(\Delta y = 1)$$

2.2 Robert梯度算子

1963年，Roberts提出了这种寻找边缘的算子。Roberts边缘算子是一个2x2的模板，采用的是对角方向相邻的两个像素之差。从图像处理的实际效果来看，边缘定位较准，对噪声敏感。在Roberts检测算子中：

$$\frac{\partial I}{\partial x} = I(i,j)-I(i+1,j+1)$$
$$\frac{\partial I}{\partial y} = I(i+1,j)-I(i,j+1)$$
则可以得到对应的x方向和y方向上的卷积核：
$$m_x = \begin{bmatrix}1&0 \\ 0&-1\end{bmatrix}, m_y=\begin{bmatrix}0&-1 \\ 1&0\end{bmatrix}$$

2.3 Prewitt梯度算子

利用周边八邻域的灰度值来估计中心梯度。其简化梯度计算公式：

$$\frac{\partial I}{\partial x} = I(i-1,j+1)+I(i,j+1)+I(i+1,j+1)-I(i-1,j-1)-I(i,j-1)-I(i+1,j-1)$$
$$\frac{\partial I}{\partial y} =-I(i+1,j-1)+I(i+1,j)+I(i+1,j+1)-I(i-1,j-1)+I(i-1,j)+I(i-1,j+1)$$
对应的卷积核为：
$$m_x = \begin{bmatrix}-1&0&+1 \\ -1&0&+1 \\ –1&0&+1\end{bmatrix}, m_y=\begin{bmatrix}-1&-1&-1 \\ 0&0&0 \\ +1&+1&+1\end{bmatrix}$$

2.4 Sobel梯度算子

相比于Prewitt梯度算子，Sobel算子提高了靠近中间位置的权重。其对应的卷积核为：

$$m_x = \begin{bmatrix}-1&0&+1 \\ -2&0&+2\\ –1&0&+1\end{bmatrix}, m_y=\begin{bmatrix}-1&-2&-1 \\ 0&0&0 \\ +1&+2&+1\end{bmatrix}$$

3 二阶微分梯度算子

由一阶微分可知，边缘即是图像的一阶导数局部最大值的地方，那么也意味着该点的二阶导数为零。二阶微分边缘检测算子就是利用图像在边缘处的阶跃性导致图像二阶微分在边缘处出现零值这一特性进行边缘检测的。
对于图像的二阶微分可以用拉普拉斯算子来表示：

$$\nabla^2I = \frac{\partial^2I}{\partial x^2}+\frac{\partial^2I}{\partial y^2}$$
我们在像素点(i,j)的3×3的邻域内，可以有如下的近似：
$$\frac{\partial^2I}{\partial x^2} = I(i,j+1)-2I(i,j)+I(i,j-1)$$
$$\frac{\partial^2I}{\partial y^2}=I(i+1,j)-2I(i,j)+I(i-1,j)$$
$$\nabla^2I = –4I(i,j)+I(i,j+1)+I(i,j-1)+I(i+1,j)+I(i-1,j)$$
对应的二阶微分卷积核为：
$$\boldsymbol{m} = \begin{bmatrix}0&1&0 \\ 1&4&1\\ 0&1&0\end{bmatrix}$$
所以二阶微分检测边缘的方法就分两步：1）用上面的Laplace核与图像进行卷积；2）对卷积后的图像，去掉那些卷积结果为0的点。

虽然上述使用二阶微分检测边缘的方法简单，但它的缺点是对噪声十分敏感，同时也没有能够提供边缘的方向信息。为了实现对噪声的抑制，Marr等提出了LOG的方法。

为了减少噪声对边缘的影响，首先图像要进行低通滤波，LOG采用了高斯函数作为低通滤波器。高斯函数为：

$$G(x,y)=\frac{1}{2\pi\sigma^2}e^{-\frac{x^2+y^2}{2\sigma^2}}$$
上面的公式中σ决定了对图像的平滑程度。高斯函数生成的滤波模板尺寸一般设定为 $6\sigma+1$（加1是会了使滤波器的尺寸为奇数）。使用高斯函数对图像进行滤波并对图像滤波结果进行二阶微分运算的过程，可以转换为先对高斯函数进行二阶微分，再利用高斯函数的二阶微分结果对图像进行卷积运算：
$$\nabla^2[G(x,y)*f(x,y)] = \nabla^2[G(x,y)]*f(x,y)$$
$$\nabla^2G(x,y) = -\frac{1}{2\pi\sigma^4}[1-\frac{x^2+y^2}{\sigma^2}]\cdot exp(-\frac{x^2+y^2}{2\sigma^2})$$

4 Canny 边缘检测

canny边缘检测实际上是一种一阶微分算子检测算法，但为什么这里拿出来说呢，因为它几乎是边缘检测算子中最为常用的一种，也是个人认为现在最优秀的边缘检测算子。Canny提出了边缘检测算子优劣评判的三条标准：

高的检测率。边缘检测算子应该只对边缘进行响应，检测算子不漏检任何边缘，也不应该将非边缘标记为边缘。
精确定位。检测到的边缘与实际边缘之间的距离要尽可能的小。
明确的响应。对每一条边缘只有一次响应，只得到一个点。
Canny边缘检测之所以优秀是因为它在一阶微分算子的基础上，增加了非最大值抑制和双阈值两项改进。利用非极大值抑制不仅可以有效地抑制多响应边缘，而且还可以提高边缘的定位精度；利用双阈值可以有效减少边缘的漏检率。

Canny边缘检测主要分四步进行：
1 去噪声；
2 计算梯度与方向角；
3 非最大值抑制；
4 滞后阈值化；
其中前两步很简单，先用一个高斯滤波器对图像进行滤波，然后用Sobel水平和竖直检测子与图像卷积，来计算梯度和方向角。

4.1 去噪

一般用高斯滤波。高斯滤波是一种线性平滑滤波，适用于消除高斯噪声，广泛应用于图像处理的减噪过程。其本质是一个低通滤波器，查看其模板可以看出。
注：理论上，高斯分布在所有定义域上都有非负值，这就需要一个无限大的卷积核。在计算高斯函数的离散近似时，在大概3σ距离之外的像素都可以看作不起作用，这些像素的计算也就可以忽略。通常，图像处理程序只需要计算 ${\displaystyle (6\sigma +1)\times (6\sigma +1)} (6 \sigma + 1) \times (6 \sigma + 1)$的矩阵就可以保证相关像素影响。对于边界上的点，通常采用复制周围的点到另一面再进行加权平均运算。
matlab实现举例：

w   = fspecial('gausian', [3 3], 0.5);
img = imfilter(I, w, 'replicate');

4.2 计算梯度与方向角

利用前面讲的Sobel算子计算。
MATLAB实现：

wx_sobel = [1 0 -1; 2 0 -2; 1 0 -1];
wy_sobel = [1 2 1; 0 0 0; -1 -2 -1];
edgeX = imfilter(tempEdge, wx_sobel, 'replicate');
edgeY = imfilter(tempEdge, wy_sobel, 'replicate');
imgEdge(:, :, i) = sqrt(edgeX .^ 2 + edgeY .^ 2);
% for speed
% imgEdge(:, :, i) = abs(edgeX) + abs(edgeY);

thetaTemp( edgeX ~= 0 ) = atan( edgeY ./ edgeX );

4.3 非极大值抑制

图像梯度幅值矩阵中的元素值越大，说明图像中该点的梯度值越大，但这不不能说明该点就是边缘（也许两边的图像刚好在增加或者减小，表现出来增长的过程）。在Canny算法中，非极大值抑制是进行边缘检测的重要步骤，通俗意义上是指寻找像素点（边缘）局部最大值，将非极大值点所对应的灰度值置为0，这样可以剔除掉一大部分非边缘的点。

先对上图作一个说明，每一个格子代表一个像素点，P(i,j)为坐标系原点，绿色的线为我的理想边缘，我们要检查A点是不是我们的边缘，那么需要看其是不是梯度方向上的最大值。蓝色线条为梯度方向（垂直于边缘方向），要确认A点是不是梯度方向的极大值，简单办法就是看A点是不是比梯度方向左右两边的点都大。

非极大值抑制原理说明：已知需要判断的点A的值和对应的梯度方向，找到A点在正负梯度方向两个邻近点B、C的值，比较A点的值是不是比B、C都大，如果不满足该条件就将A点的值置为最小值（一般是0，对于红外图像需要处理）。

问题在于B C两点不一定在0°、45°、90°、135°几个方向上，那怎么比较？

除了在梯度方向刚好在0°、45°、90°、135°几个方向时，B、C在特殊点。其他点无法直接得到，那么只能通过插值得到。如上图所示，我们要要得到B、C两点的值也就是M1和M2，我们已知M1、M2两点最邻近两点的值和x和y方向的梯度值，那么权重可以用Gx/Gy表示，用线性插值公式：
$$weight = \frac{|G_{x,y}|}{|G_{y,x}|}$$
$$M = P_0\times weight + (1-weight)\times P_1$$
则我们可以得到如上图所示的计算公式，代码如下：

% 代码计算梯度角预先加上pi/2, 因为atan计算出来的范围是[-pi/2, pi/2]
% 0 ~ pi/4
M1 = p(i-1, j+1) * weight + (1 - weight) * p(i, j+1);
M2 = p(i+1, j-1) * weight + (1 - weight) * p(i, j-1);
% pi/4 ~ pi/2
M1 = p(i-1, j+1) * weight + (1 - weight) * p(i-1, j);
M2 = p(i+1, j-1) * weight + (1 - weight) * p(i+1, j);

同理，下图也能用同样的方法进行计算：

90°到135°，135°到180°，M1和M2值计算的代码如下:

% 代码计算梯度角预先加上pi/2, 因为atan计算出来的范围是[-pi/2, pi/2]
% pi/2 ~ 3pi/2
M1 = p(i-1, j-1) * weight + (1 - weight) * p(i-1, j);
M2 = p(i+1, j+1) * weight + (1 - weight) * p(i+1, j);
% 3pi/2 ~ pi
M1 = p(i-1, j-1) * weight + (1 - weight) * p(i, j-1);
M2 = p(i+1, j+1) * weight + (1 - weight) * p(i, j+1);

4.4 滞后阈值化

由于噪声的影响，经常会在本应该连续的边缘出现断裂的问题。滞后阈值化设定两个阈值：一个为高阈值Th，一个为低阈值Tl。如果任何像素边缘算子的影响超过高阈值，将这些像素标记为边缘。所以不整个过程描述如下：

如果该像素的梯度值小于Tl，则该像素为非边缘像素；
如果该像素的梯度值大于Th，则该像素为边缘像素；

但是阈值怎么设置呢？我们先来看 MATLAB 的edge函数，其中有这么一段

function [lowThresh, highThresh] = selectThresholds(thresh, magGrad, PercentOfPixelsNotEdges, ThresholdRatio, ~)

[m,n] = size(magGrad);

% Select the thresholds
if isempty(thresh)
    % 计算统计直方图，灰阶设置为6阶
    counts=imhist(magGrad, 64);
    % 设置累积直方图的总数大于0.7总像素点数的归一化值为高阈值
    highThresh = find(cumsum(counts) > PercentOfPixelsNotEdges*m*n,...
        1,'first') / 64;
    % 低阈值为0.4*highThresh
    lowThresh = ThresholdRatio*highThresh;
elseif length(thresh)==1
    highThresh = thresh;
    if thresh>=1
        error(message('images:edge:thresholdMustBeLessThanOne'))
    end
    lowThresh = ThresholdRatio*thresh;
elseif length(thresh)==2
    lowThresh = thresh(1);
    highThresh = thresh(2);
    if (lowThresh >= highThresh) || (highThresh >= 1)
        error(message('images:edge:thresholdOutOfRange'))
    end
end

可见光图像我没有测试，我对14位的原始数据测试，直方图太奇特，无法用此种方法实现。变为固定阈值。

4.5 弱边缘检测（连通性分析）

再完成双阈值处理后，响应超过低阈值（高低阈值之间）的像素，如果与已经标记为边缘的像素4-邻接或8-邻接，则将这些像素也标记为边缘。
如果该像素的梯度值介于Tl与Th之间，需要进一步检测该像素的3×3邻域内的8个点，如果这8个点内有一个或以上的点梯度超过了Th，则该像素为边缘像素，否则不是边缘像素。
这种方法一样不适合红外灰度数据分析。

5 总结

Canny最重要的是非极大值抑制和阈值处理以及弱边缘检测。
可是如果在FPGA中做呢？如何做？

1. 关于非极大值抑制
   
   如上图所示，我们将区域变为①~⑤，直接用临近点代替M1和M2的值。
2. 阈值处理
   考虑从均值入手（考虑红外的背景问题）。

CODE

function imgOut = edgeDetecNew(imgIn, rad, Sigma)
    [M, N, num] = size(imgIn);
    imgOut = zeros(M, N, num);
    %
    % Version 2.0
    % Date:2018-04-18
    % Author: Deng
    % This is test algorithm
    % Use Canny operator
    %--------------------------------------------------------------%
    %---------------------- algorithm steps -----------------------%
    % 1. Gaussian smoothing of the original image
    % 2. Perform sobel edge detection on Gaussian-smoothed images.
    %    There is also a need for horizontal and vertical joints.
    % 3. Non-maximally suppression of joint sobel detection images
    % 4. Connect edge points and perform hysteresis thresholding
    %--------------------------------------------------------------%

    w_gauss = fspecial('gaussian', [rad rad], Sigma);
    gaussTemp = zeros(M, N, num);

    % sobel operator
    wx_sobel = [1 0 -1; 2 0 -2; 1 0 -1];
    wy_sobel = [1 2 1; 0 0 0; -1 -2 -1];

    Theta    = zeros(M, N, num);
    imgEdge  = zeros(M, N, num);
    for i = 1 : num
        temp = imgIn(:, :, i);
        gaussTemp(:, :, i) = imfilter(temp, w_gauss, 'replicate');
        figure, imshow(gaussTemp, [ ]);
        % Theta is Theta is the angle between edgeY and edgeX.
        thetaTemp = Theta(:, :, i);
        tempEdge  = gaussTemp(:, :, i);
        edgeX = imfilter(tempEdge, wx_sobel, 'replicate');
        edgeY = imfilter(tempEdge, wy_sobel, 'replicate');
        imgEdge(:, :, i) = sqrt(edgeX .^ 2 + edgeY .^ 2);
        % for speed
        % imgEdge(:, :, i) = abs(edgeX) + abs(edgeY);

        thetaTemp( edgeX ~= 0 ) = atan( edgeY ./ edgeX );

        other = thetaTemp( edgeX == 0 );
        thetaTemp( edgeY( other ) > 0 ) = pi / 2;
        thetaTemp( edgeY( other ) < 0 ) = - pi /2;
        Theta(:, :, i) = thetaTemp;
    end


    % The key two steps
    % Non-maximally suppression
    for i = 1 : num
        thetaTemp = Theta(:, :, i);
        edgeTemp  = imgEdge(:, :, i);
        Edge  = padarray(edgeTemp,  [1 1], 'replicate');
        theta = padarray(thetaTemp, [1 1], 'replicate');

        %------------------test Ir flag------------------%
        % flag  = mean( Edge(:) );
        flag = 0;
        %------------------------------------------------%

        for j = 2 : M+1
            for k = 2 : N+1
                interP = Edge(j-1 : j+1, k-1 : k+1);
                % The first interpolation point: M1
                % The second interpolation point: M2
                if theta(j, k) == -pi/2
                    M1 = interP(3, 2);
                    M2 = interP(1, 2);
                    if interP(2, 2) < M1 || interP(2, 2) < M2
                        Edge(j, k) = flag;
                    end
                elseif theta(j, k) > -pi/2 && theta(j, k) < -pi /4
                    weight = abs( edgeX(j-1, k-1) ) / abs( edgeY(j-1, k-1) );
                    M1 = weight * interP(3, 3) + (1 - weight) * interP(3, 2);
                    M2 = weight * interP(1, 1) + (1 - weight) * interP(1, 2);
                    if interP(2, 2) < M1 || interP(2, 2) < M2
                        Edge(j, k) = flag;
                    end
                elseif theta(j, k) >= -pi /4 && theta(j, k) < 0
                    weight = abs( edgeY(j-1, k-1) ) / abs( edgeX(j-1, k-1) );
                    M1 = weight * interP(3, 3) + (1 - weight) * interP(2, 3);
                    M2 = weight * interP(1, 1) + (1 - weight) * interP(2, 1);
                    if interP(2, 2) < M1 || interP(2, 2) < M2
                        Edge(j, k) = flag;
                    end
                elseif theta(j, k) == 0
                    M1 = interP(2, 3);
                    M2 = interP(2, 1);
                    if interP(2, 2) < M1 || interP(2, 2) < M2
                        Edge(j, k) = flag;
                    end
                elseif theta(j, k) > 0 && theta(j, k) < pi / 4
                    weight = abs( edgeY(j-1, k-1) ) / abs( edgeX(j-1, k-1) );
                    M1 = weight * interP(1, 3) + (1 - weight) * interP(2, 3);
                    M2 = weight * interP(3, 1) + (1 - weight) * interP(2, 1);
                    if interP(2, 2) < M1 || interP(2, 2) < M2
                        Edge(j, k) = flag;
                    end
                elseif theta(j, k) >= pi/4 && theta(j, k) < pi/2
                    weight = abs( edgeX(j-1, k-1) ) / abs( edgeY(j-1, k-1) );
                    M1 = weight * interP(1, 3) + (1 - weight) * interP(1, 2);
                    M2 = weight * interP(3, 1) + (1 - weight) * interP(3, 2);
                    if interP(2, 2) < M1 || interP(2, 2) < M2
                        Edge(j, k) = flag;
                    end
                else
                    M1 = interP(1, 2);
                    M2 = interP(3, 2);
                    if interP(2, 2) < M1 || interP(2, 2) < M2
                        Edge(j, k) = flag;
                    end
                end
            end
        end
        figure, imshow(Edge, [ ]);

        % Double threshold processing
        % The first steps: find th by use hist.

        %------------------------- Test Code -------------------------%
        HistTemp = Edge(2 : M+1, 2 : N+1);
        meanVaule = mean( HistTemp(:) );
        maxVaule = max( HistTemp(:) );

        % XX需要根据实际情况调整
        th = XX * meanVaule;
        tl = XX * th;

        HistTemp( HistTemp <=  tl ) = flag;
        HistTemp( HistTemp >= th ) = maxVaule;
        HistTemp = imfilter(HistTemp, w_gauss, 'replicate');
        figure, imshow(HistTemp, [ ]);

        for j = 1 : M-2
            for k = 1 : N-2
                connectTemp = HistTemp(j : j+2, k : k+2);
                if HistTemp(j, k) > tl && HistTemp(j, k) < th
                    if numel( connectTemp >= maxVaule ) > 1
                        HistTemp(j, k) =  maxVaule;
                    else
                        HistTemp(j, k) = flag;
                    end
                end
            end
        end
        %--------------------------------------------------------------%
        imgOut(:, :, i) = HistTemp;
    end
end

Reference

[1] https://zh.wikipedia.org/wiki/%E9%AB%98%E6%96%AF%E6%A8%A1%E7%B3%8A
[2] https://zh.wikipedia.org/wiki/Canny%E7%AE%97%E5%AD%90
[3] https://zh.wikipedia.org/wiki/%E8%BE%B9%E7%BC%98%E6%A3%80%E6%B5%8B
[4] https://blog.youkuaiyun.com/likezhaobin/article/details/6892176
[5] http://www.cnblogs.com/tiandsp/archive/2012/12/13/2817240.html