11.图像边缘检测的原理与实现

数字图像处理(19): 边缘检测算子(Roberts算子、Prewitt算子、Sobel算子 和 Laplacian算子)
数字图像处理(20): 边缘检测算子(Canny算子)

1.边缘检测介绍

1.1 边缘检测的基本原理

  边缘是图像的基本特征，所谓的边缘就是指的图像的局部不连续性。灰度或者结构等信息的突变之处称之为边缘。如灰度级的突变、颜色的突变、纹理结构的突变等。边缘是一个区域的结束，也是另一个区域的开始，利用该特征可以分割图像。
  图像的边缘有方向和幅度两种特性，边缘通常可以通过一阶导数或二阶导数检测得到。一阶导数是以最大值作为对应边缘的位置，而二阶导数则以过零点作为对应边缘的位置。
  边缘检测是一种常用的图像分割技术，常用的边缘检测算子有Roberts Cross算子、Prewitt算子、Sobel算子、Kirsch算子、Laplacian算子以及Canny算子。

1.2 边缘检测算子分类

（1）一阶导数的边缘算子
  通过模版作为卷积核与图像的每个像素点做卷积运算，然后选择合适的阈值来提取图像的边缘。常用的有Roberts算子、Sobel算子和Prewitt算子。
（2）二阶导数的边缘算子
  依据与二阶导数过零点，常见的有Laplacian算子，此类算子对噪声敏感。
（3）其他边缘算子
  前面两类算子均通过微分来检测图像边缘，还有一种就是Canny算子，其就是在满足一定约束条件下推到出来的边缘检测最优化算子。

1.3 梯度

1.3.1 图像梯度

  为了达到寻找边缘的目的，检测灰度变化可用一阶导数或者二阶导数来完成。下面讨论一阶导数。
  为了在一副图像$f$的$(x,y)$位置处寻找边缘的强度和方向，所选择的工具就是梯度，梯度用$\nabla f$来表示，并用向量来定义，定义如下所示：
$$\nabla f=grad(f)=\left[\begin{array}{c}{g}_x\\ g_y\end{array}\right]=\left[\begin{array}{c}\frac{\partial f}{\partial x}\\ \\ \frac{\partial f}{\partial x}\end{array}\right]\text{\hspace{1em}(1)}$$
  其中，梯度$\nabla f$为一个向量，他表示$f$在位置$(x,y)$处的最大变化率的方向。
  梯度的大小用$M(x,y)$表示，则：
$$M(x,y)=mag(\nabla f)=\sqrt{g_x^2+y_y^2}\text{\hspace{1em}(2)}$$
  其中，$M(x,y)$表示梯度向量方向变化率的值。
数学梯度的简单推导
  对于以函数$f(x)$在点$x$处的导数近似：将函数$f(x+\Delta x)$展开为$x$的泰勒级数，令$\Delta x=1$，且只保留该级数的线性项，则函数$f(x)$的梯度$\nabla f$计算为：
$$\nabla f=\frac{\partial f}{\partial x}=f^{{}^{\prime }}(x)=f(x+1)-f(x)\text{\hspace{1em}(3)}$$

1.3.2 梯度算子

  由上面的数学推导可知，要得到一副图像的梯度，则要求图像的每一个像素点位置处计算偏导数。我们处理的是数字量，因此需要求关于一点的邻域上的偏导数的数字近似，因此一副图像$f$，在$(x,y)$位置处的$x$和$y$方向上的梯度大小$g_x$和$g_y$分别计算为：
$$\begin{array}{c}{g}_x=\frac{\partial f(x,y)}{\partial x}=f(x+1,y)-f(x,y)\\ \\ g_y=\frac{\partial f(x,y)}{\partial x}=f(x,y+1)-f(x,y)\end{array}\text{\hspace{1em}(4)}$$
  上述公式对所有$x$和$y$的有关值可用下图的一维模版对$f(x,y)$的滤波得到。

  用于计算梯度偏导数的滤波器模版，通常称之为梯度算子、边缘算子和边缘检测算子等。
  对于不同的滤波器模版得到的梯度是不同的，这也就衍生出了很多算子，如Roberts、Prewitt、Sobel和Laplacian算子等。下面将详细介绍不同的算子。

2 Roberts算子

3 Prewitt算子

4 Sobel算子

4.1 基本原理

  Sobel算子是一种用于边缘检测的离散微分算子，他结合了高斯平滑和微分求导。该算子用于计算图像明暗程度近似值，根据图像边缘旁边明暗程度把该区域内超过某个数的特定点记为边缘。Sobel算子再Prewitt算子的基础上增加了权重的概念，认为相邻点的距离远近对当前像素点的影响是不同的，距离越近的点影响越大，从而实现图像锐化并突出边缘轮廓。
  Sobel算子根据像素点上下、左右邻点灰度的加权差，在边缘达到极值这一现象检测边缘。对噪音具有平滑作用，提供较为准确的边缘信息。因为Soble算子结合了高斯平滑和微分求导（分化），因此结果会具有较多的抗噪性，当对精度要求不高时，Sobel算子是一种较为常用的边缘检测算法。
  Soble算子的边缘定位更为准确，常用于噪声较多、灰度渐变的图像。其算法模板如下面的公式所示，其中$d_x$表示水平方向，$d_y$表示垂直方向。
$$d_x=\left[\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ -1& 0& 1\end{array}\right]d_y=\left[\begin{array}{ccc}-1& -2& -1\\ 0& 0& 0\\ 1& 2& 1\end{array}\right]\text{\hspace{1em}(4)}$$
  例如，下面给出Sobel算子的模板，在像素点P5处$x$和$y$方向上的梯度大小$g_x$和$g_y$分别计算为：

$$\begin{array}{c}{g}_x=\frac{\partial f(x,y)}{\partial x}=(P7+2P8+P9)-(P1+2P2+P3)\\ \\ g_y=\frac{\partial f(x,y)}{\partial x}=(P3+2P6+P9)-(P1+2P4+P7)\end{array}\text{\hspace{1em}(4)}$$
  图像中的每一个像素的横向以及纵向灰度值通过以下公式结合，来计算该点的灰度值大小：
$$G=\sqrt{G_x^2+G_y^2}$$
  通常为了提高效率，使用不开平方的近似值，但是这样做会损失精度，迫不得已的时候可以如下这样子：
$$G=\mid G_x\mid +\mid G_y\mid$$
  如果梯度G大于某一阈值，则认为该点$(x,y)$为边缘点。

4.2 Sobel边缘检测的FPGA实现

4.3 Sobel边缘检测的MATLAB算法实现与验证

clc;
clear all ;
GRAY = imread('../img/1920X1080_gray.bmp');
[row,col] = size(GRAY);
sobel_result = zeros(row,col);
sobel_padding = zeros(row+2,col+2);

x_mod = [-1,0,1;-2,0,2;-1,0,1];
y_mod = [-1,-2,-1;0,0,0;1,2,1];

for i = 1:row
    for j = 1:col
        sobel_padding(i+1,j+1) = GRAY(i,j);
    end
end

for i = 1:row+2
    sobel_padding(i,1) = sobel_padding(i,2);
    sobel_padding(i,col+2) = sobel_padding(i,col+1);
end

for i = 1:col+2
   sobel_padding(1,i) = sobel_padding(2,i);
   sobel_padding(row+2,i) = sobel_padding(row+1,i);
end

for i = 2:row+1
    for j = 2:col+1
        matrix11 = sobel_padding(i-1,j-1);
        matrix12 = sobel_padding(i-1,j);
        matrix13 = sobel_padding(i-1,j+1);
        
        matrix21 = sobel_padding(i,j-1);
        matrix22 = sobel_padding(i,j);
        matrix23 = sobel_padding(i,j+1);
        
        matrix31 = sobel_padding(i+1,j-1);
        matrix32 = sobel_padding(i+1,j);
        matrix33 = sobel_padding(i+1,j+1);

        matrix = [matrix11,matrix12,matrix13;matrix21,matrix22,matrix23;matrix31,matrix32,matrix33];

        x_mult = matrix.* x_mod;
        y_mult = matrix.* y_mod;

        gx1 = sum(sum(x_mult()));
        gy1 = sum(sum(x_mult()));

        gx2 = gx1*gx1;
        gy2 = gy1*gy1;
        
        sobel_result(i-1,j-1) = sqrt(gx2+gy2);
    end
end
subplot(1,2,1)
matlab_Y = uint8(floor(sobel_result));
imshow(matlab_Y),title('MATLAB  Sobel算法图像');
subplot(1,2,2)
imshow(GRAY),title('原图像');

5 Laplacian算子

6 Canny算子

7 小结