陈仁爱，凌强，徐骏，李峰

　　(中国科学技术大学 信息科学技术学院，安徽 合肥 230026)

摘要：霍夫圆变换是图像处理中人眼检测的一种常见方法，但是其处理的数据量多，处理速度慢，在移植到DSP上后难以满足实时性要求。对此，提出了一种将两阶段霍夫圆变换算法应用到 TMS320C6000系列 DSP上的实现与优化方法。首先，在算法上对霍夫圆变换使用MarrHildreth算子增强等方法进行改进以保证检测的准确率；之后根据 DSP的特点，利用C代码优化、浮点定点转换和软件流水等技术对算法进行深度优化。实验结果表明，程序的运行时间明显缩短，为视线检测的实时性实现创造了良好的条件。

关键词：C6000 DSP；霍夫变换；程序优化

0引言

　　在视频图像中快速地检测出圆形是目标跟踪、目标分类和行为理解等更高层次视频图像分析的重要基础，比如人眼检测、视线跟踪和交通视频分析等。霍夫圆变换是图像处理中识别和定位圆形的常用方法，其准确率高且与图中形状的方向无关［12］，被广泛用于运动目标轨迹的检测与识别。在视线检测领域，也有众多关于它的研究。MATHEWS R 提出了一种利用霍夫圆变换来设计鼠标的方法［3］，ZIA M A等人尝试利用眼球追踪来为残疾人设计轮椅［4］。

　　标准霍夫变换虽然具有显著的优势，但其不足也不容忽视［56］。它需消耗大量的时空资源，对于在嵌入式平台上进行视线检测这样的应用背景无法做到实时控制。目前有很多研究都仅针对算法上的改进，而没有结合具体的实现平台的特点来进行优化，不能充分利用硬件的性能优势。本文拟在德州仪器的TMS320C6000系列DSP上使用两阶段霍夫变换算法来实现圆形检测，并开展基于DSP平台的程序优化研究。

1改进的两阶段霍夫圆变换算法

　　参考文献［2］提出将霍夫变换一般化的方法，对于任何曲线，只要给出了它的函数方程，就可以利用霍夫变换的方法，将图像空间变换到霍夫参数空间，利用投票的方法求得曲线参数。对于检测圆的情形，由于圆的方程有3个未知量，变换到霍夫空间中需要一个三维的累加器，对于高清图片来说将耗费大量的内存，而且搜索极值时时间代价很大。这二者对于DSP平台都是致命的，故在移植到DSP时本文采用了两阶段霍夫圆变换算法，并对其进行性能上的改进。本节阐述了基于两阶段霍夫圆变换算法的圆检测方案设计。

　1.1两阶段霍夫圆变换算法

　　两阶段霍夫圆变换是一种针对标准霍夫圆变换的参数空间分解的方法，主要目的是为了减少原算法的空间复杂度，其输入是边缘图像［5］。

　　考虑如下圆的参数方程：

　　其中(x0,y0,r)是一组圆心和坐标参数。两阶段霍夫圆变换的第一步是对圆心参数空间累加。根据圆的一阶导数和二阶导数的特性，过圆周上任意一点的圆切线的垂线经过圆心，如图1所示。对已知边缘上的任意点做垂线，这些垂线将会在(a, b)空间汇集，形成一个热点，在(a, b)空间搜索极值即得圆心坐标。给定r的范围，在边缘点上做垂线段，得(a, b)空间。即：

　　其中，(minr,maxr)是给定的半径的范围，也是做出的垂线段的长度，A是(a, b)空间的累加器， E(i, j)是待检测图像的边缘图。在(a, b)空间搜索极值即得圆心坐标。

　　求得圆心坐标之后，在此基础上可以进行半径参数空间的累加。对每一个检测的圆，R 空间累加方式为：

　　其中，E是边缘图，r是给定的半径范围。在R空间搜索极值即可求得半径。

　　1.2圆形检测整体方案设计

　　1.1节描述了一种节省空间开销的霍夫变换方法，基于此，本节设计了一套圆检测方案，增加了预处理和利用MarrHildreth算子进行图像增强等步骤，具体如图2所示。

　　由于形状只与灰度值有关，故首先获取灰度图。之后，对于含噪声的图像，对其进行平滑处理，减少边缘检测的工作量和出错率。平滑操作使用高斯平均算子，模板大小为5，方差为1，实现时采用空域卷积的方式。为了保存更多的边缘信息，本文使用一阶及二阶边缘检测的方式来求取边缘图而不是直接对灰度图进行二值化。一阶边缘提取使用Sobel算子，包括水平Sobel算子和垂直Sobel算子，最后取二者的几何平均作为最终的Sobel图像。在使用圆的方向导数缩减参数空间时需要使用边缘图像的方向信息，本文使用Sobel处理后的图像来直接求得。

　　对于每个被检测到的边缘点(i,j)，计算其方向角度θ(i,j)：

　　其中，f(i,j) 为边缘图像像素值，Sobel(i,j)为Sobel处理后的像素值，脚标v、h分别表示垂直和水平。C语言中atan2求得的角度值在 (-π, π) 间，由于θ 与θ±π的方向相同，可以将θ左右平移到区间 (-π2, π2)中。

　　对于边缘更复杂的图像，当一阶边缘检测过粗或者错误时，使用二阶边缘检测能获得更好的检测效果。MarrHildreth算子是利用高斯滤波的一种二阶滤波方式，它先对图像进行高斯平滑，之后应用拉普拉斯运算，即：

　　其中P是待处理的图像，g(x,y)是高斯平滑滤波器。此外，MarrHildreth算子还被用于图像增强。

　　在获得边缘图像和边缘方向信息之后使用1.1节中的方法就能求得圆心和半径。

　　1.3圆形检测中的实现细节

　　由于边缘检测的不精确性和待检测圆的不确定性，两阶段霍夫圆变换方法在实现的过程中存在一些问题，本节将讨论这些问题的解决方式。

　　使用式（2）、（3）进行(a,b)空间累加时，会得到一些热点，但是由于圆的形变等原因这些热点常常是不清晰的或者弥漫开来的，故有必要将它们聚集成一个更集中的点。使用式(6)MarrHildreth算子进行图像增强，将获得更明亮的热点。对增强后的(a, b)空间进行阈值化处理，得所求圆心。对于原图中有多个圆的情形，不同的阈值将保留不同数量的圆心，数值越大，对所检测圆的要求越高，得到的圆心越少。由于圆的大小不一，得到的热点（即圆心）亮度不一：大的圆在(a, b)空间中比小的圆累积值更大。考虑在变换时添加衰减因子1k，累加公式(3)可变为：

　　实验中发现大的圆由于边缘方向估计误差更大，在圆心区域的偏差会更大，这会使它的聚集程度下降，从而与由半径大带来的优势相抵消，此时k可设为1。

　　在R空间使用式(4)时，累加过程中只利用了垂线经过圆心的边缘点，由于误差的原因导致效果不太好。考虑使用所有的边缘点进行累加。设圆心与边缘点的连线夹角是ψ，则将垂线夹角在 (ψ-∈,ψ+∈) 间的边缘点都加入计算，其中∈是允许误差，一般取π4或π8都能有较好的效果。

　　对于同心圆，在R空间累加后选择合适的阈值以保留多个峰值，能获得不同的半径值。

　　图3是一组检测结果示意图。其中图(c) 是将θ映射到(0,255) 得来。在采取合适的参数后能获得较好的检测结果。

2基于C6000 DSP的移植与优化

　　2.1移植到DSP

　　本文使用的硬件平台是TI的TMS320C64x+TM定点型DSP核，使用的开发环境为CCSv4。使用C64x+ simulator进行软件仿真，通过CCS的时钟工具测得试运行时钟周期，通过profile工具分析工程耗时分布［7］。

　2.2具体的优化步骤

　　本文的优化基于所用DSP的结构特性，尽量充分利用DSP的计算资源来缩短检测时间［89］。本节将介绍本文使用的优化方法。

　　2.2.1基于编译器的优化方法

　　CCS的编译器中设置了众多参数，选择合适的参数能减少检测耗时。

　　选择优化级别。由于本文不考虑代码体积，故使用-o3文件级优化。

　　使用-mt。假设不存在多个指针对同一个内存(块)进行读写操作。

　　使用 -mh。 允许编译器取超过数组边界num字节的值。该操作使编译器在编排软件流水时有额外的弹性，可提高流水性能。在CMD文件中定义大于num的缓存空间，避免 EDMA或其他cache冲突。

　　不使用 -g、-ss等参数。这些参数对调试很有效，但是会造成性能下降，在最终发布产品时不应使用。

　　CCS编译器选项中有部分能减少手工整定时间但是对代码性能没有影响的参数，使用这些参数有助于程序优化。使用-s［-k|-al］-o［2|3］能生成优化后的类似于C的代码，且优化器描述被嵌套在汇编代码中，使用-mw或者-mw-al输出额外的关于软件流水的信息，包括单次循环的调度方式；使用-on2-o3生成*.nfo文件，以给出高级别的优化总结和可用的优化建议。

　　2.2.2手工整定方法

　　仅使用基于编译器的优化所减少的程序耗时是有限的，本文对算法的实时性要求较高，需要在2.2.1节的基础上进行手工整定。手工整定的方式很多，本文使用的几种方法有基于编译器和优化器描述、浮点定点转换、不用除法等。

　　根据优化器描述给所有安全的指针添加restrict关键字，消除循环间的依赖，使生成的汇编文件中循环依赖项值为零；根据汇编嵌入信息添加pragma 指令MUST_ITERATE和UNROLL，告诉编译器循环的最大值、最小值、公约数和展开次数，并根据软件流水信息中各资源的使用情况在循环前使用 _nassert() 告诉编译器数据是64位对齐的，一次读取多个对齐数据以减少D单元和T通道的使用数，以能较好地平衡资源。这些措施极大地提高了软件流水的效率。

　　由于算法涉及众多的浮点运算，在定点型甚至浮点型DSP中效率都很低，故在精度允许的情况下，使用定点运算代替浮点型运算。浮点转换为定点除了手动使用Qn定标外，对于C64x+ DSP还可以使用TI的C64x+ IQmath库。IQmath集合了很多高精度且高度优化的数学函数，适合于将浮点的算法转换成定点［10］。除了提供_iq数据类型及各类型互相转换之外，IQmath还提供了各种高度优化的数学函数比如_IQcos(余弦函数)，并且支持可调节精度，在不同地方可以选择不同的定标Q值，即GLOBAL_Q可以在需要的地方指定0~31中的任意值。比如融合两个Sobel算子时求两个数的几何平均：

　　(*imageTot)［i］［j］= pow( (*imageH)［i］［j］*(*imageH)［i］［j］ +(*imageV)［i］［j］*(*imageV)［i］［j］, 0.5f);

　　可以修改为：

　　tmp=_IQpow(*imageH［i］［j］,2) +_IQpow(*imageV［i］［j］,2);//sqrt(H^2+V^2)

　　*imageTot［i］［j］ = _IQsqrt(tmp);//pow(tmp,0.5f)

　　先将数据转换成_iq类型，然后使用IQmath中的_IQpow和_IQsqrt来求幂和根号，之后再转换回需要的float型。

　　对于除法运算，DSP是利用函数调用来完成，耗费的时钟周期比乘法高出2个数量级以上。为了消除除法，可以使用移位或者查找表（lookup table）来代替。比如:

　　(int)(_abs(im［i］［j］ + shift) /maxval*255);

　　其中maxval值在0~255之间，可以提前将255/maxval计算出来，计算时查表代替；或者使用8替代*255。另外，使用 _IQdiv( _iq A, _iq B) 函数可完成_iq类型的除法。

　　对于DSP的优化方法还有很多，如使用内联函数和线性汇编等，在进一步优化过程中可使用。

　　2.3优化结果与分析

　　在综合了各种优化方式之后，重新使用profile测试各部分的时钟周期，保持其他条件不变。程序优化前后所用的时间见表1。可见经过优化，耗时明显减少。

3结论

　　本文先在PC上实现了基于霍夫变换的圆形检测算法，然后将它移植到了DSP上并进行了基于时间的优化分析，使得实时性有很大的提高。

　　若要将此算法应用到人眼检测，可以在前述的基础上继续改进，比如缩减霍夫变换中r的范围，在R空间寻找峰值时指定目标是2个等。本文所完成的只是视线跟踪中人眼检测的一部分，离实际应用还很远。本文设计还有很多不足，希望在未来的工作中继续改进。

参考文献

　　［1］ BALLARD D H. Generalizing the Hough transform to detect arbitraryshapes［J］. Pattern Recognition, 1981, 13(2): 111122.

　　［2］ YUEN H K, PRINCEN J, ILLINGWORTH J. Comparative study of Hough transform methods for circle finding［J］. Image and Vision Computing, 1990, 8(1): 7177.

　　［3］ MATHEWS R, CHANDRA N. Computer mouse using eyetracking system based on houghman circle detection algorithm with grid analysis［J］. International Journal of Computer Applications, 2012, 40(13): 1216.

　　［4］ ZIA M A, ANSARI U, JAMIL M, et al.Face and eye detection in images using skin color segmentation and circular Hough transform［C］.Robotics and Emerging Allied Technologies in Engineering (iCREATE), 2014 International Conference on. IEEE, 2014: 211213.