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算子

  Roberts算子又称为交叉微分算法，它是基于交叉差分的梯度算法，通过局部差分计算检测边缘线条。常用来处理具有陡峭的地噪声图像，当图像边缘接近正45度或者负45度时，该算法处理效果更理想。缺点是对边缘的定位不大准确，提取的边缘线条比较粗。
  Roberts算子的模板分为水平方向和垂直方向，如下图所示，从其模板可以看出，Roberts算子能较好的增强正负45度的图像边缘。
$$d_x=\left[\begin{array}{cc}-1& 0\\ 0& 1\end{array}\right]d_y=\left[\begin{array}{cc}0& -1\\ 1& 0\end{array}\right]$$
  例如，下面给出Roberts算子的模板，在像素点P5处的$x$和$y$方向上的梯度大小$g_x$和$g_y$分别计算为；

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}=(P3+2P6+P9)-(P1+2P4+P7)\\ \\ g_y=\frac{\partial f(x,y)}{\partial x}=(P7+2P8+P9)-(P1+2P2+P3)\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实现

  3X3矩阵生成后，我们需要同时计算两个矩阵结果，一个是$g_x$，一个是$g_y$，计算后再取绝对值，然后再算乘方并相加，最后将数据写入开方的CORDIC的IP核计算开方后的结果，即可得到Sobel算法的结果，如果想要进一步的边缘检测，可以将得到的数据进行阈值处理，这里我们得到Sobel算法后的图像就不做进一步的处理了,如果有需要可以自己做个阈值分割。

`timescale 1ns / 1ps

module sobel#(
    parameter DW = 8
)(
    input   wire                clk             ,
    input   wire                rst_n           ,

    input   wire                matrix_de       ,   
    input   wire    [DW-1:0]    matrix11        , 
    input   wire    [DW-1:0]    matrix12        ,  
    input   wire    [DW-1:0]    matrix13        ,
    input   wire    [DW-1:0]    matrix21        ,
    input   wire    [DW-1:0]    matrix22        ,
    input   wire    [DW-1:0]    matrix23        ,
    input   wire    [DW-1:0]    matrix31        ,
    input   wire    [DW-1:0]    matrix32        ,
    input   wire    [DW-1:0]    matrix33        ,

    output  wire                sobel_data_de   ,
    output  wire    [DW-1:0]    sobel_data
);

//              gx                      gy                  pix
//      [-1     0       1]      [1      2       1]      [P11    P12     P13]
//      [-2     0       2]      [0      0       0]      [P21    P22     P23]
//      [-1     0       1]      [-1     -2     -1]      [P31    P32     P33]

reg     [15:0]    	x_one_line          ;
reg     [15:0]    	x_two_line          ;
reg     [15:0]    	x_three_line        ;
reg     [15:0]    	x_sum_matrix        ;
reg     [15:0]    	x_abs_sum_matrix    ;
reg     [31:0]   	gx2                 ;

reg     [15:0]    	y_one_line          ;
reg     [15:0]    	y_two_line          ;
reg     [15:0]    	y_three_line        ; 
reg     [15:0]    	y_sum_matrix        ;
reg     [15:0]    	y_abs_sum_matrix    ;
reg     [31:0]  	gy2                 ;

reg     [31:0]      g2                  ;

wire    [23:0]      square_data         ;
wire                square_data_de      ;

reg                 pre_sobel_data_de   ;
reg     [DW-1:0]    pre_sobel_data      ;

reg     [4:0]       matrix_de_r;
always @(posedge clk)begin
    if(rst_n==0)begin
        matrix_de_r   <=     0;  
    end 
    else begin
        matrix_de_r   <=    {matrix_de_r[3:0],matrix_de};         
    end
end

always @(posedge clk)begin
    if(rst_n == 0)begin
        x_one_line      <=  0   ;
        x_two_line      <=  0   ;
        x_three_line    <=  0   ;

        y_one_line      <=  0   ;
        y_two_line      <=  0   ;
        y_three_line    <=  0   ;
    end
    else if(matrix_de)begin
        x_one_line      <=  (-1)*matrix11   +   0*matrix12      +   1*matrix13 ;
        x_two_line      <=  (-2)*matrix21   +   0*matrix22      +   2*matrix23 ;
        x_three_line    <=  (-1)*matrix31   +   0*matrix32      +   1*matrix33 ;

        y_one_line      <=  (1)*matrix11    +   2*matrix12      +   1*matrix13 ;
        y_two_line      <=  0*matrix21      +   0*matrix22      +   0*matrix23 ;
        y_three_line    <=  (-1)*matrix31   -   2*matrix32    	-   1*matrix33 ;
    end
end

always @(posedge clk)begin
    if(rst_n==0)begin
        x_sum_matrix    <=  0;
        y_sum_matrix    <=  0;
    end
    else if(matrix_de_r[0])begin
        x_sum_matrix    <=  x_one_line + x_two_line + x_three_line;
        y_sum_matrix    <=  y_one_line + y_two_line + y_three_line;
    end
end

always @(posedge clk)begin
    if(rst_n==0)begin
        x_abs_sum_matrix    <=  0;
    end
    else if(matrix_de_r[1])begin
        if(x_sum_matrix[15] == 1)begin
            x_abs_sum_matrix    <=  65536 - x_sum_matrix;
        end
        else begin
            x_abs_sum_matrix    <=  x_sum_matrix;
        end
    end
end

always @(posedge clk)begin
    if(rst_n==0)begin
        y_abs_sum_matrix    <=  0;
    end
    else if(matrix_de_r[1])begin
        if(y_sum_matrix[15] == 1)begin
            y_abs_sum_matrix    <=  65536 - y_sum_matrix;
        end
        else begin
            y_abs_sum_matrix    <=  y_sum_matrix;
        end
    end
end

always @(posedge clk)begin
    if(rst_n==0)begin
        gx2     <=  0;
        gy2     <=  0;
    end
    else if(matrix_de_r[2])begin
        gx2     <=  x_abs_sum_matrix * x_abs_sum_matrix;
        gy2     <=  y_abs_sum_matrix * y_abs_sum_matrix; 
    end
end

always @(posedge clk)begin
    if(rst_n==0)begin
        g2  <=  0;
    end 
    else if(matrix_de_r[3])begin
        g2  <=  gx2 + gy2;
    end
end

sqrt u1_sqrt (
    .aclk                           (clk                ), 
    .aresetn                        (rst_n              ),                                                    
    .s_axis_cartesian_tvalid        (matrix_de_r[4]     ),  
    .s_axis_cartesian_tdata         (g2                 ),   
    .m_axis_dout_tvalid             (square_data_de 	),         
    .m_axis_dout_tdata              (square_data        )           
);

always @(posedge clk)begin
	if(rst_n==0)begin
		pre_sobel_data_de 	<= 	1'b0 ;
	end
	else begin
		pre_sobel_data_de 	<= 	square_data_de 	;
	end
end

always @(posedge clk)begin
	if(rst_n==0)begin
		pre_sobel_data 	<= 	0;
	end
	else if(square_data > 255)begin
		pre_sobel_data 	<= 	255;
	end
	else begin
		pre_sobel_data 	<= 	square_data;
	end
end

assign 	sobel_data 		= 	pre_sobel_data 		;
assign 	sobel_data_de 	= 	pre_sobel_data_de	;

endmodule

  其中例化了一个XILINX官方IP核CORDIC，IP核的配置如下图：

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

clear;
clc;
close all;
sigma = 0.8 ;
A = exp(-(1+1)/(2*sigma*sigma))/(2*pi*sigma*sigma);
B = exp(-(1+0)/(2*sigma*sigma))/(2*pi*sigma*sigma);
C = exp(-(0+0)/(2*sigma*sigma))/(2*pi*sigma*sigma);
D = A*4 + B*4 + C;

gauss_double = [A,B,A;B,C,B;A,B,A];
gauss_normal = gauss_double / sum(sum(gauss_double));
gauss_integer = floor(gauss_normal/gauss_normal(1,1));

GRAY = imread('../img/gray.bmp');
[row,col] = size(GRAY);
gassin_padding  =   zeros(row+2,col+2);
gassin_result   =   zeros(row,col);

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

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

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

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

        matrix = [matrix11,matrix12,matrix13;matrix21,matrix22,matrix23;matrix31,matrix32,matrix33];
        gassin_mult = matrix.* gauss_integer;

        sum_gassin_matrix = sum(sum(gassin_mult()));
        gassin_result(i-1,j-1) = sum_gassin_matrix/16;
    end
end

a = textread('../data/gassin_filter.txt','%s');
IMdec1 = hex2dec(a);

IM1 = reshape(IMdec1,col,row);
fpga_Y = uint8(IM1)';

b = textread('../data/pre.txt','%s');

subplot(1,3,1)
matlab_Y = uint8(floor(gassin_result));
imshow(matlab_Y),title('MATLAB gassin算法图像');
subplot(1,3,2)
imshow(fpga_Y),title('FPGA gassin算法图像');
subplot(1,3,3)
imshow(GRAY),title('原图像');

sub = matlab_Y - fpga_Y;

min_sub = min(min(sub));
max_sub = max(max(sub));

5 Laplacian算子

6 Canny算子

7 小结