圆形目标的圆心定位方法

1.一种圆形目标的圆心定位方法，其特征在于：
第1步：对图像中圆形目标的圆心和半径进行粗定位，得到圆形目标的粗定位圆心(xoc，yoc)和粗定位圆半径Rc，其中，(xoc，yoc)为图像坐标系下圆心的坐标，
第2步：对图像中的圆形目标进行像素级边缘定位，具体步骤如下：
第2.1步：根据粗定位圆心和粗定位圆半径，从图像中提取一个正方形区域，称为“感兴趣区域”，该正方形区域的提取方法为：以图像中圆形目标的粗定位圆心作为正方形区域的中心点，以圆形目标的粗定位直径加上2～6个像素点作为正方形区域的边长，
第2.2步：用canny算子对图像中所提取的感兴趣区域进行像素级边缘检测，得到圆形目标的像素级边缘，再对canny算子提取的圆形目标的像素级边缘点进行边界跟踪，得到每个像素级边缘点的坐标，且按顺时针顺序存储每个像素级边缘点的坐标，
第3步：对图像中的圆形目标进行亚像素边缘定位，具体步骤如下：
第3.1步：将圆形目标的边缘像素点进行区域划分，划分方法如下：称以圆形目标的粗定位圆心为原点、以水平向右为正方向、以像素为单位长度的直线为x轴；称以圆形目标的粗定位圆心为原点、以垂直向上为正方向、以像素为单位长度的直线为y轴；以粗定位圆心为旋转点，由x轴正方向逆时针旋转的角度为θ，根据θ的大小将圆形目标的边缘像素点划分为两个部分，一部分是θ∈[0°，45°]∪[135°，225°]∪[315°，360°]所对应的边缘像素点，称为第1区域的边缘像素点，另一部分是θ∈(45°，135°)∪(225°，315°)对应的边缘像素点，称为第2区域的边缘像素点，粗定位圆心与边缘像素点的连线为该边缘像素点的梯度方向，
第3.2步：求取边缘像素点沿梯度方向的邻近点，(xp，yp)为以像素为单位的边缘点P的坐标值，具体方法如下：若边缘像素点P(xp，yp)在第1区域内，取该像素点的梯度方向直线与直线x＝xp+1的交点A1(xa1，ya1)为该边缘像素点沿梯度方向的右边第一邻近点；取该像素点的梯度方向直线与直线x＝xp+2的交点B1(xb1，yb1)为该边缘像素点沿梯度方向的右边第二邻近点；取该像素点的梯度方向直线与直线x＝xp+3的交点F1(xf1，yf1)为该边缘像素点沿梯度方向的右边第三邻近点；取该像素点的梯度方向直线与直线x＝xp-1的交点C1(xc1，yc1)为该边缘像素点沿梯度方向的左边第一邻近点；取该像素点的梯度方向直线与直线x＝xp-2的交点D1(xd1，yd1)为该边缘像素点沿梯度方向的左边第二邻近点，取该像素点的梯度方向直线与直线x＝xp-3的交点E1(xe1，ye1)为该边缘像素点沿梯度方向的左边第三邻近点；若边缘像素点P(xp，yp)在第2区域内，取该像素点的梯度方向直线与直线y＝yp+1的交点A2(xa2，ya2)为该边缘像素点沿梯度方向的上边第一邻近点；取该像素点的梯度方向直线与直线y＝yp+2的交点B2(xb2，yb2)为该边缘像素点沿梯度方向的上边第二邻近点；取该像素点的梯度方向直线与直线y＝yp+3的交点F2(xf2，yf2)为该边缘像素点沿梯度方向的上边第三邻近点；取该像素点的梯度方向直线与直线y＝yp-1的交点C2(xc2，yc2)为该边缘像素点沿梯度方向的下边第一邻近点；取该像素点的梯度方向直线与直线y＝yp-2的交点D2(xd2，yd2)为该边缘像素点沿梯度方向的下边第二邻近点；取该像素点的梯度方向直线与直线y＝yp-3的交点E2(xe2，ye2)为该边缘像素点沿梯度方向的下边第三邻近点，
第3.3步：采用线性灰度插值的方法求取以像素为单位的边缘像素点P(xp，yp)沿梯度方向的邻近点的灰度值，其中，(xp，yp)为以像素为单位的边缘像素点P的坐标值，使用f(x，y)表示坐标为(x，y)的像素点的灰度值，使用符号[]表示取整数部分，所述的用于获取邻近点灰度值的线性灰度插值方法如下：
(1)若边缘像素点P(xp，yp)在第1区域内，邻近点A1(xa1，ya1)的灰度值
f(xa1，ya1)＝(1-λ)*f(xa1，[ya1])+λ*f(xa1，[ya1]+1)
λ＝ya1-[ya1]，
邻近点B1(xb1，yb1)的灰度值
f(xb1，yb1)＝(1-λ)*f(xb1，[yb1])+λ*f(xb1，[yb1]+1)
λ＝yb1-[yb1]
邻近点F1(xf1，yf1)的灰度值
f(xf1，yf1)＝(1-λ)*f(xf1，[yf1])+λ*f(xf1，[yf1]+1)
λ＝yf1-[yf1]
邻近点C1(xc1，yc1)的灰度值
f(xc1，yc1)＝(1-λ)*f(xc1，[yc1])+λ*f(xc1，[yc1]+1)
λ＝yc1-[yc1]
邻近点D1(xd1，yd1)的灰度值
f(xd1，yd1)＝(1-λ)*f(xd1，[yd1])+λ*f(xd1，[yd1]+1)
λ＝yd1-[yd1]
邻近点E1(xe1，ye1)的灰度值
f(xe1，ye1)＝(1-λ)*f(xe1，[ye1])+λ*f(xe1，[ye1]+1)
λ＝ye1-[ye1]
(2)若边缘像素点P(xp，yp)在第2区域内，
邻近点A2(xa2，ya2)的灰度值
f(xa2，ya2)＝(1-λ)*f([xa2]，ya2)+λ*f([xa2]+1，ya2)
λ＝xa2-[xa2]
邻近点B2(xb2，yb2)的灰度值
f(xb2，yb2)＝(1-λ)*f([xb2]，yb2)+λ*f([xb2]+1，yb2)
λ＝xb2-[xb2]
邻近点F2(xf2，yf2)的灰度值
f(xf2，yf2)＝(1-λ)*f([xf2]，yf2)+λ*f([xf2]+1，yf2)
λ＝xf2-[xf2]
邻近点C2(xc2，yc2)的灰度值
f(xc2，yc2)＝(1-λ)*f([xc2]，yc2)+λ*f([xc2]+1，yc2)
λ＝xc2-[xc2]
邻近点D2(xd2，yd2)的灰度值
f(xd2，yd2)＝(1-λ)*f([xd2]，yd2)+λ*f([xd2]+1，yd2)
λ＝xd2-[xd2]
邻近点E2(xe2，ye2)的灰度值
f(xe2，ye2)＝(1-λ)*f([xe2]，ye2)+λ*f([xe2]+1，ye2)
λ＝xe2-[xe2]
第3.4步：根据第3.3步所求边缘点沿梯度方向的邻近点的灰度值以及边缘点P(xp，yp)本身的灰度值，选取灰度值的前向差分和后向差分的平均值作为边缘点及其沿梯度方向的邻近点的对应的灰度差，其中，(xp，yp)为以像素为单位的边缘像素点P的坐标值，使用f(x，y)表示坐标为(x，y)的像素点的灰度值，所述的求取灰度差的方法如下：
(1)若边缘像素点P(xp，yp)在第1区域内，
邻近点A1(xa1，ya1)对应的灰度差
$f_{a1}=\frac{|f(x_p,y_p)-f(x_{a1},y_{a1})|}{2}+\frac{|f(x_{a1},y_{a1})-f(x_{b1},y_{b1})|}{2}$
邻近点B1(xb1，yb1)对应的灰度差
$f_{b1}=\frac{|f(x_{a1},y_{a1})-f(x_{b1},y_{b1})|}{2}+\frac{|f(x_{b1},y_{b1})-f(x_{f1},y_{f1})|}{2}$
边缘点P(xp，yp)对应的灰度差
$f_{p1}=\frac{|f(x_{c1},y_{c1})-f(x_p,y_p)|}{2}+\frac{|f(x_p,y_p)-f(x_{a1},y_{a1})|}{2}$
邻近点C1(xc1，yc1)对应的灰度差
$f_{c1}=\frac{|f(x_{d1},y_{d1})-f(x_{c1},y_{c1})|}{2}+\frac{|f(x_{c1},y_{c1})-f(x_p,y_p)|}{2}$
邻近点D1(xd，yd)对应的灰度差
$f_{d1}=\frac{|f(x_{e1},y_{e1})-f(x_{d1},y_{d1})|}{2}+\frac{|f(x_{d1},y_{d1})-f(x_{c1},y_{c1})|}{2}$
(2)若边缘像素点P(xp，yp)在第2区域内，
邻近点A2(xa2，ya2)对应的灰度差
$f_{a2}=\frac{|f(x_p,y_p)-f(x_{a2},y_{a2})|}{2}+\frac{|f(x_{a2},y_{a2})-f(x_{b2},y_{b2})|}{2}$
邻近点B2(xb2，yb2)对应的灰度差
$f_{b2}=\frac{|f(x_{a2},y_{a2})-f(x_{b2},y_{b2})|}{2}+\frac{|f(x_{b2},y_{b2})-f(x_{f2},y_{f2})|}{2}$
边缘点P(xp，yp)对应的灰度差
$f_{p2}=\frac{|f(x_{c2},y_{c2})-f(x_p,y_p)|}{2}+\frac{|f(x_p,y_p)-f(x_{a2},y_{a2})|}{2}$
邻近点C2(xc2，yc2)对应的灰度差
$f_{c2}=\frac{|f(x_{d2},y_{d2})-f(x_{c2},y_{c2})|}{2}+\frac{|f(x_{c2},y_{c2})-f(x_p,y_p)|}{2}$
邻近点D2(xd2，yd2)对应的灰度差
$f_{d2}=\frac{|f(x_{e2},y_{e2})-f(x_{d2},y_{d2})|}{2}+\frac{|f(x_{d2},y_{d2})-f(x_{c2},y_{c2})|}{2}$
第3.5步：根据第3.4步中所得边缘点及其沿梯度方向的邻近点的对应的灰度差，求取亚像素边缘点P′与像素边缘点P的距离差δ为

若边缘点及其沿梯度方向的邻近点对应的灰度差中有一个为0，则令δ＝0，并存储δ＝0出现的次数Count_zero，
第3.6步：求取亚像素边缘点P′的坐标，具体方法为：根据3.5步所求亚像素边缘点P′与像素级边缘点P的距离差δ以及像素级边缘点P的梯度方向直线的斜率kgradient，分别求出P点与P′点在x方向和y方向上的坐标差δx，δy，公式如下：
${\delta }_x=\frac{\delta }{\sqrt{k_{\mathrm{grad}\mathrm{ient}}^2+1}}$
${\delta }_y=\frac{k_{\mathrm{grad}\mathrm{ient}}*\delta }{\sqrt{k_{\mathrm{grad}\mathrm{ient}}^2+1}}$
其中，像素级边缘点P的梯度方向直线的斜率kgradient可由该像素边缘点的坐标(xp，yp)和粗定位圆心坐标(xoc，yoc)获得，即：
$k_{\mathrm{grad}\mathrm{ient}}=\frac{y_{\mathrm{oc}}-y_p}{x_{\mathrm{oc}}-x_p}$
那么，对于像素点P(xp，yp)，其对应的亚像素边缘点为P′(xp+δx，yp+δy)，
第3.7步：对所有边缘像素点P(xp，yp)，依据第3.2步～第3.6步依次求取对应的亚像素边缘点为P′(xp+δx，yp+δy)，其中，(xp+δx，yp+δy)为亚像素边缘点P′的坐标，
第4步：对亚像素边缘点进行滤波，针对亚像素边缘点出现的“孤立点”和噪声，分别用曲率滤波的方法和均值滤波的方法进行处理，
第4.1步：针对亚像素边缘出现的“孤立点”，用曲率滤波的方法进行处理，具体步骤如下：
(1)求取所有亚像素边缘点的曲率，方法如下：在顺序存储的亚像素边缘点中，为亚像素边缘点P′的坐标，为P′点的紧邻前一点的坐标，为P′点的紧邻后一点的坐标，那么亚像素边缘点P′的曲率为
$k=\frac{1}{r}=\frac{1}{\sqrt{({(x_0-x_{p_2^{\prime }})}^2+{(y_0-y_{p_2^{\prime }})}^2)}}$
其中，(x0，y0)为P′、P′点的紧邻前一点、P′点的紧邻后一点所组成圆的圆心，$x_0=\frac{a-b+c}{d},$$y_0=\frac{e-f+g}{-d},$
$a=(x_{p_1^{\prime }}+x_{p_2^{\prime }})(x_{p_2^{\prime }}-x_{p_1^{\prime }})(y_{p_3^{\prime }}-y_{p_2^{\prime }})$
$b=(x_{p_2^{\prime }}+x_{p_3^{\prime }})(x_{p_3^{\prime }}-x_{p_2^{\prime }})(y_{p_2^{\prime }}-y_{p_1^{\prime }})$
$c=(y_{p_1^{\prime }}-y_{p_3^{\prime }})(y_{p_2^{\prime }}-y_{p_1^{\prime }})(y_{p_3^{\prime }}-y_{p_2^{\prime }})$
$d=2[(x_{p_2^{\prime }}-x_{p_1^{\prime }})(y_{p_3^{\prime }}-y_{p_2^{\prime }})-(x_{p_3^{\prime }}-x_{p_2^{\prime }})(y_{p_2^{\prime }}-y_{p_1^{\prime }})]$
$e=(y_{p_1^{\prime }}+y_{p_2^{\prime }})(y_{p_2^{\prime }}-y_{p_1^{\prime }})(x_{p_3^{\prime }}-x_{p_2^{\prime }})$
$f=(y_{p_2^{\prime }}+y_{p_3^{\prime }})(y_{p_3^{\prime }}-y_{p_2^{\prime }})(x_{p_2^{\prime }}-x_{p_1^{\prime }})$
$g=(x_{p_1^{\prime }}-x_{p_3^{\prime }})(x_{p_2^{\prime }}-x_{p_1^{\prime }})(x_{p_3^{\prime }}-x_{p_2^{\prime }})$
(2)根据亚像素边缘点的曲率对所有亚像素边缘点进行滤波，滤波准则如下：
首先将各亚像素边缘点的曲率进行降序排列，取曲率阈值为该序列曲率值由高到低的第3*n个元素，其中，n为孤立点的个数，其数值为第3.5步中亚像素边缘定位时δ＝0的出现次数Count_zero；然后用曲率阈值对亚像素边缘点的曲率进行分割，如果亚像素边缘点P′的曲率大于曲率阈值，且P′的曲率大于其前后紧邻点的曲率，则认为点P′是孤立点，给予滤除，
第4.2步：将滤除“孤立点”后的所有亚像素边缘点用均值方法进行滤波，具体方法为：选取亚像素边缘点前后紧邻的两点的横、纵坐标与亚像素边缘点本身的横、纵坐标的均值作为该亚像素边缘点的新的横、纵坐标，
第5步：对滤波后的亚像素边缘点用最小二乘法拟合圆，最终得到圆形目标的圆心和半径。
2.根据权利要求1所述的圆形目标的圆心定位方法，其特征在于：
第1步所述对图像中圆形目标的圆心和半径进行粗定位的方法是：
第1.1步：对图像进行除噪、阈值分割，得到图像中各像素点灰度值为255或0的二值图像，
第1.2步：对二值图像中的圆形目标进行边界提取及边界跟踪，
第1.3步：对跟踪的边界点用周线形心法拟合圆，得到该圆形目标的粗定位圆心和粗定位圆半径，并存储粗定位圆心坐标(xoc，yoc)和粗定位圆半径Rc。

技术领域\n本发明涉及一种圆形目标的亚像素边缘定位以及亚像素边缘点滤波的方法，尤其涉及一种圆形目标的圆心定位方法。\n背景技术\n在众多图像分析和处理系统中，快速而准确地获取圆心被广泛地应用在目标识别以及定位上，如在三维重构的系统中，高精度地获取标定板圆心对整个三维重构精度有至关重要的影响。目前求取圆心的常见算法有：形心法，Hough变换法，高斯曲面拟合法，最小二乘曲线拟合法。形心法要求图像灰度分布比较均匀，否则会产生较大误差；Hough变换法需要对每一个边界点逐点投票、记录，计算时间较长，而且占用计算机内存较大，在实际应用中受到限制；高斯曲面拟合是利用圆形光斑的灰度分布近似于高斯模型的原理进行曲面拟合，其定位精度高，但是如果处理圆面积较大，则运算量很大，计算时间长，不适用于实时光学检测系统中；最小二乘曲线拟合法是最为常见的方法，因为其实现简单，计算量小，在实际中得到广泛应用，但是目前最小二乘曲线拟合法定位圆中往往都是利用求取的像素级边缘点进行最小二乘法拟合圆心，使得其拟合精度受到了一定的限制，所以能否获取精确的图像中圆形的亚像素边缘，对于提高圆心的拟合精度有着重要的影响。\n目前常见的亚像素边缘定位算法有：参数拟合法，矩的方法，插值法以及利用图像其他几何特征的方法。文献“Subpixel edge localization and the interpolation of still images”(Jensen K，Anastassiou，IEEE Trans on PAMI，1995，4(3)：285～295)是插值法，是将像素级边缘定位算法用于插值放大的图像，在实际图像灰度分布规律与插值规律相一致的前提下，使边缘定位达到亚像素精度；文献“利用切线方向信息检测亚像素边缘”(章毓晋，傅倬，模式识别与人工智能，1997，10(1)：69～74)主要利用了切线方向的信息，对已知形状的目标进行亚像素边缘定位，这种方法过分依赖于粗定位的圆心和半径信息，没有充分利用图像的灰度值信息；文献“A fast subpixel edge detection method usingSobel-Zernike moments operator”(Qu Ying-Dong，Cui Cheng-Song，Chen San-Ben，LiJin-Quan，Image and Vision Computing，2005，23：11～17)中用矩的方法实现边缘亚像素定位。基于矩的方法目前应用比较广泛，但是由于要用到多个7×7模板，因此不适用于大型位图的实时处理；文献“图像的快速亚像素边缘检测方法”(刘力双，张铫，光电子激光，2005，16(8)：993～996)介绍了亚像素边缘参数拟合的原理，其采用的是Sobel算子和参数拟合的方法实现亚像素边缘定位，其中存在两个问题，第一，Sobel算子对绝大多数图像都不能实现整像素级边缘定位，导致亚像素定位不准；第二，该文献提供的方法不能确保二次曲线型的边缘定位精度。\n发明内容\n针对现有技术所存在的缺点和限制，本发明的目的在于提供一种能够提高测量系统精度的圆形目标的圆心定位方法。\n本发明设计一种首先将图像中的圆形目标进行圆心粗定位，在用canny算子获取圆形目标的像素级边缘点后，综合考虑圆的几何特征及边缘灰度分布特征求取亚像素边缘点坐标，然后对精确的亚像素边缘点进行滤波，最后用最小二乘法拟合滤波后的亚像素边缘点求得圆形目标的圆心的方法。本发明采用如下技术方案：\n一种圆形目标的圆心定位方法，其主要步骤为：\n第1步：对图像中圆形目标的圆心和半径进行粗定位，得到圆形目标的粗定位圆心(xoc，yoc)和粗定位圆半径Rc，其中，(xoc，yoc)为图像坐标系下圆心的坐标，具体方法如下：\n第1.1步：对图像进行除噪、阈值分割，得到图像中各像素点灰度值为255或0的二值图像；\n第1.2步：对二值图像中的圆形目标进行边界提取及边界跟踪；\n第1.3步：对跟踪的边界点用周线形心法拟合圆，得到该圆形目标的粗定位圆心和粗定位圆半径，并存储粗定位圆心坐标(xoc，yoc)和粗定位圆半径Rc；\n第2步：对图像中的圆形目标进行像素级边缘定位，具体步骤如下：\n第2.1步：根据粗定位圆心和粗定位圆半径，从图像中提取一个正方形区域，称为“感兴趣区域”，该正方形区域的提取方法为：以图像中圆形目标的粗定位圆心作为正方形区域的中心点，以圆形目标的粗定位直径加上2～6个像素点作为正方形区域的边长；\n第2.2步：用canny算子对图像中所提取的感兴趣区域进行像素级边缘检测，得到圆形目标的像素级边缘，再对canny算子提取的圆形目标的像素级边缘点进行边界跟踪，得到每个像素级边缘点的坐标，且按顺时针顺序存储每个像素级边缘点的坐标；\n第3步：对图像中的圆形目标进行亚像素边缘定位，具体步骤如下：\n第3.1步：将圆形目标的边缘像素点进行区域划分，划分方法如下：称以圆形目标的粗定位圆心为原点、以水平向右为正方向、以像素为单位长度的直线为x轴；称以圆形目标的粗定位圆心为原点、以垂直向上为正方向、以像素为单位长度的直线为y轴；以粗定位圆心为旋转点，由x轴正方向逆时针旋转的角度为θ，根据θ的大小将圆形目标的边缘像素点划分为两个部分，一部分是θ∈[0°，45°]∪[135°，225°]∪[315°，360°]所对应的边缘像素点，称为第1区域的边缘像素点，另一部分是θ∈(45°，135°)∪(225°，315°)对应的边缘像素点，称为第2区域的边缘像素点，粗定位圆心与边缘像素点的连线为该边缘像素点的梯度方向；\n第3.2步：求取边缘像素点沿梯度方向的邻近点，(xp，yp)为以像素为单位的边缘点P的坐标值，具体方法如下：若边缘像素点P(xp，yp)在第1区域内，取该像素点的梯度方向直线与直线x＝xp+1的交点A1(xa1，ya1)为该边缘像素点沿梯度方向的右边第一邻近点；取该像素点的梯度方向直线与直线x＝xp+2的交点B1(xb1，yb1)为该边缘像素点沿梯度方向的右边第二邻近点；取该像素点的梯度方向直线与直线x＝xp+3的交点F(xf，yf)为该边缘像素点沿梯度方向的右边第三邻近点；取该像素点的梯度方向直线与直线x＝xp-1的交点C1(xc1，yc1)为该边缘像素点沿梯度方向的左边第一邻近点；取该像素点的梯度方向直线与直线x＝xp-2的交点D1(xd1，yd1)为该边缘像素点沿梯度方向的左边第二邻近点，取该像素点的梯度方向直线与直线x＝xp-3的交点E1(xe1，ye1)为该边缘像素点沿梯度方向的左边第三邻近点；若边缘像素点P(xp，yp)在第2区域内，取该像素点的梯度方向直线与直线y＝yp+1的交点A2(xa2，ya2)为该边缘像素点沿梯度方向的上边第一紧邻点；取该像素点的梯度方向直线与直线y＝yp+2的交点B2(xb2，yb2)为该边缘像素点沿梯度方向的上边第二紧邻点；取该像素点的梯度方向直线与直线y＝yp+3的交点F2(xf2，yf2)为该边缘像素点沿梯度方向的上边第三紧邻点；取该像素点的梯度方向直线与直线y＝yp-1的交点C2(xc2，yc2)为该边缘像素点沿梯度方向的下边第一紧邻点；取该像素点的梯度方向直线与直线y＝yp-2的交点D2(xd2，yd2)为该边缘像素点沿梯度方向的上边第二紧邻点；取该像素点的梯度方向直线与直线y＝yp-3的交点E2(xe2，ye2)为该边缘像素点沿梯度方向的上边第三紧邻点；\n第3.3步：采用线性灰度插值的方法求取以像素为单位的边缘像素点P(xp，yp)沿梯度方向的邻近点的灰度值，其中，(xp，yp)为以像素为单位的边缘像素点P的坐标值，使用f(x，y)表示坐标为(x，y)的像素点的灰度值，使用符号[]表示取整数部分，所述的用于获取邻近点灰度值的线性灰度插值方法如下：\n(1)若边缘像素点P(xp，yp)在第1区域内，\n邻近点A1(xa1，ya1)的灰度值\nf(xa1，ya1)＝(1-λ)*f(xa1，[ya1])+λ*f(xa1，[ya1]+1)\nλ＝ya1-[ya1]，\n邻近点B1(xb1，yb1)的灰度值\nf(xb1，yb1)＝(1-λ)*f(xb1，[yb1])+λ*f(xb1，[yb1]+1)\nλ＝yb1-[yb1]\n邻近点F1(xf1，yf1)的灰度值\nf(xf1，yf1)＝(1-λ)*f(xf1，[yf1])+λ*f(xf1，[yf1]+1)\nλ＝yf1-[yf1]\n邻近点C1(xc1，yc1)的灰度值\nf(xc1，yc1)＝(1-λ)*f(xc1，[yc1])+λ*f(xc1，[yc1]+1)\nλ＝yc1-[yc1]\n邻近点D1(xd，yd)的灰度值\nf(xd1，yd1)＝(1-λ)*f(xd1，[yd1])+λ*f(xd1，[yd1]+1)\nλ＝yd1-[yd1]\n邻近点E1(xe1，ye1)的灰度值\nf(xe1，ye1)＝(1-λ)*f(xe1，[ye1])+λ*f(xe1，[ye1]+1)\nλ＝ye1-[ye1]\n(2)若边缘像素点P(xp，yp)在第2区域内，\n邻近点A2(xa2，ya2)的灰度值\nf(xa2，ya2)＝(1-λ)*f([xa2]，ya2)+λ*f([xa2]+1，ya2)\nλ＝xa2-[xa2]\n邻近点B2(xb2，yb2)的灰度值\nf(xb2，yb2)＝(1-λ)*f([xb2]，yb2)+λ*f([xb2]+1，yb2)\nλ＝xb2-[xb2]\n邻近点F2(xf2，yf2)的灰度值\nf(xf2，yf2)＝(1-λ)*f([xf2]，yf2)+λ*f([xf2]+1，yf2)\nλ＝xf2-[xf2]\n邻近点C2(xc2，yc2)的灰度值\nf(xc2，yc2)＝(1-λ)*f([xc2]，yc2)+λ*f([xc2]+1，yc2)\nλ＝xc2-[xc2]\n邻近点D2(xd2，yd2)的灰度值\nf(xd2，yd2)＝(1-λ)*f([xd2]，yd2)+λ*f([xd2]+1，yd2)\nλ＝xd2-[xd2]\n邻近点E2(xe2，ye2)的灰度值\nf(xe2，ye2)＝(1-λ)*f([xe2]，ye2)+λ*f([xe2]+1，ye2)\nλ＝xe2-[xe2]\n第3.4步：根据第3.3步所求边缘点沿梯度方向的邻近点的灰度值以及边缘点P(xp，yp)本身的灰度值，选取灰度值的前向差分和后向差分的平均值作为边缘点及其沿梯度方向的邻近点的对应的灰度差，其中，(xp，yp)为以像素为单位的边缘像素点P的坐标值，使用f(x，y)表示坐标为(x，y)的像素点的灰度值，所述的求取灰度差的方法如下：\n(1)若边缘像素点P(xp，yp)在第1区域内，\n邻近点A1(xa1，ya1)对应的灰度差\n$f_{a1}=\frac{|f(x_p,y_p)-f(x_{a1},y_{a1})|}{2}+\frac{|f(x_{a1},y_{a1})-f(x_{b1},y_{b1})|}{2}$\n邻近点B1(xb1，yb1)对应的灰度差\n$f_{b1}=\frac{|f(x_{a1},y_{a1})-f(x_{b1},y_{b1})|}{2}+\frac{|f(x_{b1},y_{b1})-f(x_{f1},y_{f1})|}{2}$\n边缘点P(xp，yp)对应的灰度差\n$f_{p1}=\frac{|f(x_{c1},y_{c1})-f(x_p,y_p)|}{2}+\frac{|f(x_p,y_p)-f(x_{a1},y_{a1})|}{2}$\n邻近点C1(xc1，yc1)对应的灰度差\n$f_{c1}=\frac{|f(x_{d1},y_{d1})-f(x_{c1},y_{c1})|}{2}+\frac{|f(x_{c1},y_{c1})-f(x_p,y_p)|}{2}$\n邻近点D1(xd，yd)对应的灰度差\n$f_{d1}=\frac{|f(x_{e1},y_{e1})-f(x_{d1},y_{d1})|}{2}+\frac{|f(x_{d1},y_{d1})-f(x_{c1},y_{c1})|}{2}$\n(2)若边缘像素点P(xp，yp)在第2区域内，\n邻近点A2(xa2，ya2)对应的灰度差\n$f_{a2}=\frac{|f(x_p,y_p)-f(x_{a2},y_{a2})|}{2}+\frac{|f(x_{a2},y_{a2})-f(x_{b2},y_{b2})|}{2}$\n邻近点B2(xb2，yb2)对应的灰度差\n$f_{b2}=\frac{|f(x_{a2},y_{a2})-f(x_{b2},y_{b2})|}{2}+\frac{|f(x_{b2},y_{b2})-f(x_{f2},y_{f2})|}{2}$\n边缘点P(xp，yp)对应的灰度差\n$f_{p2}=\frac{|f(x_{c2},y_{c2})-f(x_p,y_p)|}{2}+\frac{|f(x_p,y_p)-f(x_{a2},y_{a2})|}{2}$\n邻近点C2(xc2，yc2)对应的灰度差\n$f_{c2}=\frac{|f(x_{d2},y_{d2})-f(x_{c2},y_{c2})|}{2}+\frac{|f(x_{c2},y_{c2})-f(x_p,y_p)|}{2}$\n邻近点D2(xd2，yd2)对应的灰度差\n$f_{d2}=\frac{|f(x_{e2},y_{e2})-f(x_{d2},y_{d2})|}{2}+\frac{|f(x_{d2},y_{d2})-f(x_{c2},y_{c2})|}{2}$\n第3.5步：根据第3.4步中所得边缘点及其沿梯度方向的邻近点的对应的灰度差，求取亚像素边缘点P′与像素边缘点P的距离差δ为\n\n若边缘点及其沿梯度方向的邻近点对应的灰度差中有一个为0，则令δ＝0，并存储δ＝0出现的次数Count_zero；\n第3.6步：求取亚像素边缘点P′的坐标，具体方法为：根据3.5步所求亚像素边缘点P′与像素级边缘点P的距离差δ以及像素级边缘点P的梯度方向直线的斜率kgradient，分别求出P点与P′点在x方向和y方向上的坐标差δx，δy，公式如下：\n${\delta }_x=\frac{\delta }{\sqrt{k_{\mathrm{grad}\mathrm{ient}}^2+1}}$\n${\delta }_y=\frac{k_{\mathrm{grad}\mathrm{ient}}*\delta }{\sqrt{k_{\mathrm{grad}\mathrm{ient}}^2+1}}$\n其中，像素级边缘点P的梯度方向直线的斜率kgradient可由该像素边缘点的坐标(xp，yp)和粗定位圆心坐标(xoc，yoc)获得，即：\n$k_{\mathrm{grad}\mathrm{ient}}=\frac{y_{\mathrm{oc}}-y_p}{x_{\mathrm{oc}}-x_p}$\n那么，对于像素点P(xp，yp)，其对应的亚像素边缘点为P′(xp+δx，yp+δy)：\n第3.7步：对所有边缘像素点P(xp，yp)，依据第3.2步～第3.6步依次求取对应的亚像素边缘点为P′(xp+δx，yp+δy)，其中，(xp+δx，yp+δy)为亚像素边缘点P′的坐标；\n第4步：对亚像素边缘点进行滤波，针对亚像素边缘点出现的“孤立点”和噪声，分别用曲率滤波的方法和均值滤波的方法进行处理：\n第4.1步：针对亚像素边缘出现的“孤立点”，用曲率滤波的方法进行处理，具体步骤如下：\n(1)求取所有亚像素边缘点的曲率，方法如下：在顺序存储的亚像素边缘点中，为亚像素边缘点P′的坐标，为P′点的紧邻前一点的坐标，为P′点的紧邻后一点的坐标，那么亚像素边缘点P′的曲率为\n$k=\frac{1}{r}=\frac{1}{\sqrt{({(x_0-x_{p_2^{\prime }})}^2+{(y_0-y_{p_2^{\prime }})}^2)}}$\n其中，(x0，y0)为P′、P′点的紧邻前一点、P′点的紧邻后一点所组成圆的圆心，\n$x_0=\frac{a-b+c}{d},$ $y_0=\frac{e-f+g}{-d},$\n$a=(x_{p_1^{\prime }}+x_{p_2^{\prime }})(x_{p_2^{\prime }}-x_{p_1^{\prime }})(y_{p_3^{\prime }}-y_{p_2^{\prime }})$\n$b=(x_{p_2^{\prime }}+x_{p_3^{\prime }})(x_{p_3^{\prime }}-x_{p_2^{\prime }})(y_{p_2^{\prime }}-y_{p_1^{\prime }})$\n$c=(y_{p_1^{\prime }}-y_{p_3^{\prime }})(y_{p_2^{\prime }}-y_{p_1^{\prime }})(y_{p_3^{\prime }}-y_{p_2^{\prime }})$\n$d=2[(x_{p_2^{\prime }}-x_{p_1^{\prime }})(y_{p_3^{\prime }}-y_{p_2^{\prime }})-(x_{p_3^{\prime }}-x_{p_2^{\prime }})(y_{p_2^{\prime }}-y_{p_1^{\prime }})]$\n$e=(y_{p_1^{\prime }}+y_{p_2^{\prime }})(y_{p_2^{\prime }}-y_{p_1^{\prime }})(x_{p_3^{\prime }}-x_{p_2^{\prime }})$\n$f=(y_{p_2^{\prime }}+y_{p_3^{\prime }})(y_{p_3^{\prime }}-y_{p_2^{\prime }})(x_{p_2^{\prime }}-x_{p_1^{\prime }})$\n$g=(x_{p_1^{\prime }}-x_{p_3^{\prime }})(x_{p_2^{\prime }}-x_{p_1^{\prime }})(x_{p_3^{\prime }}-x_{p_2^{\prime }})$\n(2)根据亚像素边缘点的曲率对所有亚像素边缘点进行滤波，滤波准则如下：\n首先将各亚像素边缘点的曲率进行降序排列，取曲率阈值为该序列曲率值由高到低的第3*n个元素，其中，n为孤立点的个数，其数值为第3.5步中亚像素边缘定位时δ＝0的出现次数Count_zero；然后用曲率阈值对亚像素边缘点的曲率进行分割，如果亚像素边缘点P的曲率大于曲率阈值，且P的曲率大于其前后紧邻点的曲率，则认为点P是孤立点，给予滤除；\n第4.2步：将滤除“孤立点”后的所有亚像素边缘点用均值方法进行滤波，具体方法为：选取亚像素边缘点前后紧邻的两点的横、纵坐标与亚像素边缘点本身的横、纵坐标的均值作为该亚像素边缘点的新的横、纵坐标；\n第5步：对滤波后的亚像素边缘点用最小二乘法拟合圆，最终得到圆形目标的圆心和半径。\n与现有技术相比，本发明具有如下优点：\n(1)与传统的二维图像边缘定位方法相比，本发明提出的亚像素边缘定位的方法基于粗定位圆心所确定的感兴趣区域进行，而且在拟合边缘灰度差时将高斯曲面拟合简化为沿梯度方向的高斯曲线拟合，简化了计算，提高算法的快速性。\n(2)本发明提出的亚像素边缘定位的方法在考虑每个边缘像素点的梯度方向时，不拘泥于传统的边缘检测算子确定的几个特殊的方向，而是根据圆的特征考虑了每个边缘像素点的具体梯度方向，并且利用了线性灰度插值获取边缘像素点沿梯度方向的邻近点的灰度值，进而使边缘点的定位更加准确。\n(3)本发明针对亚像素边缘点中出现的“孤立点”和噪声，提出了一种简单实用的滤波方法，该方法不仅原理和算法简单，而且滤波效果很好。\n(4)本发明提出的圆形目标的圆心定位方法与已有的方法相比，有更高的精度和更好的稳定性，并且具有很强的实用性。\n附图说明\n图1是圆形目标的圆心定位方法具体步骤的流程图。\n图2是感兴趣区域提取图。\n图3是非极大值抑制示意图，其中图3(a)为根据梯度方向划分的四个扇区图；3(b)为3*3邻域图。\n图4(a)是边缘像素点P在第1区域时的定义像素点沿梯度方向邻近点的原理图；图4(b)是边缘像素点P在第2区域时的定义像素点沿梯度方向邻近点的原理图。\n图5(a)是当边缘像素点P在第1区域时，边缘像素点P沿梯度方向的邻近点定义图；图5(b)是当边缘像素点P在第2区域时，边缘像素点P沿梯度方向的邻近点定义图。\n图6是亚像素边缘点求取原理图。\n具体实施方式\n下面结合附图对本发明的具体实施方式作进一步描述。根据上述方法，在Windows操作系统中通过VC++6.0平台用C++编程实现了圆形目标的圆心定位操作。\n应用本方法进行圆形目标的圆心定位主要包括圆心粗定位、像素级边缘定位、亚像素边缘定位、亚像素边缘点滤波、最小二乘拟合求得圆心五个操作步骤，具体步骤的流程图如图1所示，应用本方法进行圆形目标的圆心定位的具体步骤如下(步骤中未标明的单位都以像素为单位)：\n第1步：对图像中的圆形目标的圆心和半径进行粗定位，得到圆形目标的粗定位圆心(xoc，yoc)和粗定位圆半径Rc，其中，(xoc，yoc)为图像坐标系下圆心的坐标，具体步骤如下：\n第1.1步：对图像进行除噪、阈值分割，得到图像中各像素点灰度值为255或0的二值图像。其中，阈值分割采用迭代阈值法。迭代阈值法能自动计算出比较合适的分割阈值，其计算步骤如下：\n(1)选择阈值T，通常选择图像的平均灰度值作为初始阈值；\n(2)通过初始阈值T，将图像的平均灰度值分成为两组R1和R2；\n(3)计算两组平均灰度值μ1和μ2；\n(4)重新选择阈值T，新的T定义为：T＝(μ1+μ2)/2；\n循环(2)到(4)，直到两组的平均灰度值μ1和μ2不再发生改变，因此得到所需要的阈值，再根据所得到的阈值对图像进行二值化。\n第1.2步：对二值图像中的圆形目标进行边界提取及边界跟踪。在二值图像中，灰度值为255的像素即为圆形目标，对圆形目标进行边界提取的方法如下：如果原图中有一点为白，而相邻的8个像素全部为黑，则将该像素点的灰度值设为0。经过边界提取后，对该图中的目标轮廓线进行跟踪并存储，可得到圆形目标的边界点及其坐标。\n第1.3步：对跟踪的边界点用周线形心法拟合圆，得到该圆形目标的粗定位圆心和粗定位圆半径，并存储粗定位圆心坐标(xoc，yoc)和粗定位圆半径Rc。周线形心法是一种最为简单的亚像素定位算法，算法如下：设跟踪的边界点坐标为(i，j)，跟踪边界点总数为N，根据周线形心法，圆形目标的圆心坐标为：\n$\left\{\begin{array}{c}{x}_{\mathrm{oc}}=\frac{\mathrm{\Sigma i}}{N}\\ y_{\mathrm{oc}}=\frac{\mathrm{\Sigma j}}{N}\end{array}\right.$\n圆形目标的半径为各边界点到圆心距离的平均值，即：\n$R_c=\frac{\Sigma \sqrt{{(i-x_{\mathrm{oc}})}^2+{(j-y_{\mathrm{oc}})}^2}}{N}$\n第2步：对图像中的圆形目标进行像素级边缘定位，具体步骤如下：\n第2.1步：根据粗定位圆心和粗定位圆半径，从图像中提取一个正方形区域，称为“感兴趣区域”，该正方形区域的提取方法为：以图像中圆形目标的粗定位圆心作为正方形区域的中心点，以圆形目标的粗定位直径加上2～6个像素点作为正方形区域的边长。如图2所示，图中圆代表图像中的圆形目标，正方形区域即为感兴趣区域，其中，O为粗定位圆心，Rc为粗定位半径，d一般取2～6个像素。后续步骤只针对感兴趣区域进行处理。\n第2.2步：用canny算子对图像中所提取的感兴趣区域进行像素级边缘检测，得到圆形目标的像素级边缘，再对canny算子提取的圆形目标的像素级边缘点进行边界跟踪，得到每个像素级边缘点的坐标，且按顺时针顺序存储每个像素级边缘点的坐标。其中，Canny边缘检测对信噪比与定位乘积进行测度，属于先平滑后再求导的方法，其主要步骤和基本原理如下：\n(1)用高斯滤波器平滑图像，高斯平滑函数如下：\n$H(x,y)=e^{-\frac{x^2+y^2}{2{\sigma }^2}}$\n其中，σ为高斯函数的标准偏差，x，y分别是图像中像素点的横、纵坐标。\n对图像进行高斯滤波，即对图像fI(x，y)用高斯函数H(x，y)进行卷积，得到滤波后的平滑图像G(x，y)：\nG(x，y)＝fI(x，y)*H(x，y)\n(2)用一阶偏导数的有限差分来计算梯度的幅值和方向。一阶差分卷积模板如下：\n$H_1=\left|\begin{array}{cc}-1& -1\\ 1& 1\end{array}\right|$ $H_2=\left|\begin{array}{cc}1& -1\\ 1& -1\end{array}\right|$\n将平滑后的图像G(x，y)用以上两个模板进行卷积，得到水平和垂直方向上的梯度图像分别为进而由可得到梯度的幅值和方向推算式如下：\n\n\n(3)对梯度幅值进行非极大值抑制(non-maxima suppression，NMS)：图3是非极大值抑制示意图。根据梯度的方向，四个扇区的标号为0到3，如图3(a)所示，对应3*3邻域的四种可能组合，3*3邻域如图3(b)所示。在每一个像素点上，将邻域的中心象素Pcp与沿着梯度线的两个像素相比，如果Pcp的梯度值不比沿梯度线的两个相邻像素梯度值大，则该点不为边缘点。\n(4)用双阈值算法检测和连接边缘。双阈值算法对非极大值抑制图像作用两个阈值τ1和τ2，且2τ1≈τ2，从而可以得到两个阈值边缘图像N1[i，j]和N2[i，j]。由于N2[i，j]使用高阈值得到，因而含有很少的假边缘，但有间断。双阈值法要在N2[i，j]中把边缘连接成轮廓，当到达轮廓的端点时，该算法就在N1[i，j]的8邻域点中寻找可以连接到轮廓上的边缘，这样，算法不断地在N1[i，j]中收集边缘，直到将N2[i，j]连接起来为止。\n第3步：对图像中的圆形目标进行亚像素边缘定位，具体步骤如下：\n第3.1步：将圆形目标的边缘像素点进行区域划分，划分方法如下：称以圆形目标的粗定位圆心为原点、以水平向右为正方向、以像素为单位长度的直线为x轴；称以圆形目标的粗定位圆心为原点、以垂直向上为正方向、以像素为单位长度的直线为y轴；以粗定位圆心为旋转点，由x轴正方向逆时针旋转的角度为θ，根据θ的大小将圆形目标的边缘像素点划分为两个部分，一部分是θ∈[0°，45°]∪[135°，225°]∪[315°，360°]所对应的边缘像素点，称为第1区域的边缘像素点，另一部分是θ∈(45°，135°)∪(225°，315°)对应的边缘像素点，称为第2区域的边缘像素点，粗定位圆心与边缘像素点的连线为该边缘像素点的梯度方向。\n第3.2步：求取边缘像素点沿梯度方向的邻近点，具体原理如下：若边缘像素点P(xp，yp)在第1区域内，图4(a)是定义像素点沿梯度方向邻近点的原理图，(xp，yp)为以像素为单位的边缘点P的坐标值，LG为P点的梯度方向直线，直线y1、y2、y3、y4为P点附近以像素点为单位的四条横线，分别与直线LG相交于点B′1、A′1、C′1、D′1，直线x1，x2，x3，x4为P点附近以像素点为单位的四条纵线，分别与直线LG相交于点D1、C1、A1、B1，选取{B′1、A′1、C′1、D′1}和{B1、A1、C1、D1}两组点作为边缘点沿梯度方向的候选邻近点，由于P点在第1区域，梯度方向以x方向为主，即LG斜率的绝对值小于1，因此，候选邻近点{B′1、A′1、C′1、D′1}到P点的距离明显要大于候选邻近点{B1、A1、C1、D1}到P点的距离，若选用{B′1、A′1、C′1、D′1}作为P点沿梯度方向的邻近点，会导致边缘定位不够精确，因此，选用{B1、A1、C1、D1}作为其邻近点；若边缘像素点P(xp，yp)在第2区域内，图4(b)是定义像素点沿梯度方向邻近点的原理图，(xp，yp)为以像素为单位的边缘点P的坐标值，LG为P点的梯度方向直线，直线y1、y2、y3、y4为P点附近以像素点为单位的四条横线，分别与直线LG相交于点B2、A2、C2、D2，直线x1，x2，x3，x4为P点附近以像素点为单位的四条纵线，分别与直线LG相交于点D′2、C′2、A′2、B′2，选取{B′2、A′2、C′2、D′2}和{B2、A2、C2、D2}两组点作为边缘点沿梯度方向的候选邻近点。由于P点在第2区域，梯度方向以y方向为主，即LG的斜率的绝对值大于1，因此，候选邻近点{B′2、A′2、C′2、D′2}到P点的距离明显要大于候选邻近点{B2、A2、C2、D2}到P点的距离，若选用{B′2、A′2、C′2、D′2}作为P点沿梯度方向的邻近点，会导致边缘定位不够精确，因此，选用{B2、A2、C2、D2}作为其邻近点。因此，针对第3.1步划分的两个不同区域的像素点，其沿梯度方向的邻近点的定义如下：(xp，yp)为以像素为单位的边缘点P的坐标值，若边缘像素点P(xp，yp)在第1区域内，图5(a)是像素点沿梯度方向的邻近点定义图，LG为边缘像素点P(xp，yp)的梯度方向直线，L1代表直线x＝xp-3，L2代表直线x＝xp-2，L3代表直线x＝xp-1，L4代表直线x＝xp+1，L5代表直线x＝xp+2，L6代表直线x＝xp+3。取该像素点的梯度方向直线与直线x＝xp+1的交点A1(xa1，ya1)为该边缘像素点沿梯度方向的右边第一邻近点；取该像素点的梯度方向直线与直线x＝xp+2的交点B1(xb1，yb1)为该边缘像素点沿梯度方向的右边第二邻近点；取该像素点的梯度方向直线与直线x＝xp+3的交点F(xf，yf)为该边缘像素点沿梯度方向的右边第三邻近点；取该像素点的梯度方向直线与直线x＝xp-1的交点C1(xc1，yc1)为该边缘像素点沿梯度方向的左边第一邻近点；取该像素点的梯度方向直线与直线x＝xp-2的交点D1(xd1，yd1)为该边缘像素点沿梯度方向的左边第二邻近点，取该像素点的梯度方向直线与直线x＝xp-3的交点E1(xe1，ye1)为该边缘像素点沿梯度方向的左边第三邻近点；若边缘像素点P(xp，yp)在第2区域内，图5(b)是像素点沿梯度方向的邻近点定义图，LG为边缘像素点P(xp，yp)的梯度方向直线，取该像素点的梯度方向直线与直线y＝yp+1的交点A2(xa2，ya2)为该边缘像素点沿梯度方向的上边第一紧邻点；取该像素点的梯度方向直线与直线y＝yp+2的交点B2(xb2，yb2)为该边缘像素点沿梯度方向的上边第二紧邻点；取该像素点的梯度万向直线与直线y＝yp+3的交点F2(xf2，yf2)为该边缘像素点沿梯度方向的上边第三紧邻点；取该像素点的梯度方向直线与直线y＝yp-1的交点C2(xc2，yc2)为该边缘像素点沿梯度方向的下边第一紧邻点；取该像素点的梯度方向直线与直线y＝yp-2的交点D2(xd2，yd2)为该边缘像素点沿梯度方向的上边第二紧邻点；取该像素点的梯度方向直线与直线y＝yp-3的交点E2(xe2，ye2)为该边缘像素点沿梯度方向的上边第三紧邻点。\n第3.3步：采用线性灰度插值的方法求取以像素为单位的边缘像素点P(xp，yp)沿梯度方向的邻近点的灰度值，具体方法为：首先找出与该像素点的沿梯度方向的邻近点距离最近的两个整像素点，然后对这两个整像素点进行线性灰度插值，插值结果为该边缘像素点的沿梯度方向的邻近点的灰度值，其中，(xp，yp)为以像素为单位的边缘像素点P的坐标值，使用f(x，y)表示坐标为(x，y)的像素点的灰度值，使用符号[]表示取整数部分，所述的用于获取邻近点灰度值的线性灰度插值方法如下：\n(1)若边缘像素点P(xp，yp)在第1区域内，\n邻近点A1(xa1，ya1)的灰度值\nf(xa1，ya1)＝(1-λ)*f(xa1，[ya1])+λ*f(xa1，[ya1]+1)\nλ＝ya1－[ya1]，\n邻近点B1(xb1，yb1)的灰度值\nf(xb1，yb1)＝(1-λ)*f(xb1，[yb1])+λ*f(xb1，[yb1]+1)\nλ＝yb1-[yb1]\n邻近点F1(xf1，yf1)的灰度值\nf(xf1，yf1)＝(1-λ)*f(xf1，[yf1])+λ*f(xf1，[yf1]+1)\nλ＝yf1-[yf1]\n邻近点C1(xc1，yc1)的灰度值\nf(xc1，yc1)＝(1-λ)*f(xc1，[yc1])+λ*f(xc1，[yc1]+1)\nλ＝yc1-[yc1]\n邻近点D1(xd，yd)的灰度值\nf(xd1，yd1)＝(1-λ)*f(xd1，[yd1])+λ*f(xd1，[yd1]+1)\nλ＝yd1-[yd1]\n邻近点E1(xe1，ye1)的灰度值\nf(xe1，ye1)＝(1-λ)*f(xe1，[ye1])+λ*f(xe1，[ye1]+1)\nλ＝ye1-[ye1]\n(2)若边缘像素点P(xp，yp)在第2区域内，\n邻近点A2(xa2，ya2)的灰度值\nf(xa2，ya2)＝(1-λ)*f([xa2]，ya2)+λ*f([xa2]+1，ya2)\nλ＝xa2-[xa2]\n邻近点B2(xb2，yb2)的灰度值\nf(xb2，yb2)＝(1-λ)*f([xb2]，yb2)+λ*f([xb2]+1，yb2)\nλ＝xb2-[xb2]\n邻近点F2(xf2，yf2)的灰度值\nf(xf2，yf2)＝(1-λ)*f([xf2]，yf2)+λ*f([xf2]+1，yf2)\nλ＝xf2-[xf2]\n邻近点C2(xc2，yc2)的灰度值\nf(xc2，yc2)＝(1-λ)*f([xc2]，yc2)+λ*f([xc2]+1，yc2)\nλ＝xc2-[xc2]\n邻近点D2(xd2，yd2)的灰度值\nf(xd2，yd2)＝(1-λ)*f([xd2]，yd2)+λ*f([xd2]+1，yd2)\nλ＝xd2-[xd2]\n邻近点E2(xe2，ye2)的灰度值\nf(xe2，ye2)＝(1-λ)*f([xe2]，ye2)+λ*f([xe2]+1，ye2)\nλ＝xe2-[xe2]\n第3.4步：根据第3.3步所求边缘点沿梯度方向的邻近点的灰度值以及边缘点P(xp，yp)本身的灰度值，选取灰度值的前向差分和后向差分的平均值作为边缘点及其沿梯度方向的邻近点的对应的灰度差，其中，(xp，yp)为以像素为单位的边缘像素点P的坐标值，使用f(x，y)表示坐标为(x，y)的像素点的灰度值，所述的求取灰度差的方法如下：\n(1)若边缘像素点P(xp，yp)在第1区域内，\n邻近点A1(xa1，ya1)对应的灰度差\n$f_{a1}=\frac{|f(x_p,y_p)-f(x_{a1},y_{a1})|}{2}+\frac{|f(x_{a1},y_{a1})-f(x_{b1},y_{b1})|}{2}$\n邻近点B1(xb1，yb1)对应的灰度差\n$f_{b1}=\frac{|f(x_{a1},y_{a1})-f(x_{b1},y_{b1})|}{2}+\frac{|f(x_{b1},y_{b1})-f(x_{f1},y_{f1})|}{2}$\n边缘点P(xp，yp)对应的灰度差\n$f_{p1}=\frac{|f(x_{c1},y_{c1})-f(x_p,y_p)|}{2}+\frac{|f(x_p,y_p)-f(x_{a1},y_{a1})|}{2}$\n邻近点C1(xc1，yc1)对应的灰度差\n$f_{c1}=\frac{|f(x_{d1},y_{d1})-f(x_{c1},y_{c1})|}{2}+\frac{|f(x_{c1},y_{c1})-f(x_p,y_p)|}{2}$\n邻近点D1(xd，yd)对应的灰度差\n$f_{d1}=\frac{|f(x_{e1},y_{e1})-f(x_{d1},y_{d1})|}{2}+\frac{|f(x_{d1},y_{d1})-f(x_{c1},y_{c1})|}{2}$\n(2)若边缘像素点P(xp，yp)在第2区域内，\n邻近点A2(xa2，ya2)对应的灰度差\n$f_{a2}=\frac{|f(x_p,y_p)-f(x_{a2},y_{a2})|}{2}+\frac{|f(x_{a2},y_{a2})-f(x_{b2},y_{b2})|}{2}$\n邻近点B2(xb2，yb2)对应的灰度差\n$f_{b2}=\frac{|f(x_{a2},y_{a2})-f(x_{b2},y_{b2})|}{2}+\frac{|f(x_{b2},y_{b2})-f(x_{f2},y_{f2})|}{2}$\n边缘点P(xp，yp)对应的灰度差\n$f_{p2}=\frac{|f(x_{c2},y_{c2})-f(x_p,y_p)|}{2}+\frac{|f(x_p,y_p)-f(x_{a2},y_{a2})|}{2}$\n邻近点C2(xc2，yc2)对应的灰度差\n$f_{c2}=\frac{|f(x_{d2},y_{d2})-f(x_{c2},y_{c2})|}{2}+\frac{|f(x_{c2},y_{c2})-f(x_p,y_p)|}{2}$\n邻近点D2(xd2，yd2)对应的灰度差\n$f_{d2}=\frac{|f(x_{e2},y_{e2})-f(x_{d2},y_{d2})|}{2}+\frac{|f(x_{d2},y_{d2})-f(x_{c2},y_{c2})|}{2}$\n第3.5步：求取亚像素点边缘点与像素边缘点的距离差，求取方法及说明如下：根据方形孔径原理，某一像素的灰度值输出可以表示为\n$f(i,j)={\int }_{j-0.5}^{j+0.5}{\int }_{i-0.5}^{i+0.5}g(x,y)\mathrm{dxdy}$\n式中，g(x，y)是连续图像的光强分布，f(i，j)是像素感光面上各部分光强综合作用的结果，采样结果是以灰度值为数值的离散矩阵。由于光学元器件的卷积作用以及光学衍射作用，在物空间剧变的灰度经光学成像成为渐变的形式，边缘在图像中表征为一种灰度分布，图像边缘灰度值变化应当是高斯分布，高斯曲线顶点对应的位置即为真实边缘点的位置，而对于圆形图像，像素级边缘点的梯度方向为圆心到该点的直线方向，仅需对像素级边缘点及其沿梯度方向的邻近点进行高斯曲线拟合即可求出亚像素边缘点的位置，这样可以将二维的高斯曲面拟合转化为一维的高斯曲线拟合。如图6所示，P(xp，yp)点为像素级边缘点，LG为边缘像素点P(xp，yp)的梯度方向直线，若P在第1区域内，B1、A1、C1、D1为P点沿梯度方向的邻近点，若P在第2区域内，B2、A2、C2、D2为P点沿梯度方向的邻近点，高斯曲线顶点M所对应的P′点应为其真实边缘点位置，P点与P′点的距离差为δ。一维高斯曲线的表达式为：\n$\tilde{y}=\frac{1}{\sqrt{2\pi }{\sigma }_G}\exp \left(\frac{-{(\delta -\mu )}^2}{2{{\sigma }_G}^2}\right)---\left(1\right)$\n式中，μ为高斯函数的均值，σG为高斯函数的标准差。\n为方便计算，对上式两边取对数，并令${\tilde{y}}^{*}=\ln \tilde{y},$上式可转化为：\n${\tilde{y}}^{*}=m_{11}{\delta }^2+m_{12}\delta +m_{13}---\left(2\right)$\n根据方形孔径采样定理，像素灰度差值为\n$\tilde{y}*\left(n\right)={\int }_{n-0.5}^{n+0.5}(m_{11}{\delta }^2+m_{12}\delta +m_{13})\mathrm{d\delta }---\left(3\right)$\n令边缘点P的序号为0，若P在第1区域内，P点对应的灰度差为fp1，其邻近点D1，C1，A1，B1序号分别表示为-2、-1、1和2，对应的灰度差为fd1、fc1、fa1和fb1，根据式(3)，有：\n$f_{d1}={\int }_{-2.5}^{-1.5}(m_{11}{\delta }^2+m_{12}\delta +m_{13})\mathrm{d\delta }=\frac{49}{12}{m}_{11}-2m_{12}+m_{13}---\left(4\right)$\n$f_{c1}={\int }_{-1.5}^{-0.5}(m_{11}{\delta }^2+m_{12}\delta +m_{13})\mathrm{d\delta }=\frac{13}{12}{m}_{11}-m_{12}+m_{13}---\left(5\right)$\n$f_{p1}={\int }_{-0.5}^{0.5}(m_{11}{\delta }^2+m_{12}\delta +m_{13})\mathrm{d\delta }=\frac{1}{12}{m}_{11}+m_{13}---\left(6\right)$\n$f_{a1}={\int }_{0.5}^{1.5}(m_{11}{\delta }^2+m_{12}\delta +m_{13})\mathrm{d\delta }=\frac{13}{12}{m}_{11}+m_{12}+m_{13}---\left(7\right)$\n$f_{b1}={\int }_{1.5}^{2.5}(m_{11}{\delta }^2+m_{12}\delta +m_{13})\mathrm{d\delta }=\frac{49}{12}{m}_{11}+2m_{12}+m_{13}---\left(8\right)$\n根据式(4)～(8)，联立方程组，用最小二乘法可求得m11，m12，m13关于fd1、fc1、fp1、fa1和fb1的表达式，将其代入抛物线顶点坐标值δ＝-m12/2m11，得：\n$\delta =-\frac{-0.2f_{d1}-0.1f_{c1}+0.1f_{a1}-0.2f_{b1}}{2(0.1429f_{d1}-0.0714f_{c1}-0.1429f_{p1}-0.0714f_{a1}+0.1429f_{b1}}$\n注意式(1)到式(2)的过程中，经过了取自然对数的操作，因此，像素灰度差值也应取对数，所以P点与P′点的距离差为：\n$\delta =\frac{0.1\ln f_{d1}+0.05\ln f_{c1}-0.05\ln f_{a1}+0.1\ln f_{b1}}{0.1429\ln f_{d1}-0.0714\ln f_{c1}-0.1429\ln f_{p1}-0.0714\ln f_{a1}+0.1429\ln f_{b1}}$\n若P在第2区域内，P点对应的灰度差为fp2，其邻近点D2，C2，A2，B2序号分别表示为-2、-1、1和2，对应的灰度差为fd2、fc2、fa2和fb2，根据式(3)，有：\n$f_{d2}={\int }_{-2.5}^{-1.5}(m_{11}{\delta }^2+m_{12}\delta +m_{13})\mathrm{d\delta }=\frac{49}{12}{m}_{11}-2m_{12}+m_{13}---\left(9\right)$\n$f_{c2}={\int }_{-1.5}^{-0.5}(m_{11}{\delta }^2+m_{12}\delta +m_{13})\mathrm{d\delta }=\frac{13}{12}{m}_{11}-m_{12}+m_{13}---\left(10\right)$\n$f_{p2}={\int }_{-0.5}^{0.5}(m_{11}{\delta }^2+m_{12}\delta +m_{13})\mathrm{d\delta }=\frac{1}{12}{m}_{11}+m_{13}---\left(11\right)$\n$f_{a2}={\int }_{0.5}^{1.5}(m_{11}{\delta }^2+m_{12}\delta +m_{13})\mathrm{d\delta }=\frac{13}{12}{m}_{11}+m_{12}+m_{13}---\left(12\right)$\n$f_{b2}={\int }_{1.5}^{2.5}(m_{11}{\delta }^2+m_{12}\delta +m_{13})\mathrm{d\delta }=\frac{49}{12}{m}_{11}+2m_{12}+m_{13}---\left(13\right)$\n根据式(9)～(13)，联立方程组，用最小二乘法可求得m11，m12，m13关于fd2、fc2、fp2、fa2和fb2的表达式，将其代入抛物线顶点坐标值δ＝-m12/2m11，得：\n$\delta =-\frac{-0.2f_{d2}-0.1f_{c2}+0.1f_{a2}-0.2f_{b2}}{2(0.1429f_{d2}-0.0714f_{c2}-0.1429f_{p2}-0.0714f_{a2}+0.1429f_{b2}}$\n注意式(1)到式(2)的过程中，经过了取自然对数的操作，因此，像素灰度差值也应取对数，所以P点与P′点的距离差为：\n$\delta =\frac{0.1\ln f_{d2}+0.05\ln f_{c2}-0.05\ln f_{a2}+0.1\ln f_{b2}}{0.1429\ln f_{d2}-0.0714\ln f_{c2}-0.1429\ln f_{p2}-0.0714\ln f_{a2}+0.1429\ln f_{b2}}$\n综上所述，根据第3.4步中所得边缘点及其沿梯度方向的邻近点的对应的灰度差，求取亚像素边缘点P′与像素边缘点P的距离差δ为\n\n若边缘点及其沿梯度方向的邻近点对应的灰度差中有一个为0，则令δ＝0，并存储δ＝0出现的次数Count_zero。\n第3.6步：求取亚像素边缘点P′的坐标，具体方法为：根据3.5步所求亚像素边缘点P′与像素级边缘点P的距离差δ以及像素级边缘点P的梯度方向直线的斜率kgradient，分别求出P点与P′点在x方向和y方向上的坐标差δx，δy，公式如下：\n${\delta }_x=\frac{\delta }{\sqrt{k_{\mathrm{grad}\mathrm{ient}}^2+1}}$\n${\delta }_y=\frac{k_{\mathrm{grad}\mathrm{ient}}*\delta }{\sqrt{k_{\mathrm{grad}\mathrm{ient}}^2+1}}$\n其中，像素级边缘点P的梯度方向直线的斜率kgradient可由该像素边缘点的坐标(xp，yp)和粗定位圆心坐标(xoc，yoc)获得，即：\n$k_{\mathrm{grad}\mathrm{ient}}=\frac{y_{\mathrm{oc}}-y_p}{x_{\mathrm{oc}}-x_p}$\n那么，对于像素点P(xp，yp)，其对应的亚像素边缘点为P′(xp+δx，yp+δy)，\n第3.7步：对所有边缘像素点P(xp，yp)，依据第3.2步～第3.6步依次求取对应的亚像素边缘点为P′(xp+δx，yp+δy)，其中，(xp+δx，yp+δy)为亚像素边缘点P′的坐标。\n第4步：对亚像素边缘点进行滤波，针对亚像素边缘点出现的“孤立点”和噪声，分别用曲率滤波的方法和均值滤波的方法进行处理。\n第4.1步：针对亚像素边缘出现的“孤立点”，用曲率滤波的方法进行处理，具体步骤如下：\n(1)求取所有亚像素边缘点的曲率，方法如下：在顺序存储的亚像素边缘点中，为亚像素边缘点P′的坐标，为P′点的紧邻前一点的坐标，为\nP′点的紧邻后一点的坐标，那么亚像素边缘点P′的曲率为\n$k=\frac{1}{r}=\frac{1}{\sqrt{({(x_0-x_{p_2^{\prime }})}^2+{(y_0-y_{p_2^{\prime }})}^2)}}$\n其中，(x0，y0)为P′、P′点的紧邻前一点、P′点的紧邻后一点所组成圆的圆心，\n$x_0=\frac{a-b+c}{d},$ $y_0=\frac{e-f+g}{-d},$\n$a=(x_{p_1^{\prime }}+x_{p_2^{\prime }})(x_{p_2^{\prime }}-x_{p_1^{\prime }})(y_{p_3^{\prime }}-y_{p_2^{\prime }})$\n$b=(x_{p_2^{\prime }}+x_{p_3^{\prime }})(x_{p_3^{\prime }}-x_{p_2^{\prime }})(y_{p_2^{\prime }}-y_{p_1^{\prime }})$\n$c=(y_{p_1^{\prime }}-y_{p_3^{\prime }})(y_{p_2^{\prime }}-y_{p_1^{\prime }})(y_{p_3^{\prime }}-y_{p_2^{\prime }})$\n$d=2[(x_{p_2^{\prime }}-x_{p_1^{\prime }})(y_{p_3^{\prime }}-y_{p_2^{\prime }})-(x_{p_3^{\prime }}-x_{p_2^{\prime }})(y_{p_2^{\prime }}-y_{p_1^{\prime }})]$\n$e=(y_{p_1^{\prime }}+y_{p_2^{\prime }})(y_{p_2^{\prime }}-y_{p_1^{\prime }})(x_{p_3^{\prime }}-x_{p_2^{\prime }})$\n$f=(y_{p_2^{\prime }}+y_{p_3^{\prime }})(y_{p_3^{\prime }}-y_{p_2^{\prime }})(x_{p_2^{\prime }}-x_{p_1^{\prime }})$\n$g=(x_{p_1^{\prime }}-x_{p_3^{\prime }})(x_{p_2^{\prime }}-x_{p_1^{\prime }})(x_{p_3^{\prime }}-x_{p_2^{\prime }})$\n(2)根据所求亚像素边缘点的曲率对所有亚像素边缘点进行滤波，滤波准则如下：\n首先将各亚像素边缘点的曲率进行降序排列，取曲率阈值为该序列曲率值由高到低的第3*n个元素，n为孤立点的个数，而造成“孤立点”的原因是因为在第3.5步中，边缘像素点及其沿梯度方向的邻近点所对应的灰度差有一个或多个为0，因此，孤立点的个数n数值为第3.5步中亚像素边缘定位时δ＝0的出现次数Count_zero；然后用曲率阈值对亚像素边缘点的曲率进行分割，如果亚像素边缘点P的曲率大于曲率阈值，且P的曲率大于其前后紧邻点的曲率，则认为点P是孤立点，给予滤除。\n第4.2步：将滤除“孤立点”后的所有亚像素边缘点用均值方法进行滤波，具体方法为：选取亚像素边缘点前后紧邻的两点的横、纵坐标与亚像素边缘点本身的横、纵坐标的均值作为该亚像素边缘点的新的横、纵坐标。\n第5步：对滤波后的亚像素边缘点用最小二乘法拟合圆，最终得到圆形目标的圆心和半径。最小二乘拟合圆的算法如下：\n设滤波后的亚像素边缘点个数为N_sub，二次曲线的一般方程为\nx2+2Bexy+Cey2+2Dex+2Eey+Fe＝0\n利用N_sub个边缘点进行曲线拟合时，其均方差和为\n$e^2=\underset{i=1}{\overset{N\_\mathrm{sub}}{\Sigma }}{(x_i^2+2B_e{x}_i{y}_i+C_e{y}_i^2+2D_e{x}_i+2E_e{y}_i+F_e)}^2$\n对上式关于Be，Ce，De，Ee，Fe分别取偏导数，并令每个式子为零，可以得到一个包含5个方程和5个未知数的静定方程组。用矩阵求逆或高斯列主元消去等方法可以求得最终的中心坐标(xcf，ycf)为：\n$x_{\mathrm{cf}}=\frac{B_e{E}_e-C_e{D}_e}{C_e-{B_e}^2}$\n$y_{\mathrm{cf}}=\frac{B_e{D}_e-E_e}{C_e-{B_e}^2}$