基于最小二乘的圆及椭圆检测算法

基于最小二乘原理的圆及椭圆检测算法*孔兵王昭谭玉山（西安交通大学激光与红外应用研究所，西安，710049,[e-mail]）摘要在光学测量中，圆或椭圆检测检测是经常用到的一项关键技术。检测算法的精度、速度直接影响了光学测量的精度及速度，而传统的检测算法如重心法、Hough变换法等在检测精度或速度上存在不足之处。本文首次提出的基于最小二乘原理的圆及椭圆检测算法，达到亚像素级的定位精度，而且还具有很快的计算速度，可适用于实时的光学测量。关键词最小二乘算法图像处理圆椭圆中图分类号TP391TheCircleandEllipseDetectionAlgorithmBasedonLeastSquareMethodKONGBingWANGZhaoTANYu-Shan(InstituteofLaser&InfraredTechnologyApplication,Xi’anJiaotongUniversity,Xi’an,710049)AbstractThecircleandellipsedetectionisthekeytechnique,whichisalwaysusedintheopticalmeasurement.Theprecisionandspeedofthedetectionalgorithminfluencethoseoftheopticalmeasurementsystemdirectly.Thetraditionalalgorithmssuchasgravitymodel,Houghtransformareunsatisfactoryinsomeconditions.Thecircleandellipsedetectionalgorithmbasedontheleastsquaremethodisfirstlyreferredinthispaper.Theorientationprecisionisintheorderofinferiorpixels,andthespeedisfastwiththealgorithm.Thealgorithmissuitableforthereal-timeopticalmeasurement.KeywordsLSM(leastsquaremethod),digitalimageprocessing,circle,ellipse0引言激光光斑中心检测在激光扫描三角法、激光准直仪、激光光斑分析仪等光学测量、检测手段中是一项关键技术[1,2]，检测算法的精度、速度直接影响了光学测量的精度及速度。传统的光斑中心检测算法有重心法、中值法，以及Hough变换法[1]。前两种算法要求光斑图像分布比较均匀，否则将会参生较大误差。后一种算法需逐点投票、记录，所用时间较多，而且精度也不够高。然而由于在实际光学测量中，由于存在的散斑、被测物面反射特性不均匀以及光学系统的影响，导致光斑信号强度分布极不均匀，而且测量中一般对实时性要求较高，采用上述算法均有其不足之处。本文首次提出了基于最小二乘原理的圆及椭圆检测算法，可以同时检测光斑中心及半径（或长、短轴），达到亚像素级的定位精度，而且还具有毫秒级的计算速度，可适用于实时的光学测量。1传统圆检测算法1．1重心法以图1为例进行分析，假设光斑图像处于二维平面坐标系中，大小为NM，光斑图像*受自然科学基金(60077031)及西安交通大学在职博士基金资助是经过预处理后得到二值图像（下同），图中较亮的区域即代表了激光光斑，可表示为01),(jig)()(背景光斑（1）重心法计算的光斑中心),(00yx为MiNjMiNjMiNjMiNjjigjigiyjigjigjx1111011110),(),(),(),(（2）以时间复杂度来考虑算法的速度，假设光斑直径为n，以下均做相同的假设，（2）式是在光斑区域内求和，因此时间复杂度为)(2nO[3]。该算法简单明了，计算速度较快，在光斑光强比较均匀的情况下（对应的图1中光斑形状比较规则）不失为一种好的算法。但是该算法受光斑形状影响比较大，而且只能获取光斑的中心不能检测半径，在某些需要计算光斑半径的测量中不能适用。1．2Hough变换法采用Hough变换检测任意曲线的原理如下[4]：检测曲线的参数方程记为),;,,,(121yxaaafann（3）其中，naaa,,,21为方程参数，yx、为空间图像点坐标。对于图像中任一空间点),(00yx，可由（3）式变换为参数空间),,,(21naaa中的一条曲线。对图像曲线上n个点进行上述变换，在参数空间得到n条曲线，由（3）式可知这n条曲线必定经过同一点02010,,,naaa，根据参数空间的此点坐标便可确定图像空间域中的曲线l。直线、圆的参数方程分别为：sincosyx（4）22)()(byaxr（5）Hough变换是将空间域内每个轮廓点带入参数方程（3），根据计算结果对参数空间),,,(21naaa中的量化点按就近原则进行投票，得票最多的点既为所求图像空间域中曲线对应的参数空间点。由（5）式，圆的参数空间为),,(rba，其中),(ba表示圆心，r表示半径，因此采用Hough变换可以检测出激光光斑的中心及半径。Hough变换需要对参数空间离散化，限制了检测精度，另外参数空间得票最多的点未必唯一，选择不同的点得到的图像空间曲线差异比较大。圆的Hough变换由于对每一个边界点都需要在三维参数空间内逐点投票、记录，时间复杂度为)(4nO，计算时间比较长，而且占用计算机内存比较大，因此在实用中受到了限制。2基于最小二乘原理的圆及椭圆检测算法2．1圆检测算法xy图1光斑图像基于圆拟合的激光光斑中心检测算法，根据最小二乘原理（残差平方和最小）用圆来逼近激光光斑轮廓。圆的方程为：222)()(rbyax（6）在此，取残差为：222)()(rbyaxiii（7）其中，Ei，E表示所有边界的集合，),(iiyx表示图像边界点坐标。残差平方和函数为：22222)()(EiiiEiirbyaxQ（8）根据最小二乘原理[5]，应有0rQbQaQ（9）即0)2()()(20))(2()()(20))(2()()(2222222222rrbyaxrQbyrbyaxbQaxrbyaxaQEiiiiEiiiiEiii（10）将（10）式简化整理得：02202202232222223222222222yyxrybybyaxyayxyxrxbxybxaxaxyxrbybaxa（11）其中各参数可用下式表示：EiEiniminmyxyx1（12）对（11）式消掉二次项后整理为：)(21)()()(21)()(322222232222yyxyyyxbyyaxyyxxyxyxxxbxyyxaxx（13）由上式便可推出参数ba,的表达式，结合（11）式得圆参数为：222222222232222322222222322222232222)(2))((2))(())(()(2))((2))(())((yxbybaxarxyyxyyxxxyyxxyxyxxxxxyyxyyyxbxyyxyyxxxyyxyyxyyyxyyxyxyxxxa（14）由（14）式可以看出，根据最小二乘原理的圆拟合推导出的光斑中心（及半径）检测算法虽然形式复杂，但仅对边界点循环一次就可计算出各参数，时间复杂度为)(nO，较为复杂的开根方运算只是在计算出中心参数ba,后计算半径时仅计算一次，因此整个算法的计算速度将会很快。2．2椭圆检测算法椭圆方程标准形式为：1)()(2222byax（15）将其变形为：dbycax22)()(（16）其中222,cd同样取残差dbycaxi22)()(（17）残差平方和函数为：2222)()(EiiiEiidbycaxQ（18）由最小二乘原理有0dQcQbQaQ（19）即0)1()()(20)()()(20))(2()()(20))(2()()(2222222222EiiiiEiiiiEiiiiEiiidbycaxdQbydbycaxcQbycdbycaxbQaxdbycaxaQ（20）将（10）式简化整理得：00003032332303431230202223220242122010121321014112100002032000401200decebcecbeeaeaedecebcecbeeaeaedecebcecbeeaeaedecebcecbeeaeae（21）其中，223433343223123022422332221203141321221110204032020100222222221yxeyeyexyeyeyxeyeyexyeyexexyexyexexexeyeyexee（22）将（21）式消掉二次项a2，及三次项cb2得000323334312223242112131411cfbcffafcfbcffafcfbcffaf（23）其中，jiijijeeeef0000)4,3,2,1,0;3,2,1(ji（24）记bcz，（23）式化为：342414333231232221131211fffzcafffffffff（25）由上式可解得zca,,的值，也就计算出b的值czb/，进而由（21）式计算d的值0002032000401200/)(ecebcecbeeaeaed（26）由dc,也就可以计算椭圆长、短轴的值cdd/（27）同样，仅对边界点循环一次就可计算出椭圆各参数，时间复杂度为)(nO，算法具有很快的计算速度。3通过迭代法进一步提高检测精度获得圆或椭圆的参数后，便可根据（7）式或（17）式获得各边界点残差，只要去掉残差较大的边界点，保留残差较小的边界点，重新按照（14）式或（25）、（26）、（27）式计算圆或椭圆的参数，经多次迭代便可进一步提高检测精度。在此可设定平均残差平方和作为阈值，平均残差平方和由下式求取EiQQ1(28)其中，残差平方和由（8）式或（18）式求取。4实验生成一幅人工图像，图中的圆边界加入了干扰，原始中心为（344，288），如图2中小十字A所示，半径r为199。在PentiumII266MHz计算机上分别对重心法及本方法进行了比较，重心法获取的中心坐标位置为（336.7,289.6），如图2中小十字线B所示，不能检测半径，所用计算时间为11ms。基于圆拟合的方法采用上述思想进行了两次迭代，检测的中心坐标为(341.0,287.9)，如图2中灰十字线C所示，半径r为195.2，用灰色标记出检测出的圆，计算时间为3.3ms。生成另一幅人工图像，图中的椭圆边界同样加入了干扰，原始中心为（244，220），在图3中以小十字所示，长短轴分别为112,197ba。图a为没有经过迭代直接根据最小二乘法检测出的椭圆，图b为经过一次迭代检测出的椭圆，图c为经过两次迭代检测出