一种组合GAUSS-filter、SOBEL、NMS、OTSU 4种算法的图像边缘检测的FPGA实现

张 灏，王素珍，郑 宇，王 伟，李潍志，王 鹏，任贵珊，马家麟

(青岛大学 电子信息学院，山东 青岛 266071)

1 引 言

边缘检测是分离图像目标主体和背景的主要依据。通过检测图像上灰度变化显著的像素点集合来确定图像的边缘。图像的边缘检测在许多领域，比如医学影像、地质勘探、卫星遥感、航天、交通、计算机视觉等方面都有重要的应用。目前对图像清晰度的要求越来越高，图像处理的信息量巨大。一般串行方式已经不能满足对图像大量信息的实时处理。大规模可编程片上系统(Field Programmable Gate Array，FPGA)芯片以其独特的并行可编程处理信号的方式，成为处理图像的重要器件。

边缘检测的常用方法是SOBEL算法。文献[1-5]成功地在FPGA芯片上应用SOBEL算法实现了图像的边缘检测。但是，其边缘检测的门限阈值需要反复的实验确定，且图像有噪声的影响，同时边缘轮廓也不够细化。文献[6]引入高斯滤波(GAUSS-filter)，对图像先去噪声处理，再送入SOBEL算法进行边缘检测，减少了图像边缘背景噪声的干扰。文献[7]把经过SOBEL算法检测出的边缘图像送到非极大值抑制(Non-maximum Suppression，NMS)算法模块处理，细化了图像轮廓的边缘。文献[8-10]将大津法OTSU加入SOBEL算法中，使图像边缘的门限阈值能够自适应图像特点调整。

本文在FPGA器件上，综合应用GAUSS-filter、SOBEL、NMS、OTSU等多种算法，进行图像边缘检测。将原始采集到的彩色图像转换成灰度图像，先对图像进行高斯平滑处理，消除背景干扰性噪声；应用SOBEL算法计算灰度梯度值，得到图像边缘；再应用NMS算法处理灰度，对比出最大的灰度像素集合；最后应用OTSU算法自适应图像特点确定图像门限阈值，提高图像边缘检测的适应性。实验结果证明本文方法对不同图像都有很好的边缘检测效果。在运行速度上，比软件实现时间少很多，提高了系统处理的速度。该实时图像边缘检测方法能够嵌入式地应用到许多领域中。

2 基于FPGA的实时边缘检测系统

本系统采用XILINX公司的ARTIX7系列的xc7a100tfgg484-2芯片作为边缘检测的主要器件，同时在FPGA主控芯片外围配置双倍速率同步动态随机存储器DDR3 SDRAM(Double Data Rate SDRAM)缓存图像。基于FPGA的实时边缘检测系统如图1所示。采用乒乓操作方法将来自CCD摄像机的分辨率为640×480的图像先缓存到DDR3存储器中，在图像请求处理信号的作用下，图像数据流将从DDR3存储中被读出，进行图像的边缘检测处理。原始采集到的图像为彩色图，将图像由彩色空间转换到灰度空间；转换后的灰度图像混有噪声，在进行图像边缘检测之前，先将图像送到高斯滤波器中处理，消除图像中的噪声。

图1 整体结构框图

高斯滤波处理后的图像送到SOBEL算法模块中，进行梯度的计算，包括X、Y方向的梯度计算，进而计算出像素的梯度值，这些梯度像素点对应的轮廓比较粗糙。为一步对图像边缘进行细化处理，要对梯度值进行非极大值抑制处理，即判断当前像素点的梯度值是否是梯度方向上的极大值，如果是，则保留该梯度值；否则，则将其归零。图像的边缘检测需要确定阈值，采用OTSU算法，求取目标与背景之间的灰度类间方差来确定阈值。图像的梯度值与阈值进行比较，得到图像边缘，将检测到的边缘轮廓送到显示器进行显示。

3 算法原理

3.1 SOBEL算法

SOBEL算法是通过SOBEL算子矩阵与图像邻域矩阵进行卷积，计算X方向和Y方向的梯度分量，来获取图像的边缘信息。图2为边缘检测模板。

(a)3×3 neighborhood window (b)Edge detection openator

图2 边缘检测模板

Fig.2 Edge detection template

图2(a)为中心像素点的3×3邻域窗口，图2(b)为SOBEL算法水平方向、垂直方向的梯度算子。将像素邻域与X、Y方向梯度算子进行卷积运算，分别可以得到水平方向的梯度值Gx、垂直方向的梯度值Gy为

(1)

3×3邻域的梯度卷积值为

(2)

其中G为当前邻域内的梯度值。因为FPGA不容易执行开平方的操作，为了方便FPGA内部电路的计算，此处将邻域的梯度卷积值简化表示为

G=|Gx|+|Gy|.

(3)

3.2 NMS算法

用GAUSS filter和SOBEL算法进行边缘检测可以得到粗线条的图像轮廓。对图像轮廓细化处理的算法是非极大值抑制NMS。NMS是将邻域中心像素点的梯度值与它所在梯度方向上的梯度值的线性插值相比较，判断其是否为最大值，如果是，则保留该像素点梯度值；否则，直接将邻域中心像素点的梯度值归为0。这种方法能够去除图像中存在的伪边缘像素点[11]。

非极大值抑制必须要进行梯度方向的判断，梯度方向计算公式一般为θ(i,j)=arctan(Gx/Gy)，使用FPGA实现反正切函数相对较难。在论文中是通过判断梯度值在X方向Gx的分量和在Y方向的Gy分量之间大小，以及两者的正负来判断当前梯度方向在哪一个“方向象限”，这样去处理当前的梯度方向会减少一定的计算量。

(a)θ(i,j)∈[0°,45°) (b)θ(i,j)∈[45°,90°)

(c)θ(i,j)∈[90°,135°) (d)θ(i,j)∈[135°,180°)

如图3所示，其中Gi为3×3邻域中的像素点的梯度值，mi为中心像素点G4的梯度方向与邻域的交点。图中只画出[0°,180°]上4个“方向象限”内的梯度方向图，[180°,360°]上梯度方向由此类推。梯度方向上的值可以由邻域内相应的梯度值进行线性插值求得。例如：图3(b)的m4点的值可以由G6和G7进行线性插值[12-13]求出。G4所有梯度方向上的线性插值公式如表1所示，其中x和y分别代表|Gx|、|Gy|。

表1 线性插值公式表

从表1中可以看出，通过线性插值公式可得到[0°,360°]梯度方向上的mi值，NMS方法是G4与mi的值进行比较，如果G4比mi的值大，则保留G4;否则将G4置为0。这样经过NMS处理后的任意像素点的梯度值可表示为公式(4)，其中KnmsG4为非极大值抑制结果。

(4)

3.3 OTSU算法

设置一幅图像的灰度阈值就能进行边缘检测。如果设置阈值过大，会丢失一些边缘信息；如果设置的阈值较小，那么会出现一些假的边缘，使得轮廓不准确。为了使边缘最大程度地体现目标信息，就得恰当地确定灰度阈值。本文使用OTSU算法自适应图像灰度来确定当前的阈值。OTSU算法通过计算背景和目标灰度之间的类间方差来区分图像边缘，类间方差越大，证明背景和目标之间的差别越大[14-16]。

假设w(D0)为背景D0出现的概率，w(D1)为目标D1出现的概率，u(D0)和u(D1)分别为D0和D1的灰度均值。w(D0)、w(D1)、u(D0)以及u(D1)计算公式如公式(5)所示：

(5)

其中：P(i)为灰度级i出现的概率，T为背景和目标的分界灰度级，L为最大灰度级。则整幅图像的背景与目标之间的类间方差g为：

g=w(D0)*w(D1)*(u(D0)-u(D1))2，

(6)

其中：w(D1)=1-w(D0)，将公式(5)中变量代入公式(6)可得：

(7)

其中N(i)为灰度级i的统计数。不同的图像对应不同的轮廓边缘，其对应的最大类间方差不同，当每幅图像的g达到最大时，此时的灰度级i即为我们所需要的最佳阈值。

4 FPGA实现

4.1 基于FPGA的SOBEL算法实现

在FPGA芯片上实现SOBEL算法的过程中，产生3×3邻域领域滑动窗口是关键技术。本文采用XILINX公司的SHIFT_RAM模块进行3×3邻域窗口设计，其深度设置为一幅图像的一行像素点数，宽度为8位。

图4 水平方向和垂直方向梯度计算

如图4所示，构建3×3邻域需要2个SHIFT_RAM缓存模块，每行3个D触发器，共3行。3×3邻域窗口滑窗式地遍历整幅图像实现SOBEL算法，从3×3邻域窗输出的数据对应图2(a)中的像素灰度值。图4右是实现公式(1)的结构框图，计算输出水平方向Gx和垂直方向Gy。图5是实现公式(3)的电路结构图，判断X、Y方向梯度值正负，再分别相加计算SOBEL算法的梯度值。

图5 SOBEL算子梯度值计算

4.2 基于FPGA的NMS算法实现

4.2.1 NMS公式改进

由表1可看出，在进行线性插值计算时，会涉及除法，因为FPGA计算除法需要消耗较多的资源，本文将除法转化为乘法，实现插值计算。将“方向象限”线性插值公式两边同时乘以分母，变除法为乘法运算。例如线性插值m1的公式可表示为xm1=yG2+(x-y)G5。在m1的计算中，x>y，取Mbig=max(x,y)，Msmall=min(x,y)，则有：

Mbig*m1=Msmall*G2+(Mbig-Msmall)*G5，

(8)

同理，在m2～m8的计算中，通过比较x、y之间的大小关系，将线性插值公式总结为式(9)：

Mbig*mi=Msmall*GA+(Mbig-Msmall)*GB，

(9)

其中GA和GB分别为相应的3×3邻域内的两个梯度值。相应地在非极大值抑制判断表达式中进行改变，可得改进后的判断表达式：

(10)

4.2.2 方向象限编码

在FPGA实现公式(9)的计算中，首先通过判断Gx和Gy的大小以及正负来确定邻域内的中心像素点的梯度方向，为此用b2b1b0对8个“方向象限”方向进行编码，用info_data表示，便于程序运行时进行查找表计算。可得到表2所示内容。

表2 方向象限编码

如表2所示，“方向象限”编码信号info_data的b2位为Gy的符号位，b1位为Gx的符号位，b0位是两者之间的大小关系，如果|Gx|>|Gy|，则为0，反之，则为1。例如，第5、6“方向象限”，其Gx、Gy均为负值，则info_data的b2、b1位都为1，但确定梯度方向具体在二者中哪个，则要根据Gx和Gy的大小进行判断，第5“方向象限”|Gx|>|Gy|，“方向象限”编码为110，第6“方向象限”|Gy|>|Gx|，“方向象限”编码为111。根据表2中的“方向象限”信息看出，每个“方向象限”的Gx和Gy大小情况不同，为统一计算，将其中较大的值赋值给Mbig，较小的值赋值给Msmall。图6为“方向象限”信息仿真图，out_x_data=Mbig为梯度方向上梯度值的较大值，out_x_data=Msmall为梯度方向上梯度值的较小值，out_info=info_data为“方向象限”编码信息。

图6 “方向象限”信息仿真图

4.2.3 NMS内插公式的FPGA实现

图7为非极大值实现框图，上半部分为非极大值抑制部分，下半部分为3×3邻域生成部分。要进行非极大值抑制计算，需要构建3×3邻域窗口，窗口构建仿照SOBEL算法模块实现。进入窗口的数据是经过“方向象限”信息处理后的Mbig、Msmall、out_info与Gi拼接的数据。经3×3邻域窗口缓存后，分别提取Mbig、Msmall、out_info与Gi对应的邻域中心点的值data_max，data_min,out_info_4以及G4，将这些数据送入非极大值抑制模块进行非极大值抑制计算。将“方向象限”编码值out_info_4作为LUT查找表的地址，进而寻找到线性插值的3×3邻域梯度值数据GA和GB。例如，当3×3滑动窗口邻域的中心梯度值的out_info为001时，插值公式对应梯度值数据分别为G2、G1、G6、G7，则梯度插值索引值{c3,c2,c1,c0}={2,1,6,7}，根据索引值也可找到梯度值数据。图中减法器用来计算Mbig与Msmall的差值Mbig-Msmall，乘法器主要负责公式(9)中的乘法操作，加法器实现式公式(9)左端两数据的相加，比较器COM、与门&、二选一数据选择器MUX则是用来实现公式(10)非极大值的判断，使用D触发器实现数据延时使数据之间进行对齐。DOUT为非极大值抑制结果，且nmsG4=DOUT。

图7 非极大值实现原理图

4.3 OTSU全局阈值分割

要实现OTSU全局阈值分割，首先是要统计当前灰度级所对应的像素点数；对于灰度级的统计，采用直方图统计方法来进行。送入到直方图统计的数据是经过非极大值抑制后的数据DOUT。直方图统计主要由双口RAM及检测控制器组成，其中双口RAM位宽设置为8位，深度为32位，双口RAM用于存储每个灰度级所对应的统计值。直方图统计的思想如下：利用每个像素点所对应的灰度值作为双口RAM的存储及读取地址，当有当前图像的梯度值输入时，双口RAM中相应地址内的统计值会被读出与递增计数器数据进行相加，直到输入的梯度值不是该地址时，会将此梯度值对应的统计值缓存到RAM，并进行下一个梯度值的统计。

图8 直方图时序设计图

直方图统计的时序设计图如图8所示。对进入直方图统计模块的灰度值数据DIN延迟得到DIN_1、DIN_2。在论文中，控制器中会有一个递增计数器CNT，递增计数器使能信号INC_EN拉高的条件是DIN_1与DIN_2相等，使得CNT递增，用来计算连续灰度值的个数；否则该计数器对应的复位使能信号RST_CNT拉高，CNT置为1，不再计数。对于涉及到的双口RAM，B端口的地址address_b为DIN_1，A端口的地址address_a为DIN_2，两者之间存在时钟延迟是为了让B口读出的统计值q_b刚好与A口地址对应的递增计数器CNT对齐。A端口的写使能WR_A与B端口的读使能RD_B拉高的条件是DIN_1与DIN_2数据不相等，此时需要将前一灰度值已统计好的数据CNT_VLE写入RAM，CNT_VLE为q_b与CNT的和，同时将当前灰度值的统计值q_b从RAM中读出，进行当前灰度值的统计计算。例如，当DIN为222时，延迟后在DIN_1、DIN_2上分别出现灰度值222，因为灰度值222与前一灰度值76不相等，需要将前一灰度值76的统计数据CNT_VLE写入RAM，将上一次灰度值222出现时已统计好的数据q_b读出，同时在INC_EN拉高的控制下q_b与CNT相加，将统计值CNT_VLE从38变到42。

图9 最大类间方差实现图

直方图统计完成后，每一个灰度级所对应的统计值N(i)就存储在双口RAM中对应的地址中，背景和目标之间的最大类间方差需要使用这些统计值进行计算。在类间方差计算时，按照地址从0～255的顺序将梯度的统计值读出，将梯度值与梯度值的统计值进行对齐，方便两者的相乘运算。图9为OTSU算法的公式(7)的计算框图。为了在FPGA中实现OTSU算法过程中省去两个除法器，此处去掉公式(7)右半部分的除数N×N，这样处理不会影响类间方差之间的比较，节约了资源。图10为OTSU计算阈值仿真图。

图10 OTSU阈值仿真图

5 测试结果分析

如图11所示为在FPGA芯片上实现的未进行高斯滤波的SOEBL算法、进行高斯滤波后的不同阈值下的SOBEL算法、滤波后SOBEL和OTSU算法结合、滤波后的SOBEL和NMS算法结合、综合GAUSS-filter、SOBEL、NMS、OTSU算法后的图像边缘检测比较结果。

(a)实际拍摄的人物图片处理效果

(b)实际拍摄的物体结构图片处理效果

(c)丽达图片处理效果

图11(a)是对实际拍摄的人物图片处理的效果；图11(b)是对实际拍摄的物体结构图片处理的效果；图11(c)是对画面层次丰富的丽达图片处理的效果。其中A图为原始灰度图;B图为未进行高斯滤波的SOEBL算法边缘检测效果图，因为未进行高斯滤波，所以有很多噪声存在图像中；C图、D图、E图为只用SOBEL算法实现的不同阈值下的边缘检测效果，其中C图是因为阈值过小，导致过多非图像边缘信息显示出来的效果；D图为阈值过大，图像边缘不能全部显示的效果，轮廓不连续；E图是经过多次实验找出的恰当阈值下所显示图像边缘。由此可见，单一的SOBEL算法需要反复测试来确定一个恰当的阈值，而且不同的图像需要不同的阈值，不便于对图像边缘的快速检测。F图为SOBEL和OTSU算法结合的图像边缘检测效果，克服了多次寻找阈值的缺点；G图为SOEBL算子与NMS算法的结合，将E图像边缘进行细化,但对图像缺少自适应性；H图为综合GAUSS-filter、SOBEL、NMS、OTSU算法后的图像边缘检测结果。比较图11中的3组图片的边缘检测效果可得，综合应用4种算法的H图不仅能够画出图像边缘，而且还能对边缘进行细化，对不同图像都有好的边缘检测效果。

表3 运行时间对比

经过多次实验，得出论文中使用FPGA实现边缘检测的效果与软件实现效果差别不大，但是在实现时间上两者有较大差别。表3是图像处理软、硬件运行的平均时间比较，从表3看出，使用FPGA实现边缘检测比软件实现时间少很多。

文中提到的硬件实现，是FPGA设计软件上的功能验证，将图像处理过程在硬件仿真工具modelsim上实现；提取txt文件，然后在软件工作环境QT上进行图像处理效果比较。

FPGA硬件实现时间仅与工作时钟频率有关。文中软件实现与硬件实现速度的比较是指图像处理过程分别在硬件仿真工具modelsim上运行和在软件工作环境QT运行时间的比较。

在实验中我们使用的PC机配置参数为：64位 win7操作系统，运行内存为6 G，主频为2.7 GHz。FPGA的晶振为200 MHz，其中图像工作频率为40 MHz。

6 结 论

在FPGA上综合应用GAUSS-filter、SOBEL、NMS、OTSU等多种算法，进行图像边缘检测。将原始采集到的彩色图像转换成灰度图像，先对图像进行高斯平滑处理，消除干扰性噪声；应用SOBEL算法计算图像梯度值，再应用NMS算法处理灰度，挑出轮廓边缘梯度值最大的灰度像素集合，最后应用OTSU算法自适应地确定图像门限阈值进行边缘检测。单一SOBEL算法需要人为设置阈值，阈值容易出现过大或者过小的现象，导致图像出现假边缘或者不连续边缘；SOBEL和OTSU算法结合又不能对图像进行细化处理。4种算法综合的检测不仅准确找到图像边缘，而且还能对边缘进行细化，弥补了两种边缘检测算法的不足。实验结果显示对不同的图像都有很好的检测效果。文中所提出的在FPGA上实现的图像边缘检测方法具有准切、实时、高速、嵌入式等特点，能够应用到许多领域中。