CN101334263A

CN101334263A - Method for locating the center of a circular target

Info

Publication number: CN101334263A
Application number: CNA2008100227931A
Authority: CN
Inventors: 达飞鹏; 张虎
Original assignee: Southeast University
Current assignee: Nantong Outpace Building Material Equipment Co ltd; Southeast University
Priority date: 2008-07-22
Filing date: 2008-07-22
Publication date: 2008-12-31
Anticipated expiration: 2028-07-22
Also published as: CN101334263B

Abstract

A center positioning method of a circular target mainly involves many target recognition and target positioning. Firstly, the simple center of circumference method is used to roughly locate the center of the image, and the key square area is extracted according to the information of the center of the circle and the radius of the rough positioning circle as In the area of interest, use the canny operator to perform pixel-level edge positioning on the circular target in the interested area, and then perform sub-pixel positioning on the circular target according to the geometric characteristics of the circle and its grayscale information to obtain accurate sub-pixel edge points The coordinates of the sub-pixel edge points, and then for the "isolated points" and noise appearing in the sub-pixel edge points, use the curvature filtering method and the mean value filtering method to filter respectively, and finally use the least square method to fit the circle to the filtered sub-pixel edge points, so as to obtain The final center and radius of the circle. This method not only effectively improves the accuracy of the positioning of the center of the circle, but also improves the robustness of the positioning of the center of the circle, thereby improving the measurement accuracy and stability of the measurement system.

Description

Method for locating the center of a circular target

技术领域 technical field

本发明涉及一种圆形目标的亚像素边缘定位以及亚像素边缘点滤波的方法，尤其涉及一种圆形目标的圆心定位方法。The invention relates to a sub-pixel edge positioning method of a circular target and a sub-pixel edge point filtering method, in particular to a circle center positioning method of a circular target.

背景技术 Background technique

在众多图像分析和处理系统中，快速而准确地获取圆心被广泛地应用在目标识别以及定位上，如在三维重构的系统中，高精度地获取标定板圆心对整个三维重构精度有至关重要的影响。目前求取圆心的常见算法有：形心法，Hough变换法，高斯曲面拟合法，最小二乘曲线拟合法。形心法要求图像灰度分布比较均匀，否则会产生较大误差；Hough变换法需要对每一个边界点逐点投票、记录，计算时间较长，而且占用计算机内存较大，在实际应用中受到限制；高斯曲面拟合是利用圆形光斑的灰度分布近似于高斯模型的原理进行曲面拟合，其定位精度高，但是如果处理圆面积较大，则运算量很大，计算时间长，不适用于实时光学检测系统中；最小二乘曲线拟合法是最为常见的方法，因为其实现简单，计算量小，在实际中得到广泛应用，但是目前最小二乘曲线拟合法定位圆中往往都是利用求取的像素级边缘点进行最小二乘法拟合圆心，使得其拟合精度受到了一定的限制，所以能否获取精确的图像中圆形的亚像素边缘，对于提高圆心的拟合精度有着重要的影响。In many image analysis and processing systems, quickly and accurately obtaining the center of the circle is widely used in target recognition and positioning. For example, in a 3D reconstruction system, obtaining the center of the calibration plate with high precision is of great importance to the accuracy of the entire 3D reconstruction. important influence. At present, the common algorithms for calculating the center of the circle are: centroid method, Hough transform method, Gaussian surface fitting method, and least squares curve fitting method. The centroid method requires that the gray distribution of the image is relatively uniform, otherwise there will be a large error; the Hough transform method needs to vote and record each boundary point point by point, which takes a long time to calculate and takes up a large amount of computer memory. Limitation; Gaussian surface fitting is to use the principle that the gray distribution of the circular spot is similar to the Gaussian model for surface fitting. Its positioning accuracy is high, but if the area of the processing circle is large, the amount of calculation is large, and the calculation time is long. It is suitable for real-time optical detection systems; the least squares curve fitting method is the most common method, because it is simple to implement and has a small amount of calculation, it is widely used in practice, but the current least squares curve fitting method is often used in the positioning circle Using the obtained pixel-level edge points to fit the center of the circle with the least squares method limits the fitting accuracy to a certain extent. Therefore, whether the accurate sub-pixel edge of the circle in the image can be obtained is of great importance for improving the fitting accuracy of the center of the circle. important influence.

目前常见的亚像素边缘定位算法有：参数拟合法，矩的方法，插值法以及利用图像其他几何特征的方法。文献“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-Zemike 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算子对绝大多数图像都不能实现整像素级边缘定位，导致亚像素定位不准；第二，该文献提供的方法不能确保二次曲线型的边缘定位精度。At present, the common sub-pixel edge location algorithms include: parameter fitting method, moment method, interpolation method and methods using other geometric features of the image. The document "Subpixel edge localization and the interpolation of still images" (Jensen K, Anastassiou, IEEE Trans on PAMI, 1995, 4(3): 285-295) is an interpolation method, which uses the pixel-level edge localization algorithm for interpolation and amplification. Image, under the premise that the actual image grayscale distribution law is consistent with the interpolation law, the edge positioning can reach sub-pixel accuracy; the literature "Using tangent direction information to detect sub-pixel edges" (Zhang Yujin, Fu Zhuo, Pattern Recognition and Artificial Intelligence, 1997 , 10(1): 69～74) mainly used the information of tangent direction to locate the sub-pixel edge of the target of known shape. Degree value information; literature "A fast subpixel edge detection method using Sobel-Zemike moments operator" (Qu Ying-Dong, Cui Cheng-Song, Chen San-Ben, LiJin-Quan, Image and Vision Computing, 2005, 23: 11～17 ) uses the moment method to realize the edge sub-pixel positioning. Moment-based methods are currently widely used, but they are not suitable for real-time processing of large bitmaps due to the use of multiple 7×7 templates; Laser, 2005, 16(8): 993～996) introduced the principle of sub-pixel edge parameter fitting, which uses Sobel operator and parameter fitting method to realize sub-pixel edge positioning, and there are two problems, the first First, the Sobel operator cannot achieve integer-pixel-level edge positioning for most images, resulting in inaccurate sub-pixel positioning; second, the method provided in this document cannot ensure the edge positioning accuracy of the quadratic curve.

发明内容 Contents of the invention

针对现有技术所存在的缺点和限制，本发明的目的在于提供一种能够提高测量系统精度的圆形目标的圆心定位方法。In view of the shortcomings and limitations of the prior art, the purpose of the present invention is to provide a method for locating the center of a circular target that can improve the accuracy of the measurement system.

本发明设计一种首先将图像中的圆形目标进行圆心粗定位，在用canny算子获取圆形目标的像素级边缘点后，综合考虑圆的几何特征及边缘灰度分布特征求取亚像素边缘点坐标，然后对精确的亚像素边缘点进行滤波，最后用最小二乘法拟合滤波后的亚像素边缘点求得圆形目标的圆心的方法。本发明采用如下技术方案：The present invention designs a method of roughly locating the center of the circular target in the image, and after obtaining the pixel-level edge points of the circular target with the canny operator, comprehensively consider the geometric characteristics of the circle and the edge gray distribution characteristics to obtain sub-pixels Edge point coordinates, then filter the precise sub-pixel edge points, and finally use the least squares method to fit the filtered sub-pixel edge points to obtain the center of the circular target. The present invention adopts following technical scheme:

一种圆形目标的圆心定位方法，其主要步骤为：A method for locating the center of a circular target, the main steps of which are:

第1步：对图像中圆形目标的圆心和半径进行粗定位，得到圆形目标的粗定位圆心(x_oc，y_oc)和粗定位圆半径R_c，其中，(x_oc，y_oc)为图像坐标系下圆心的坐标，具体方法如下：Step 1: Perform rough positioning on the center and radius of the circular target in the image, and obtain the rough positioning center (x _oc , y _oc ) and rough positioning circle radius R _c of the circular target, where (x _oc , y _oc ) is the coordinates of the center of the circle in the image coordinate system, the specific method is as follows:

第1.1步：对图像进行除噪、阈值分割，得到图像中各像素点灰度值为255或0的二值图像；Step 1.1: Carry out noise removal and threshold segmentation on the image to obtain a binary image with a gray value of 255 or 0 for each pixel in the image;

第1.2步：对二值图像中的圆形目标进行边界提取及边界跟踪；Step 1.2: Perform boundary extraction and boundary tracking on the circular target in the binary image;

第1.3步：对跟踪的边界点用周线形心法拟合圆，得到该圆形目标的粗定位圆心和粗定位圆半径，并存储粗定位圆心坐标(x_oc，y_oc)和粗定位圆半径R_c；Step 1.3: Fit the circle with the centroid method of the tracked boundary point, obtain the rough positioning circle center and the rough positioning circle radius of the circular target, and store the coarse positioning circle center coordinates (x _oc , y _oc ) and the rough positioning circle Radius R _c ;

第2步：对图像中的圆形目标进行像素级边缘定位，具体步骤如下：Step 2: Perform pixel-level edge positioning on the circular target in the image, the specific steps are as follows:

第2.1步：根据粗定位圆心和粗定位圆半径，从图像中提取一个正方形区域，称为“感兴趣区域”，该正方形区域的提取方法为：以图像中圆形目标的粗定位圆心作为正方形区域的中心点，以圆形目标的粗定位直径加上2～6个像素点作为正方形区域的边长；Step 2.1: According to the center of the rough positioning circle and the radius of the rough positioning circle, extract a square area from the image, called the "region of interest". The center point of the area, the rough positioning diameter of the circular target plus 2 to 6 pixels as the side length of the square area;

第2.2步：用canny算子对图像中所提取的感兴趣区域进行像素级边缘检测，得到圆形目标的像素级边缘，再对canny算子提取的圆形目标的像素级边缘点进行边界跟踪，得到每个像素级边缘点的坐标，且按顺时针顺序存储每个像素级边缘点的坐标；Step 2.2: Use the canny operator to perform pixel-level edge detection on the region of interest extracted in the image to obtain the pixel-level edge of the circular target, and then perform boundary tracking on the pixel-level edge points of the circular target extracted by the canny operator , get the coordinates of each pixel-level edge point, and store the coordinates of each pixel-level edge point in clockwise order;

第3步：对图像中的圆形目标进行亚像素边缘定位，具体步骤如下：Step 3: Perform sub-pixel edge positioning on the circular target in the image, the specific steps are as follows:

第3.1步：将圆形目标的边缘像素点进行区域划分，划分方法如下：称以圆形目标的粗定位圆心为原点、以水平向右为正方向、以像素为单位长度的直线为x轴；称以圆形目标的粗定位圆心为原点、以垂直向上为正方向、以像素为单位长度的直线为y轴；以粗定位圆心为旋转点，由x轴正方向逆时针旋转的角度为θ，根据θ的大小将圆形目标的边缘像素点划分为两个部分，一部分是θ∈[0°，45°]∪[135°，225°]∪[315°，360°]所对应的边缘像素点，称为第1区域的边缘像素点，另一部分是θ∈(45°，135°)∪(225°，315°)对应的边缘像素点，称为第2区域的边缘像素点，粗定位圆心与边缘像素点的连线为该边缘像素点的梯度方向；Step 3.1: Divide the edge pixels of the circular target into regions, and the division method is as follows: the x-axis is the straight line with the rough positioning center of the circular target as the origin, the positive direction horizontally to the right, and the unit length in pixels ; It is called the y-axis with the coarse positioning center of the circular target as the origin, the vertical upward as the positive direction, and the length in units of pixels as the y-axis; the coarse positioning center as the rotation point, and the angle of counterclockwise rotation from the positive direction of the x-axis is θ, according to the size of θ, the edge pixels of the circular target are divided into two parts, one part is θ∈[0°, 45°]∪[135°, 225°]∪[315°, 360°] corresponding to The edge pixels are called the edge pixels of the first area, and the other part is the edge pixels corresponding to θ∈(45°, 135°)∪(225°, 315°), which are called the edge pixels of the second area. The line connecting the center of the coarse positioning circle and the edge pixel is the gradient direction of the edge pixel;

第3.2步：求取边缘像素点沿梯度方向的邻近点，(x_p，y_p)为以像素为单位的边缘点P的坐标值，具体方法如下：若边缘像素点P(x_p，y_p)在第1区域内，取该像素点的梯度方向直线与直线x＝x_p+1的交点A₁(x_a1，y_a1)为该边缘像素点沿梯度方向的右边第一邻近点；取该像素点的梯度方向直线与直线x＝x_p+2的交点B₁(x_b1，y_b1)为该边缘像素点沿梯度方向的右边第二邻近点；取该像素点的梯度方向直线与直线x＝x_p+3的交点F(x_f，y_f)为该边缘像素点沿梯度方向的右边第三邻近点；取该像素点的梯度方向直线与直线x＝x_p-1的交点C₁(x_c1，y_c1)为该边缘像素点沿梯度方向的左边第一邻近点；取该像素点的梯度方向直线与直线x＝x_p-2的交点D₁(x_d1，y_d1)为该边缘像素点沿梯度方向的左边第二邻近点，取该像素点的梯度方向直线与直线x＝x_p-3的交点E₁(x_e1，y_e1)为该边缘像素点沿梯度方向的左边第三邻近点；若边缘像素点P(x_p，y_p)在第2区域内，取该像素点的梯度方向直线与直线y＝y_p+1的交点A₂(x_a2，y_a2)为该边缘像素点沿梯度方向的上边第一紧邻点；取该像素点的梯度方向直线与直线y＝y_p+2的交点B₂(x_b2，y_b2)为该边缘像素点沿梯度方向的上边第二紧邻点；取该像素点的梯度方向直线与直线y＝y_p+3的交点F₂(x_f2，y_f2)为该边缘像素点沿梯度方向的上边第三紧邻点；取该像素点的梯度方向直线与直线y＝y_p-1的交点C₂(x_c2，y_c2)为该边缘像素点沿梯度方向的下边第一紧邻点；取该像素点的梯度方向直线与直线y＝y_p-2的交点D₂(x_d2，y_d2)为该边缘像素点沿梯度方向的上边第二紧邻点；取该像素点的梯度方向直线与直线y＝y_p-3的交点E₂(x_e2，y_e2)为该边缘像素点沿梯度方向的上边第三紧邻点；Step 3.2: Calculate the adjacent points of the edge pixel point along the gradient direction, (x _p , y _p ) is the coordinate value of the edge point P in units of pixels, the specific method is as follows: if the edge pixel point P(x _p , y _p ) In the first area, take the intersection point A ₁ (x _a1 , y _a1 ) of the gradient direction line of the pixel point and the line x=x _p +1 as the first adjacent point on the right side of the edge pixel point along the gradient direction; Take the intersection point B ₁ (x _b1 , y _b1 ) of the gradient direction straight line of the pixel point and the straight line x=x _p +2 as the second adjacent point on the right side of the edge pixel point along the gradient direction; take the gradient direction straight line of the pixel point The point of intersection F(x _f , y _f ) with the straight line x=x _p +3 is the third adjacent point on the right side of the edge pixel along the gradient direction; take the gradient direction straight line of the pixel point and the straight line x=x _p -1 Intersection point C ₁ (x _c1 , y _c1 ) is the first adjacent point on the left side of the edge pixel along the gradient direction; take the intersection point D ₁ (x _d1 , _y _d1 ) is the second adjacent point on the left side of the edge pixel along the gradient direction, and the intersection point E ₁ (x _e1 , y _e1 ) of the gradient direction line of the pixel point and the line x=x _p -3 is the edge of the edge pixel point The third adjacent point on the left side of the gradient direction; if the edge pixel point P(x _p , y _p ) is in the second area, _take the intersection point A ₂ (x _a2 , y _a2 ) is the first adjacent point above the edge pixel along the gradient direction; take the intersection point B ₂ (x _b2 , y _b2 ) of the gradient direction line of the pixel point and the line y=y _p +2 as the edge pixel The second point next to the upper side of the point along the gradient direction; take the intersection point F ₂ (x _f2 , y _f2 ) of the gradient direction line of the pixel point and the straight line y=y _p +3 as the third point of the upper side of the edge pixel point along the gradient direction Close to the point; take the intersection point C ₂ (x _c2 , y _c2 ) of the gradient direction straight line of the pixel point and the straight line y=y _p -1 as the first adjacent point of the edge pixel point along the gradient direction; take the pixel point The intersection point D ₂ (x _d2 , y _d2 ) of the gradient direction straight line and the straight line y=y _p -2 is the second adjacent point on the upper edge of the edge pixel along the gradient direction; take the gradient direction straight line of the pixel point and the straight line y=y The intersection point E ₂ (x _e2 , y _e2 ) of _p -3 is the third immediate point above the edge pixel along the gradient direction;

第3.3步：采用线性灰度插值的方法求取以像素为单位的边缘像素点P(x_p，y_p)沿梯度方向的邻近点的灰度值，其中，(x_p，y_p)为以像素为单位的边缘像素点P的坐标值，使用f(x，y)表示坐标为(x，y)的像素点的灰度值，使用符号[]表示取整数部分，所述的用于获取邻近点灰度值的线性灰度插值方法如下：Step 3.3: Use the method of linear grayscale interpolation to obtain the grayscale value of the adjacent points of the edge pixel point P(x _p , y _p ) along the gradient direction in units of pixels, where (x _p , y _p ) is The coordinate value of the edge pixel point P in pixels, use f(x, y) to represent the gray value of the pixel point whose coordinates are (x, y), use the symbol [] to represent the integer part, the above is used for The linear gray interpolation method to obtain the gray value of adjacent points is as follows:

(1)若边缘像素点P(x_p，y_p)在第1区域内，(1) If the edge pixel point P(x _p , y _p ) is in the first area,

邻近点A₁(x_a1，y_a1)的灰度值Gray value of adjacent point A ₁ (x _a1 , y _a1 )

f(x_a1，y_a1)＝(1-λ)*f(x_a1，[y_a1])+λ*f(x_a1，[y_a1]+1)f(x _a1 , y _a1 )=(1-λ)*f(x _a1 , [y _a1 ])+λ*f(x _a1 , [y _a1 ]+1)

λ＝y_a1-[y_a1]，λ=y _a1 -[y _a1 ],

邻近点B₁(x_b1，y_b1)的灰度值Gray value of adjacent point B ₁ (x _b1 , y _b1 )

f(x_b1，y_b1)＝(1-λ)*f(x_b1，[y_b1])+λ*f(x_b1，[y_b1]+1)f(x _b1 , y _b1 )=(1-λ)*f(x _b1 ,[y _b1 ])+λ*f(x _b1 ,[y _b1 ]+1)

λ＝y_b1-[y_b1]λ=y _b1 -[y _b1 ]

邻近点F₁(x_f1，y_f1)的灰度值Gray value of adjacent point F ₁ (x _f1 , y _f1 )

f(x_f1，y_f1)＝(1-λ)*f(x_f1，[y_f1])+λ*f(x_f1，[y_f1]+1)f(x _f1 , y _f1 )=(1-λ)*f(x _f1 , [y _f1 ])+λ*f(x _f1 , [y _f1 ]+1)

λ＝y_f1-[y_f1]λ=y _f1 -[y _f1 ]

邻近点C₁(x_c1，y_c1)的灰度值Gray value of adjacent point C ₁ (x _c1 , y _c1 )

f(x_c1，y_c1)＝(1-λ)*f(x_c1，[y_c1])+λ*f(x_c1，[y_c1]+1)f(x _c1 , y _c1 )=(1-λ)*f(x _c1 ,[y _c1 ])+λ*f(x _c1 ,[y _c1 ]+1)

λ＝y_c1-[y_c1]λ=y _c1 -[y _c1 ]

邻近点D₁(x_d，y_d)的灰度值Gray value of adjacent point D ₁ (x _d , y _d )

f(x_d1，y_d1)＝(1-λ)*f(x_d1，[y_d1])+λ*f(x_d1，[y_d1]+1)f(x _d1 , y _d1 )=(1-λ)*f(x _d1 , [y _d1 ])+λ*f(x _d1 , [y _d1 ]+1)

λ＝y_d1-[y_d1]λ=y _d1 -[y _d1 ]

邻近点E₁(x_e1，y_e1)的灰度值Gray value of adjacent point E ₁ (x _e1 , y _e1 )

f(x_e1，y_e1)＝(1-λ)*f(x_e1，[y_e1])+λ*f(x_e1，[y_e1]+1)f(x _e1 , y _e1 )=(1-λ)*f(x _e1 ,[y _e1 ])+λ*f(x _e1 ,[y _e1 ]+1)

λ＝y_e1-[y_e1]λ=y _e1 -[y _e1 ]

(2)若边缘像素点P(x_p，y_p)在第2区域内，(2) If the edge pixel point P(x _p , y _p ) is in the second area,

邻近点A₂(x_a2，y_a2)的灰度值Gray value of adjacent point A ₂ (x _a2 , y _a2 )

f(x_a2，y_a2)＝(1-λ)*f([x_a2]，y_a2)+λ*f([x_a2]+1，y_a2)f(x _a2 , y _a2 )=(1-λ)*f([x _a2 ], y _a2 )+λ*f([x _a2 ]+1, y _a2 )

λ＝x_a2-[x_a2]λ=x _a2 -[x _a2 ]

邻近点B₂(x_b2，y_b2)的灰度值Gray value of adjacent point B ₂ (x _b2 , y _b2 )

f(x_b2，y_b2)＝(1-λ)*f([x_b2]，y_b2)+λ*f([x_b2]+1，y_b2)f(x _b2 , y _b2 )=(1-λ)*f([x _b2 ], y _b2 )+λ*f([x _b2 ]+1, y _b2 )

λ＝x_b2-[x_b2]λ=x _b2 -[x _b2 ]

邻近点F₂(x_f2，y_f2)的灰度值Gray value of adjacent point F ₂ (x _f2 , y _f2 )

f(x_f2，y_f2)＝(1-λ)*f([x_f2]，y_f2)+λ*f([x_f2]+1，y_f2)f(x _f2 , y _f2 )=(1-λ)*f([x _f2 ], y _f2 )+λ*f([x _f2 ]+1, y _f2 )

λ＝x_f2-[x_f2]λ=x _f2 -[x _f2 ]

邻近点C₂(x_c2，y_c2)的灰度值Gray value of adjacent point C ₂ (x _c2 , y _c2 )

f(x_c2，y_c2)＝(1-λ)*f([x_c2]，y_c2)+λ*f([x_c2]+1，y_c2)f(x _c2 , y _c2 )=(1-λ)*f([x _c2 ], y _c2 )+λ*f([x _c2 ]+1, y _c2 )

λ＝x_c2-[x_c2]λ=x _c2 -[x _c2 ]

邻近点D₂(x_d2，y_d2)的灰度值Gray value of adjacent point D ₂ (x _d2 , y _d2 )

f(x_d2，y_d2)＝(1-λ)*f([x_d2]，y_d2)+λ*f([x_d2]+1，y_d2)f(x _d2 , y _d2 )=(1-λ)*f([x _d2 ], y _d2 )+λ*f([x _d2 ]+1, y _d2 )

λ＝x_d2-[x_d2]λ=x _d2 -[x _d2 ]

邻近点E₂(x_e2，y_e2)的灰度值Gray value of adjacent point E ₂ (x _e2 , y _e2 )

f(x_e2，y_e2)＝(1-λ)*f([x_e2]，y_e2)+λ*f([x_e2]+1，y_e2)f(x _e2 , y _e2 )=(1-λ)*f([x _e2 ], y _e2 )+λ*f([x _e2 ]+1, y _e2 )

λ＝x_e2-[x_e2]λ=x _e2 -[x _e2 ]

第3.4步：根据第3.3步所求边缘点沿梯度方向的邻近点的灰度值以及边缘点P(x_p，y_p)本身的灰度值，选取灰度值的前向差分和后向差分的平均值作为边缘点及其沿梯度方向的邻近点的对应的灰度差，其中，(x_p，y_p)为以像素为单位的边缘像素点P的坐标值，使用f(x，y)表示坐标为(x，y)的像素点的灰度值，所述的求取灰度差的方法如下：Step 3.4: According to the gray value of the adjacent points of the edge point along the gradient direction obtained in step 3.3 and the gray value of the edge point P(x _p , y _p ) itself, select the forward difference and backward direction of the gray value The average value of the difference is taken as the corresponding gray level difference between the edge point and its adjacent points along the gradient direction, where (x _p , y _p ) is the coordinate value of the edge pixel point P in pixels, using f(x, y) represents the gray value of the pixel whose coordinates are (x, y), and the method for obtaining the gray difference is as follows:

邻近点A₁(x_a1，y_a1)对应的灰度差The gray level difference corresponding to the adjacent point A ₁ (x _a1 , y _a1 )

${f f}_{a a 11} = = \frac{| | f f (({x x}_{p p},, {y the y}_{p p})) - - f f (({x x}_{a a 11},, {y the y}_{a a 11})) | |}{22} + + \frac{| | f f (({x x}_{a a 11},, {y the y}_{a a 11})) - - f f (({x x}_{b b 11},, {y the y}_{b b 11})) | |}{22}$

邻近点B₁(x_b1，y_b1)对应的灰度差The gray level difference corresponding to the adjacent point B ₁ (x _b1 , y _b1 )

${f f}_{b b 11} = = \frac{| | f f (({x x}_{a a 11},, {y the y}_{a a 11})) - - f f (({x x}_{b b 11},, {y the y}_{b b 11})) | |}{22} + + \frac{| | f f (({x x}_{b b 11},, {y the y}_{b b 11})) - - f f (({x x}_{f f 11},, {y the y}_{f f 11})) | |}{22}$

边缘点P(x_p，y_p)对应的灰度差The gray level difference corresponding to the edge point P(x _p , y _p )

${f f}_{p p 11} = = \frac{| | f f (({x x}_{c c 11},, {y the y}_{c c 11})) - - f f (({x x}_{p p},, {y the y}_{p p})) | |}{22} + + \frac{| | f f (({x x}_{p p},, {y the y}_{p p})) - - f f (({x x}_{a a 11},, {y the y}_{a a 11})) | |}{22}$

邻近点C₁(x_c1，y_c1)对应的灰度差The gray level difference corresponding to the adjacent point C ₁ (x _c1 , y _c1 )

${f f}_{c c 11} = = \frac{| | f f (({x x}_{d d 11},, {y the y}_{d d 11})) - - f f (({x x}_{c c 11},, {y the y}_{c c 11})) | |}{22} + + \frac{| | f f (({x x}_{c c 11},, {y the y}_{c c 11})) - - f f (({x x}_{p p},, {y the y}_{p p})) | |}{22}$

邻近点D₁(x_d，y_d)对应的灰度差The gray level difference corresponding to the adjacent point D ₁ (x _d , y _d )

${f f}_{d d 11} = = \frac{| | f f (({x x}_{e e 11},, {y the y}_{e e 11})) - - f f (({x x}_{d d 11},, {y the y}_{d d 11})) | |}{22} + + \frac{| | f f (({x x}_{d d 11},, {y the y}_{d d 11})) - - f f (({x x}_{c c 11},, {y the y}_{c c 11})) | |}{22}$

邻近点A₂(x_a2，y_a2)对应的灰度差The gray level difference corresponding to the adjacent point A ₂ (x _a2 , y _a2 )

${f f}_{a a 22} = = \frac{| | f f (({x x}_{p p},, {y the y}_{p p})) - - f f (({x x}_{a a 22},, {y the y}_{a a 22})) | |}{22} + + \frac{| | f f (({x x}_{a a 22},, {y the y}_{a a 22})) - - f f (({x x}_{b b 22},, {y the y}_{b b 22})) | |}{22}$

邻近点B₂(x_b2，y_b2)对应的灰度差The gray level difference corresponding to the adjacent point B ₂ (x _b2 , y _b2 )

${f f}_{b b 22} = = \frac{| | f f (({x x}_{a a 22},, {y the y}_{a a 22})) - - f f (({x x}_{b b 22},, {y the y}_{b b 22})) | |}{22} + + \frac{| | f f (({x x}_{b b 22},, {y the y}_{b b 22})) - - f f (({x x}_{f f 22},, {y the y}_{f f 22})) | |}{22}$

${f f}_{p p 22} = = \frac{| | f f (({x x}_{c c 22},, {y the y}_{c c 22})) - - f f (({x x}_{p p},, {y the y}_{p p})) | |}{22} + + \frac{| | f f (({x x}_{p p},, {y the y}_{p p})) - - f f (({x x}_{a a 22},, {y the y}_{a a 22})) | |}{22}$

邻近点C₂(x_c2，y_c2)对应的灰度差The gray level difference corresponding to the adjacent point C ₂ (x _c2 , y _c2 )

${f f}_{c c 22} = = \frac{| | f f (({x x}_{d d 22},, {y the y}_{d d 22})) - - f f (({x x}_{c c 22},, {y the y}_{c c 22})) | |}{22} + + \frac{| | f f (({x x}_{c c 22},, {y the y}_{c c 22})) - - f f (({x x}_{p p},, {y the y}_{p p})) | |}{22}$

邻近点D₂(x_d2，y_d2)对应的灰度差The gray level difference corresponding to the adjacent point D ₂ (x _d2 , y _d2 )

${f f}_{d d 22} = = \frac{| | f f (({x x}_{e e 22},, {y the y}_{e e 22})) - - f f (({x x}_{d d 22},, {y the y}_{d d 22})) | |}{22} + + \frac{| | f f (({x x}_{d d 22},, {y the y}_{d d 22})) - - f f (({x x}_{c c 22},, {y the y}_{c c 22})) | |}{22}$

第3.5步：根据第3.4步中所得边缘点及其沿梯度方向的邻近点的对应的灰度差，求取亚像素边缘点P′与像素边缘点P的距离差δ为Step 3.5: Calculate the distance difference δ between the sub-pixel edge point P′ and the pixel edge point P according to the corresponding gray level difference between the edge point and its adjacent points along the gradient direction obtained in step 3.4.

若边缘点及其沿梯度方向的邻近点对应的灰度差中有一个为0，则令δ＝0，并存储δ＝0出现的次数Count_zero；If one of the grayscale differences corresponding to the edge point and its adjacent points along the gradient direction is 0, then let δ=0, and store the number of occurrences of δ=0 Count_zero;

第3.6步：求取亚像素边缘点P′的坐标，具体方法为：根据3.5步所求亚像素边缘点P′与像素级边缘点P的距离差δ以及像素级边缘点P的梯度方向直线的斜率k_gradient，分别求出P点与P′点在x方向和y方向上的坐标差δ_x，δ_y，公式如下：Step 3.6: Find the coordinates of the sub-pixel edge point P′, the specific method is: according to the distance difference δ between the sub-pixel edge point P′ and the pixel-level edge point P obtained in step 3.5 and the gradient direction straight line of the pixel-level edge point P The slope k _gradient , respectively calculate the coordinate difference δ _x , δ _y between point P and point P′ in the x direction and y direction, the formula is as follows:

${δ δ}_{x x} = = \frac{δ δ}{\sqrt{{k k}_{grad grad ient ient}^{22} + + 11}}$

${δ δ}_{y the y} = = \frac{{k k}_{grad grad ient ient} * * δ δ}{\sqrt{{k k}_{grad grad ient ient}^{22} + + 11}}$

其中，像素级边缘点P的梯度方向直线的斜率k_gradient可由该像素边缘点的坐标(x_p，y_p)和粗定位圆心坐标(x_oc，y_oc)获得，即：Among them, the slope k _gradient of the gradient direction line of the pixel-level edge point P can be obtained from the coordinates (x _p , y _p ) of the pixel edge point and the coordinates (x _oc , y _oc ) of the coarse positioning circle, namely:

${k k}_{grad grad ient ient} = = \frac{{y the y}_{oc oc} - - {y the y}_{p p}}{{x x}_{oc oc} - - {x x}_{p p}}$

那么，对于像素点P(x_p，y_p)，其对应的亚像素边缘点为P′(x_p+δ_x，y_p+δ_y)；Then, for a pixel point P(x _p , y _p ), its corresponding sub-pixel edge point is P′(x _p +δ _x , y _p +δ _y );

第3.7步：对所有边缘像素点P(x_p，y_p)，依据第3.2步～第3.6步依次求取对应的亚像素边缘点为P′(x_p+δ_x，y_p+δ_y)，其中，(x_p+δ_x，y_p+δ_y)为亚像素边缘点P′的坐标；Step 3.7: For all edge pixel points P(x _p , y _p ), according to steps 3.2 to 3.6, obtain the corresponding sub-pixel edge points as P′(x _p +δ _x , y _p +δ _y ), where (x _p +δ _x , y _p +δ _y ) is the coordinate of the sub-pixel edge point P′;

第4步：对亚像素边缘点进行滤波，针对亚像素边缘点出现的“孤立点”和噪声，分别用曲率滤波的方法和均值滤波的方法进行处理；Step 4: Filter the sub-pixel edge points, and use the curvature filtering method and the mean value filtering method to process the "isolated points" and noises that appear at the sub-pixel edge points;

第4.1步：针对亚像素边缘出现的“孤立点”，用曲率滤波的方法进行处理，具体步骤如下：Step 4.1: For the "isolated point" that appears on the edge of the sub-pixel, use the method of curvature filtering to process it. The specific steps are as follows:

(1)求取所有亚像素边缘点的曲率，方法如下：在顺序存储的亚像素边缘点中，(x′_p2，y′_p2)为亚像素边缘点P′的坐标，(x′_p1，y′_p1)为P′点的紧邻前一点的坐标，(x′_p3，y′_p3)为P′点的紧邻后一点的坐标，那么亚像素边缘点P′的曲率为(1) Calculate the curvature of all sub-pixel edge points, the method is as follows: in the sequentially stored sub-pixel edge points, (x′ _p2 , y′ _p2 ) is the coordinate of the sub-pixel edge point P′, (x′ _p1 , y′ _p1 ) is the coordinate of the point immediately before P′, (x′ _p3 , y′ _p3 ) is the coordinate of the point immediately after P′, then the curvature of the sub-pixel edge point P′ is

$k k = = \frac{11}{r r} = = \frac{11}{\sqrt{(({(({x x}_{00} - - {x x}_{{p p}_{22}^{' '}}))}^{22} + + {(({y the y}_{00} - - {y the y}_{{p p}_{22}^{' '}}))}^{22}))}}$

其中，(x₀，y₀)为P′、P′点的紧邻前一点、P′点的紧邻后一点所组成圆的圆心，Among them, (x ₀ , y ₀ ) is the center of the circle formed by P', the point immediately before P', and the point immediately after P',

${x x}_{00} = = \frac{a a - - b b + + c c}{d d},,$ ${y the y}_{00} = = \frac{e e - - f f + + g g}{- - d d},,$

a＝(x′_p1+x′_p2)(x′_p2-x′_p1)(y′_p3-y′_p2)a=(x′ _p1 +x′ _p2 )(x′ _p2 -x′ _p1 )(y′ _p3 -y′ _p2 )

b＝(x′_p2+x′_p3)(x′_p3-x′_p2)(y′_p2-y′_p1)b=(x′ _p2 +x′ _p3 )(x′ _p3 -x′ _p2 )(y′ _p2 -y′ _p1 )

c＝(y′_p1-y′_p3)(y′_p2-y′_p1)(y′_p3-y′_p2)c=(y′ _p1 -y′ _p3 )(y′ _p2 -y′ _p1 )(y′ _p3 -y′ _p2 )

d＝2[(x′_p2-x′_p1)(y′_p3-y′_p2)-(x′_p3-x′_p2)(y′_p2-y′_p1)]d=2[(x′ _p2 -x′ _p1 )(y′ _p3 -y′ _p2 )-(x′ _p3 -x′ _p2 )(y′ _p2 -y′ _p1 )]

e＝(y′_p1+y′_p2)(y′_p2-y′_p1)(x′_p3-x′_p2)e=(y′ _p1 +y′ _p2 )(y′ _p2 -y′ _p1 )(x′ _p3 -x′ _p2 )

f＝(y′_p2+y′_p3)(y′_p3-y′_p2)(x′_p2-x′_p1)f=(y′ _p2 +y′ _p3 )(y′ _p3 -y′ _p2 )(x′ _p2 -x′ _p1 )

g＝(x′_p1-x′_p3)(x′_p2-x′_p1)(x′_p3-x′_p2)g=(x′ _p1 -x′ _p3 )(x′ _p2 -x′ _p1 )(x′ _p3 -x′ _p2 )

(2)根据亚像素边缘点的曲率对所有亚像素边缘点进行滤波，滤波准则如下：(2) Filter all sub-pixel edge points according to the curvature of the sub-pixel edge points, and the filtering criteria are as follows:

首先将各亚像素边缘点的曲率进行降序排列，取曲率阈值为该序列曲率值由高到低的第3*n个元素，其中，n为孤立点的个数，其数值为第3.5步中亚像素边缘定位时δ＝0的出现次数Count_zero；然后用曲率阈值对亚像素边缘点的曲率进行分割，如果亚像素边缘点P的曲率大于曲率阈值，且P的曲率大于其前后紧邻点的曲率，则认为点P是孤立点，给予滤除；First, arrange the curvature of each sub-pixel edge point in descending order, and take the curvature threshold as the 3*nth element of the curvature value of the sequence from high to low, where n is the number of isolated points, and its value is in step 3.5 Count_zero the number of occurrences of δ=0 during sub-pixel edge positioning; then use the curvature threshold to segment the curvature of the sub-pixel edge point, if the curvature of the sub-pixel edge point P is greater than the curvature threshold, and the curvature of P is greater than the curvature of its immediate front and rear points , the point P is considered to be an isolated point and is filtered out;

第4.2步：将滤除“孤立点”后的所有亚像素边缘点用均值方法进行滤波，具体方法为：选取亚像素边缘点前后紧邻的两点的横、纵坐标与亚像素边缘点本身的横、纵坐标的均值作为该亚像素边缘点的新的横、纵坐标；Step 4.2: Filter all sub-pixel edge points after filtering out "isolated points" with the mean value method, the specific method is: select the horizontal and vertical coordinates of the two points immediately before and after the sub-pixel edge point and the sub-pixel edge point itself The mean value of the horizontal and vertical coordinates is used as the new horizontal and vertical coordinates of the sub-pixel edge point;

第5步：对滤波后的亚像素边缘点用最小二乘法拟合圆，最终得到圆形目标的圆心和半径。Step 5: Use the least square method to fit the circle to the filtered sub-pixel edge points, and finally get the center and radius of the circular target.

与现有技术相比，本发明具有如下优点：Compared with prior art, the present invention has following advantage:

(1)与传统的二维图像边缘定位方法相比，本发明提出的亚像素边缘定位的方法基于粗定位圆心所确定的感兴趣区域进行，而且在拟合边缘灰度差时将高斯曲面拟合简化为沿梯度方向的高斯曲线拟合，简化了计算，提高算法的快速性。(1) Compared with the traditional two-dimensional image edge positioning method, the sub-pixel edge positioning method proposed by the present invention is based on the region of interest determined by the rough positioning center, and the Gaussian surface is simulated when fitting the edge gray level difference The combination is simplified as a Gaussian curve fitting along the gradient direction, which simplifies the calculation and improves the speed of the algorithm.

(2)本发明提出的亚像素边缘定位的方法在考虑每个边缘像素点的梯度方向时，不拘泥于传统的边缘检测算子确定的几个特殊的方向，而是根据圆的特征考虑了每个边缘像素点的具体梯度方向，并且利用了线性灰度插值获取边缘像素点沿梯度方向的邻近点的灰度值，进而使边缘点的定位更加准确。(2) When considering the gradient direction of each edge pixel in the sub-pixel edge location method proposed by the present invention, it does not stick to several special directions determined by the traditional edge detection operator, but considers it according to the characteristics of the circle. The specific gradient direction of each edge pixel, and using linear grayscale interpolation to obtain the grayscale value of the adjacent points of the edge pixel along the gradient direction, so as to make the positioning of the edge point more accurate.

(3)本发明针对亚像素边缘点中出现的“孤立点”和噪声，提出了一种简单实用的滤波方法，该方法不仅原理和算法简单，而且滤波效果很好。(3) The present invention proposes a simple and practical filtering method aimed at "isolated points" and noise appearing in sub-pixel edge points. The method is not only simple in principle and algorithm, but also has a good filtering effect.

(4)本发明提出的圆形目标的圆心定位方法与已有的方法相比，有更高的精度和更好的稳定性，并且具有很强的实用性。(4) Compared with existing methods, the method for locating the center of a circular target proposed by the present invention has higher precision and better stability, and has strong practicability.

附图说明 Description of drawings

图1是圆形目标的圆心定位方法具体步骤的流程图。Fig. 1 is a flow chart of specific steps of a method for locating the center of a circular target.

图2是感兴趣区域提取图。Figure 2 is a map of region of interest extraction.

图3是非极大值抑制示意图，其中图3(a)为根据梯度方向划分的四个扇区图；3(b)为3*3邻域图。Figure 3 is a schematic diagram of non-maximum value suppression, in which Figure 3(a) is a four-sector map divided according to the gradient direction; 3(b) is a 3*3 neighborhood map.

图4(a)是边缘像素点P在第1区域时的定义像素点沿梯度方向邻近点的原理图；图4(b)是边缘像素点P在第2区域时的定义像素点沿梯度方向邻近点的原理图。Figure 4(a) is a schematic diagram of the adjacent points of the defined pixel point along the gradient direction when the edge pixel point P is in the first area; Figure 4(b) is the defined pixel point along the gradient direction when the edge pixel point P is in the second area Schematic of proximity points.

图5(a)是当边缘像素点P在第1区域时，边缘像素点P沿梯度方向的邻近点定义图；图5(b)是当边缘像素点P在第2区域时，边缘像素点P沿梯度方向的邻近点定义图。Figure 5(a) is a definition diagram of the adjacent points of the edge pixel point P along the gradient direction when the edge pixel point P is in the first area; Figure 5(b) is when the edge pixel point P is in the second area, the edge pixel point The graph is defined by the neighbors of P along the gradient direction.

图6是亚像素边缘点求取原理图。Fig. 6 is a schematic diagram of sub-pixel edge point calculation.

具体实施方式 Detailed ways

下面结合附图对本发明的具体实施方式作进一步描述。根据上述方法，在Windows操作系统中通过VC++6.0平台用C++编程实现了圆形目标的圆心定位操作。The specific embodiments of the present invention will be further described below in conjunction with the accompanying drawings. According to the above method, the center positioning operation of the circular target is realized by programming in C++ through the VC++6.0 platform in the Windows operating system.

应用本方法进行圆形目标的圆心定位主要包括圆心粗定位、像素级边缘定位、亚像素边缘定位、亚像素边缘点滤波、最小二乘拟合求得圆心五个操作步骤，具体步骤的流程图如图1所示，应用本方法进行圆形目标的圆心定位的具体步骤如下(步骤中未标明的单位都以像素为单位)：Applying this method to locate the center of a circular target mainly includes five operational steps: coarse center positioning, pixel-level edge positioning, sub-pixel edge positioning, sub-pixel edge point filtering, and least squares fitting to obtain the center of the circle. The flow chart of the specific steps As shown in Figure 1, the specific steps of applying this method to locate the center of a circular target are as follows (unmarked units in the steps are all in pixels):

第1步：对图像中的圆形目标的圆心和半径进行粗定位，得到圆形目标的粗定位圆心(x_oc，y_oc)和粗定位圆半径R_c，其中，(x_oc，y_oc)为图像坐标系下圆心的坐标，具体步骤如下：Step 1: Perform rough positioning on the center and radius of the circular target in the image, and obtain the rough positioning circle center (x _oc , y _oc ) and rough positioning circle radius R _c of the circular target, where (x _oc , y _oc ) is the coordinates of the center of the circle in the image coordinate system, and the specific steps are as follows:

第1.1步：对图像进行除噪、阈值分割，得到图像中各像素点灰度值为255或0的二值图像。其中，阈值分割采用迭代阈值法。迭代阈值法能自动计算出比较合适的分割阈值，其计算步骤如下：Step 1.1: Perform denoising and threshold segmentation on the image to obtain a binary image with a gray value of 255 or 0 for each pixel in the image. Among them, the threshold segmentation adopts the iterative threshold method. The iterative threshold method can automatically calculate a more appropriate segmentation threshold, and the calculation steps are as follows:

(1)选择阈值T，通常选择图像的平均灰度值作为初始阈值；(1) Select the threshold T, usually choose the average gray value of the image as the initial threshold;

(2)通过初始阈值T，将图像的平均灰度值分成为两组R₁和R₂；(2) by the initial threshold T, the average gray value of the image is divided into two groups R ₁ and R ₂ ;

(3)计算两组平均灰度值μ₁和μ₂；(3) Calculate the average gray value μ ₁ and μ ₂ of the two groups;

(4)重新选择阈值T，新的T定义为：T＝(μ₁+μ₂)/2；(4) Reselect the threshold T, and the new T is defined as: T=(μ ₁ +μ ₂ )/2;

循环(2)到(4)，直到两组的平均灰度值μ₁和μ₂不再发生改变，因此得到所需要的阈值，再根据所得到的阈值对图像进行二值化。Cycle (2) to (4) until the average gray value μ ₁ and μ ₂ of the two groups no longer change, so the required threshold is obtained, and then the image is binarized according to the obtained threshold.

第1.2步：对二值图像中的圆形目标进行边界提取及边界跟踪。在二值图像中，灰度值为255的像素即为圆形目标，对圆形目标进行边界提取的方法如下：如果原图中有一点为白，而相邻的8个像素全部为黑，则将该像素点的灰度值设为0。经过边界提取后，对该图中的目标轮廓线进行跟踪并存储，可得到圆形目标的边界点及其坐标。Step 1.2: Perform boundary extraction and boundary tracking on the circular target in the binary image. In a binary image, a pixel with a grayscale value of 255 is a circular target. The method of boundary extraction for a circular target is as follows: If one point in the original image is white, and the adjacent 8 pixels are all black, Then the gray value of the pixel is set to 0. After the boundary extraction, the target contour line in the figure is tracked and stored, and the boundary points and coordinates of the circular target can be obtained.

第1.3步：对跟踪的边界点用周线形心法拟合圆，得到该圆形目标的粗定位圆心和粗定位圆半径，并存储粗定位圆心坐标(x_oc，y_oc)和粗定位圆半径R_c。周线形心法是一种最为简单的亚像素定位算法，算法如下：设跟踪的边界点坐标为(i，j)，跟踪边界点总数为N，根据周线形心法，圆形目标的圆心坐标为：Step 1.3: Fit the circle with the centroid method of the tracked boundary point, obtain the rough positioning circle center and the rough positioning circle radius of the circular target, and store the coarse positioning circle center coordinates (x _oc , y _oc ) and the rough positioning circle Radius R _c . The contour centroid method is the simplest sub-pixel positioning algorithm. The algorithm is as follows: Set the coordinates of the tracking boundary point as (i, j), and the total number of tracking boundary points is N. According to the contour centroid method, the circle center coordinates of the circular target for:

$\{\begin{matrix} {x x}_{oc oc} = = \frac{Σi Σ i}{N N} \\ {y the y}_{oc oc} = = \frac{Σj Σj}{N N} \end{matrix}$

圆形目标的半径为各边界点到圆心距离的平均值，即：The radius of the circular target is the average distance from each boundary point to the center of the circle, namely:

${R R}_{c c} = = \frac{Σ Σ \sqrt{{((i i - - {x x}_{oc oc}))}^{22} + + {((j j - - {y the y}_{oc oc}))}^{22}}}{N N}$

第2.1步：根据粗定位圆心和粗定位圆半径，从图像中提取一个正方形区域，称为“感兴趣区域”，该正方形区域的提取方法为：以图像中圆形目标的粗定位圆心作为正方形区域的中心点，以圆形目标的粗定位直径加上2～6个像素点作为正方形区域的边长。如图2所示，图中圆代表图像中的圆形目标，正方形区域即为感兴趣区域，其中，O为粗定位圆心，R_c为粗定位半径，d一般取2～6个像素。后续步骤只针对感兴趣区域进行处理。Step 2.1: According to the center of the rough positioning circle and the radius of the rough positioning circle, extract a square area from the image, called the "region of interest". For the center point of the area, the rough positioning diameter of the circular target plus 2 to 6 pixels is used as the side length of the square area. As shown in Figure 2, the circle in the figure represents the circular target in the image, and the square area is the region of interest, where O is the center of the rough positioning circle, R _c is the radius of the rough positioning, and d generally takes 2 to 6 pixels. Subsequent steps only target regions of interest.

第2.2步：用canny算子对图像中所提取的感兴趣区域进行像素级边缘检测，得到圆形目标的像素级边缘，再对canny算子提取的圆形目标的像素级边缘点进行边界跟踪，得到每个像素级边缘点的坐标，且按顺时针顺序存储每个像素级边缘点的坐标。其中，Canny边缘检测对信噪比与定位乘积进行测度，属于先平滑后再求导的方法，其主要步骤和基本原理如下：Step 2.2: Use the canny operator to perform pixel-level edge detection on the region of interest extracted in the image to obtain the pixel-level edge of the circular target, and then perform boundary tracking on the pixel-level edge points of the circular target extracted by the canny operator , get the coordinates of each pixel-level edge point, and store the coordinates of each pixel-level edge point in clockwise order. Among them, Canny edge detection measures the signal-to-noise ratio and the positioning product, which belongs to the method of first smoothing and then derivation. The main steps and basic principles are as follows:

(1)用高斯滤波器平滑图像，高斯平滑函数如下：(1) Smooth the image with a Gaussian filter, and the Gaussian smoothing function is as follows:

$H h ((x x,, y the y)) = = {e e}^{- - \frac{{x x}^{22} + + {y the y}^{22}}{{22 σ σ}^{22}}}$

其中，σ为高斯函数的标准偏差，x，y分别是图像中像素点的横、纵坐标。对图像进行高斯滤波，即对图像f_I(x，y)用高斯函数H(x，y)进行卷积，得到滤波后的平滑图像G(x，y)：Among them, σ is the standard deviation of the Gaussian function, and x and y are the horizontal and vertical coordinates of the pixels in the image, respectively. Gaussian filtering is performed on the image, that is, the image f _I (x, y) is convolved with the Gaussian function H(x, y) to obtain a filtered smooth image G(x, y):

G(x，y)＝f_I(x，y)*H(x，y)G(x, y) = f _I (x, y)*H(x, y)

(2)用一阶偏导数的有限差分来计算梯度的幅值和方向。一阶差分卷积模板如下：(2) Use the finite difference of the first-order partial derivatives to calculate the magnitude and direction of the gradient. The first-order differential convolution template is as follows:

${H h}_{11} = = |\begin{matrix} - - 11 & - - 11 \\ 11 & 11 \end{matrix}|$ ${H h}_{22} = = |\begin{matrix} 11 & - - 11 \\ 11 & - - 11 \end{matrix}|$

将平滑后的图像G(x，y)用以上两个模板进行卷积，得到水平和垂直方向上的梯度图像分别为进而由

可得到梯度的幅值

和方向

推算式如下：The smoothed image G(x, y) is convolved with the above two templates, and the gradient images in the horizontal and vertical directions are obtained as And then by

The magnitude of the gradient can be obtained

and direction

The calculation formula is as follows:

(3)对梯度幅值进行非极大值抑制(non-maxima suppression，NMS)：图3是非极大值抑制示意图。根据梯度的方向，四个扇区的标号为0到3，如图3(a)所示，对应3*3邻域的四种可能组合，3*3邻域如图3(b)所示。在每一个像素点上，将邻域的中心象素P_cp与沿着梯度线的两个像素相比，如果P_cp的梯度值不比沿梯度线的两个相邻像素梯度值大，则该点不为边缘点。(3) Perform non-maxima suppression (NMS) on the gradient amplitude: FIG. 3 is a schematic diagram of non-maxima suppression. According to the direction of the gradient, the four sectors are labeled from 0 to 3, as shown in Figure 3(a), corresponding to four possible combinations of 3*3 neighborhoods, as shown in Figure 3(b) . At each pixel point, compare the central pixel P _cp of the neighborhood with the two pixels along the gradient line, if the gradient value of P _cp is not greater than the gradient value of the two adjacent pixels along the gradient line, then the Point is not an edge point.

(4)用双阈值算法检测和连接边缘。双阈值算法对非极大值抑制图像作用两个阈值τ1和τ2，且2τ1≈τ2，从而可以得到两个阈值边缘图像N₁[i，j]和N₂[i，j]。由于N₂[i，j]使用高阈值得到，因而含有很少的假边缘，但有间断。双阈值法要在N₂[i，j]中把边缘连接成轮廓，当到达轮廓的端点时，该算法就在N₁[i，j]的8邻域点中寻找可以连接到轮廓上的边缘，这样，算法不断地在N₁[i，j]中收集边缘，直到将N₂[i，j]连接起来为止。(4) Detect and connect edges with a double-threshold algorithm. The double-threshold algorithm applies two thresholds τ1 and τ2 to the non-maximum suppressed image, and 2τ1≈τ2, so that two threshold edge images N ₁ [i, j] and N ₂ [i, j] can be obtained. Since N ₂ [i, j] is obtained using a high threshold, it contains few false edges, but there are discontinuities. The double-threshold method needs to connect the edges into a contour in N ₂ [i, j]. When the end point of the contour is reached, the algorithm searches for edges that can be connected to the contour in the 8 neighborhood points of N ₁ [i, j]. , so that the algorithm keeps collecting edges in N ₁ [i, j] until N ₂ [i, j] is connected.

第3.1步：将圆形目标的边缘像素点进行区域划分，划分方法如下：称以圆形目标的粗定位圆心为原点、以水平向右为正方向、以像素为单位长度的直线为x轴；称以圆形目标的粗定位圆心为原点、以垂直向上为正方向、以像素为单位长度的直线为y轴；以粗定位圆心为旋转点，由x轴正方向逆时针旋转的角度为θ，根据θ的大小将圆形目标的边缘像素点划分为两个部分，一部分是θ∈[0°，45°]∪[135°，225°]∪[315°，360°]所对应的边缘像素点，称为第1区域的边缘像素点，另一部分是θ∈(45°，135°)∪(225°，315°)对应的边缘像素点，称为第2区域的边缘像素点，粗定位圆心与边缘像素点的连线为该边缘像素点的梯度方向。Step 3.1: Divide the edge pixels of the circular target into regions, and the division method is as follows: the x-axis is the straight line with the rough positioning center of the circular target as the origin, the positive direction horizontally to the right, and the unit length in pixels ; It is called the y-axis with the coarse positioning center of the circular target as the origin, the vertical upward as the positive direction, and the length in units of pixels as the y-axis; the coarse positioning center as the rotation point, and the angle of counterclockwise rotation from the positive direction of the x-axis is θ, according to the size of θ, the edge pixels of the circular target are divided into two parts, one part is θ∈[0°, 45°]∪[135°, 225°]∪[315°, 360°] corresponding to The edge pixels are called the edge pixels of the first area, and the other part is the edge pixels corresponding to θ∈(45°, 135°)∪(225°, 315°), which are called the edge pixels of the second area. The line connecting the center of the coarse positioning circle and the edge pixel point is the gradient direction of the edge pixel point.

第3.2步：求取边缘像素点沿梯度方向的邻近点，具体原理如下：若边缘像素点P(x_p，y_p)在第1区域内，图4(a)是定义像素点沿梯度方向邻近点的原理图，(x_p，y_p)为以像素为单位的边缘点P的坐标值，L_G为P点的梯度方向直线，直线y₁、y₂、y₃、y₄为P点附近以像素点为单位的四条横线，分别与直线L_G相交于点B′₁、A′₁、C′₁、D′₁，直线x₁，x₂，x₃，x₄为P点附近以像素点为单位的四条纵线，分别与直线L_G相交于点D₁、C₁、A₁、B₁，选取{B′₁、A′₁、C′₁、D′₁}和{B₁、A₁、C₁、D₁}两组点作为边缘点沿梯度方向的候选邻近点，由于P点在第1区域，梯度方向以x方向为主，即L_G的斜率的绝对值小于1，因此，候选邻近点{B′₁、A′₁、C′₁、D′₁}到P点的距离明显要大于候选邻近点{B₁、A₁、C₁、D₁}到P点的距离，若选用{B′₁、A′₁、C′₁、D′₁}作为P点沿梯度方向的邻近点，会导致边缘定位不够精确，因此，选用{B₁、A₁、C₁、D₁}作为其邻近点；若边缘像素点P(x_p，y_p)在第2区域内，图4(b)是定义像素点沿梯度方向邻近点的原理图，(x_p，y_p)为以像素为单位的边缘点P的坐标值，L_G为P点的梯度方向直线，直线y₁、y₂、y₃、y₄为P点附近以像素点为单位的四条横线，分别与直线I_G相交于点B₂、A₂、C₂、D₂，直线x₁，x₂，x₃，x₄为P点附近以像素点为单位的四条纵线，分别与直线L_G相交于点D′₂、C′₂、A′₂、B′₂，选取{B′₂、A′₂、C′₂、D′₂}和{B₂、A₂、C₂、D₂}两组点作为边缘点沿梯度方向的候选邻近点。由于P点在第2区域，梯度方向以y方向为主，即L_G的斜率的绝对值大于1，因此，候选邻近点{B′₂、A′₂、C′₂、D′₂}到P点的距离明显要大于候选邻近点{B₂、A₂、C₂、D₂}到P点的距离，若选用{B′₂、A′₂、C′₂、D′₂}作为P点沿梯度方向的邻近点，会导致边缘定位不够精确，因此，选用{B₂、A₂、C₂、D₂}作为其邻近点。因此，针对第3.1步划分的两个不同区域的像素点，其沿梯度方向的邻近点的定义如下：(x_p，y_p)为以像素为单位的边缘点P的坐标值，若边缘像素点P(x_p，y_p)在第1区域内，图5(a)是像素点沿梯度方向的邻近点定义图，L_G为边缘像素点P(x_p，y_p)的梯度方向直线，L₁代表直线x＝x_p-3，L₂代表直线x＝x_p-2，L₃代表直线x＝x_p-1，L₄代表直线x＝x_p+1，L₅代表直线x＝x_p+2，L₆代表直线x＝x_p+3 。取该像素点的梯度方向直线与直线x＝x_p+1的交点A₁(x_a1，y_a1)为该边缘像素点沿梯度方向的右边第一邻近点；取该像素点的梯度方向直线与直线x＝x_p+2的交点B₁(x_b1，y_b1)为该边缘像素点沿梯度方向的右边第二邻近点；取该像素点的梯度方向直线与直线x＝x_p+3的交点F(x_f，y_f)为该边缘像素点沿梯度方向的右边第三邻近点；取该像素点的梯度方向直线与直线x＝x_p-1的交点C₁(x_c1，y_c1)为该边缘像素点沿梯度方向的左边第一邻近点；取该像素点的梯度方向直线与直线x＝x_p-2的交点D₁(x_d1，y_d1)为该边缘像素点沿梯度方向的左边第二邻近点，取该像素点的梯度方向直线与直线x＝x_p-3的交点E₁(x_e1，y_e1)为该边缘像素点沿梯度方向的左边第三邻近点；若边缘像素点P(x_p，y_p)在第2区域内，图5(b)是像素点沿梯度方向的邻近点定义图，L_G为边缘像素点P(x_p，y_p)的梯度方向直线，取该像素点的梯度方向直线与直线y＝y_p+1的交点A₂(x_a2，y_a2)为该边缘像素点沿梯度方向的上边第一紧邻点；取该像素点的梯度方向直线与直线y＝y_p+2的交点B₂(x_b2，y_b2)为该边缘像素点沿梯度方向的上边第二紧邻点；取该像素点的梯度方向直线与直线y＝y_p+3的交点F₂(x_f2，y_f2)为该边缘像素点沿梯度方向的上边第三紧邻点；取该像素点的梯度方向直线与直线y＝y_p-1的交点C₂(x_c2，y_c2)为该边缘像素点沿梯度方向的下边第一紧邻点；取该像素点的梯度方向直线与直线y＝y_p-2的交点D₂(x_d2，y_d2)为该边缘像素点沿梯度方向的上边第二紧邻点；取该像素点的梯度方向直线与直线y＝y_p-3的交点E₂(x_e2，y_e2)为该边缘像素点沿梯度方向的上边第三紧邻点。Step 3.2: Calculate the adjacent points of the edge pixel along the gradient direction, the specific principle is as follows: if the edge pixel P(x _p , y _p ) is in the first area, Figure 4(a) is to define the pixel along the gradient direction The schematic diagram of adjacent points, (x _p , y _p ) is the coordinate value of the edge point P in pixels, L _G is the gradient direction straight line of point P, and the straight lines y ₁ , y ₂ , y ₃ , and y ₄ are P The four horizontal lines in pixel units near the point respectively intersect with the straight line L _G at points B′ ₁ , A′ ₁ , C′ ₁ , D′ ₁ , and the straight lines x ₁ , x ₂ , x ₃ , and x ₄ are P The four vertical lines in pixel units near the point respectively intersect with the straight line L _G at points D ₁ , C ₁ , A ₁ , and B ₁ , and select {B′ ₁ , A′ ₁ , C′ ₁ , D′ ₁ } and {B ₁ , A ₁ , C ₁ , D ₁ } two groups of points are used as candidate adjacent points of the edge point along the gradient direction. Since point P is in the first area, the gradient direction is dominated by the x direction, that is, the slope of _LG The absolute value is less than 1, therefore, the distance between the candidate neighboring points {B′ ₁ , A′ ₁ , C′ ₁ , D′ ₁ } and point P is obviously greater than that of the candidate neighboring points {B ₁ , A ₁ , C ₁ , D ₁ } to point P. If {B′ ₁ , A′ ₁ , C′ ₁ , D′ ₁ } are selected as the adjacent points of point P along the gradient direction, the edge positioning will not be accurate enough. Therefore, {B ₁ , A ₁ , C ₁ , D ₁ } are their neighboring points; if the edge pixel point P(x _p , y _p ) is in the second area, Figure 4(b) is a schematic diagram of defining the neighboring points of the pixel point along the gradient direction, (x _p , y _p ) is the coordinate value of the edge point P in units of pixels, L _G is the gradient direction straight line of point P, and the straight lines y ₁ , y ₂ , y ₃ , and y ₄ are the pixels near point P. The four horizontal lines of the unit intersect with the straight line I _G at points B ₂ , A ₂ , C ₂ , and D ₂ , and the straight lines x ₁ , x ₂ , x ₃ , and x ₄ are the four vertical lines in units of pixels near point P. Lines intersect with straight line L _G at points D′ ₂ , C′ ₂ , A′ ₂ , B′ ₂ , and select {B′ ₂ , A′ ₂ , C′ ₂ , D′ ₂ } and {B ₂ , A ₂ , C ₂ , D ₂ } two groups of points are used as candidate adjacent points of the edge point along the gradient direction. Since point P is in the second area, the gradient direction is dominated by the y direction, that is, the absolute value of the slope of L _G is greater than 1. Therefore, the candidate adjacent points {B′ ₂ , A′ ₂ , C′ ₂ , D′ ₂ } to The distance of point P is obviously greater than the distance from candidate neighboring points {B ₂ , A ₂ , C ₂ , D ₂ } to point P. If {B′ ₂ , A′ ₂ , C′ ₂ , D′ ₂ } is selected as P The adjacent points along the gradient direction will lead to inaccurate edge positioning. Therefore, {B ₂ , A ₂ , C ₂ , D ₂ } are selected as their adjacent points. Therefore, for the pixel points in the two different areas divided in step 3.1, the definition of the adjacent points along the gradient direction is as follows: (x _p , y _p ) is the coordinate value of the edge point P in units of pixels, if the edge pixel The point P(x _p , y _p ) is in the first area, Figure 5(a) is the definition map of the adjacent points of the pixel point along the gradient direction, L _G is the gradient direction straight line of the edge pixel point P(x _p , y _p ) , L ₁ represents the straight line x=x _p -3, L ₂ represents the straight line x=x _p -2, L ₃ represents the straight line x=x _p -1, L ₄ represents the straight line x=x _p +1, L ₅ represents the straight line x =x _p +2, L ₆ represents the straight line x=x _p +3 . Take the intersection point A ₁ (x _a1 , y _a1 ) of the gradient direction straight line of the pixel point and the straight line x=x _p +1 as the first adjacent point on the right side of the edge pixel point along the gradient direction; take the gradient direction straight line of the pixel point The intersection point B ₁ (x _b1 , y _b1 ) with the straight line x=x _p +2 is the second adjacent point on the right side of the edge pixel along the gradient direction; take the gradient direction straight line of the pixel point and the straight line x=x _p +3 _The intersection point F(x _f , y _f ) of the edge pixel point along the gradient direction is the third adjacent point on the right; take the intersection point C ₁ (x _c1 , y _c1 ) is the first adjacent point on the left side of the edge pixel point along the gradient direction; take the intersection point D ₁ (x _d1 , y _d1 ) of the gradient direction line of the pixel point and the line x=x _p -2 as the edge pixel point The second adjacent point on the left side of the gradient direction, take the intersection point E ₁ (x _e1 , y _e1 ) of the gradient direction line of the pixel point and the line x=x _p -3 as the third adjacent point on the left side of the edge pixel point along the gradient direction ; If the edge pixel point P(x _p , y _p ) is in the second area, Figure 5(b) is the definition map of the adjacent points of the pixel point along the gradient direction, L _G is the edge pixel point P(x _p , y _p ) Gradient direction straight line of the pixel point, take the intersection point A ₂ (x _a2 , y _a2 ) of the gradient direction straight line of the pixel point and the straight line y=y _p +1 as the first adjacent point on the upper edge of the edge pixel point along the gradient direction; take the pixel The intersection point B ₂ (x _b2 , y _b2 ) of the gradient direction straight line of the point and the straight line y=y _p +2 is the second adjacent point on the upper side of the edge pixel along the gradient direction; take the gradient direction straight line and the straight line y of the pixel point The intersection point F ₂ (x _f2 , y _f2 ) of =y _p +3 is the third adjacent point above the edge pixel along the gradient direction; take the intersection point C of the gradient direction line of the pixel point and the line y=y _p -1 ₂ (x _c2 , y _c2 ) is the first adjacent point on the lower edge of the edge pixel along the gradient direction; take the intersection point D ₂ (x _d2 , y _d2 ) of the gradient direction line of the pixel point and the line y=y _p -2 is the second adjacent point on the upper side of the edge pixel point along the gradient direction; take the intersection point E ₂ (x _e2 , y _e2 ) of the gradient direction line of the pixel point and the line y=y _p -3 as the edge pixel point along the gradient direction The third immediate point above .

第3.3步：采用线性灰度插值的方法求取以像素为单位的边缘像素点P(x_p，y_p)沿梯度方向的邻近点的灰度值，具体方法为：首先找出与该像素点的沿梯度方向的邻近点距离最近的两个整像素点，然后对这两个整像素点进行线性灰度插值，插值结果为该边缘像素点的沿梯度方向的邻近点的灰度值，其中，(x_p，y_p)为以像素为单位的边缘像素点P的坐标值，使用f(x，y)表示坐标为(x，y)的像素点的灰度值，使用符号[]表示取整数部分，所述的用于获取邻近点灰度值的线性灰度插值方法如下：Step 3.3: Use the method of linear grayscale interpolation to obtain the grayscale value of the adjacent point of the edge pixel point P(x _p , y _p ) along the gradient direction in units of pixels. The specific method is: first find out the The adjacent points along the gradient direction of the point are the nearest two integer pixels, and then linear grayscale interpolation is performed on these two integer pixels, and the interpolation result is the grayscale value of the adjacent points along the gradient direction of the edge pixel point, Among them, (x _p , y _p ) is the coordinate value of the edge pixel point P in units of pixels, f(x, y) is used to represent the gray value of the pixel point whose coordinates are (x, y), and the symbol [] Indicates the integer part, the linear gray interpolation method used to obtain the gray value of adjacent points is as follows:

λ＝y_a1-[y_a1]，λ=y _a1 -[y _a1 ],

λ＝y_b1-[y_b1]λ=y _b1 -[y _b1 ]

λ＝y_f1-[y_f1]λ=y _f1 -[y _f1 ]

λ＝y_c1-[y_c1]λ=y _c1 -[y _c1 ]

λ＝y_d1-[y_d1]λ=y _d1 -[y _d1 ]

λ＝y_e1-[y_e1]λ=y _e1 -[y _e1 ]

λ＝x_a2-[x_a2]λ=x _a2 -[x _a2 ]

λ＝x_b2-[x_b2]λ=x _b2 -[x _b2 ]

λ＝x_f2-[x_f2]λ=x _f2 -[x _f2 ]

λ＝x_c2-[x_c2]λ=x _c2 -[x _c2 ]

λ＝x_d2-[x_d2]λ=x _d2 -[x _d2 ]

f(x_e2，y_e2)＝(1-λ)*f([_e2]，y_e2)+λ*f([_e2]+1，y_e2)f(x _e2 , y _e2 )=(1-λ)*f([ _e2 ], y _e2 )+λ*f([ _e2 ]+1, y _e2 )

λ＝x_e2-[x_e2]λ=x _e2 -[x _e2 ]

第3.5步：求取亚像素点边缘点与像素边缘点的距离差，求取方法及说明如下：根据方形孔径原理，某一像素的灰度值输出可以表示为Step 3.5: Calculate the distance difference between the sub-pixel edge point and the pixel edge point. The calculation method and description are as follows: According to the square aperture principle, the gray value output of a certain pixel can be expressed as

$f f ((i i,, j j)) = = {&Integral; &Integral;}_{j j - - 0.5 0.5}^{j j + + 0.5 0.5} {&Integral; &Integral;}_{i i - - 0.5 0.5}^{i i + + 0.5 0.5} g g ((x x,, y the y)) dxdy dxdy$

式中，g(x，y)是连续图像的光强分布，f(i，j)是像素感光面上各部分光强综合作用的结果，采样结果是以灰度值为数值的离散矩阵。由于光学元器件的卷积作用以及光学衍射作用，在物空间剧变的灰度经光学成像成为渐变的形式，边缘在图像中表征为一种灰度分布，图像边缘灰度值变化应当是高斯分布，高斯曲线顶点对应的位置即为真实边缘点的位置，而对于圆形图像，像素级边缘点的梯度方向为圆心到该点的直线方向，仅需对像素级边缘点及其沿梯度方向的邻近点进行高斯曲线拟合即可求出亚像素边缘点的位置，这样可以将二维的高斯曲面拟合转化为一维的高斯曲线拟合。如图6所示，P(x_p，y_p)点为像素级边缘点，L_G为边缘像素点P(x_p，y_p)的梯度方向直线，若P在第1区域内，B₁、A₁、C₁、D₁为P点沿梯度方向的邻近点，若P在第2区域内，B₂、A₂、C₂、D₂为P点沿梯度方向的邻近点，高斯曲线顶点M所对应的P′点应为其真实边缘点位置，P点与P′点的距离差为δ。一维高斯曲线的表达式为：In the formula, g(x, y) is the light intensity distribution of the continuous image, f(i, j) is the result of the comprehensive effect of the light intensity of each part on the photosensitive surface of the pixel, and the sampling result is a discrete matrix with gray values. Due to the convolution effect of optical components and the effect of optical diffraction, the gray level that changes drastically in the object space becomes a gradual form through optical imaging, and the edge is represented as a gray level distribution in the image, and the change of the gray level value at the edge of the image should be a Gaussian distribution , the position corresponding to the vertex of the Gaussian curve is the position of the real edge point, and for a circular image, the gradient direction of the pixel-level edge point is the straight line direction from the center of the circle to the point, and only the pixel-level edge point and its gradient along the gradient direction The position of the sub-pixel edge point can be obtained by Gaussian curve fitting of adjacent points, so that the two-dimensional Gaussian surface fitting can be transformed into one-dimensional Gaussian curve fitting. As shown in Figure 6, P(x _p , y _p ) is a pixel-level edge point, L _G is the gradient direction straight line of edge pixel point P(x _p , y _p ), if P is in the first area, B ₁ , A ₁ , C ₁ , D ₁ are adjacent points of point P along the gradient direction. If P is in the second area, B ₂ , A ₂ , C ₂ , D ₂ are adjacent points of point P along the gradient direction. The point P' corresponding to the vertex M should be the position of the real edge point, and the distance difference between point P and point P' is δ. The expression of a one-dimensional Gaussian curve is:

$\overset{~ ~}{y the y} = = \frac{11}{\sqrt{22 π π} {σ σ}_{G G}} exp exp ((\frac{- - {((δ δ - - μ μ))}^{22}}{22 {σ σ}_{G G}^{22}})) - - - - - - ((11))$

式中，μ为高斯函数的均值，σ_G为高斯函数的标准差。In the formula, μ is the mean value of the Gaussian function, and σ _G is the standard deviation of the Gaussian function.

为方便计算，对上式两边取对数，并令 ${\tilde{y}}^{*} = \ln \tilde{y},$ 上式可转化为：For the convenience of calculation, logarithms are taken on both sides of the above formula, and let ${\tilde{the y}}^{*} = \ln \tilde{the y},$ The above formula can be transformed into:

$\overset{~ ~ * *}{y the y} = = {m m}_{1111} {δ δ}^{22} + + {m m}_{1212} δ δ + + {m m}_{1313} - - - - - - ((22))$

根据方形孔径采样定理，像素灰度差值为According to the square aperture sampling theorem, the pixel gray level difference is

$\overset{~ ~}{y the y} * * ((n no)) = = {&Integral; &Integral;}_{n no - - 0.5 0.5}^{n no + + 0.5 0.5} (({m m}_{1111} {δ δ}^{22} + + {m m}_{1212} δ δ + + {m m}_{1313})) dδ dδ - - - - - - ((33))$

令边缘点P的序号为0，若P在第1区域内，P点对应的灰度差为f_p1，其邻近点D₁，C₁，A₁，B₁序号分别表示为-2、-1、1和2，对应的灰度差为f_d1、f_c1、f_a1和f_b1，根据式(3)，有：Let the serial number of the edge point P be 0, if P is in the first area, the gray level difference corresponding to P point is f _p1 , and the serial numbers of its adjacent points D ₁ , C ₁ , A ₁ , B ₁ are expressed as -2, - 1, 1 and 2, the corresponding gray level difference is f _d1 , f _c1 , f _a1 and f _b1 , according to formula (3), we have:

${f f}_{d d 11} = = {&Integral; &Integral;}_{- - 2.5 2.5}^{- - 1.5 1.5} (({m m}_{1111} {δ δ}^{22} + + {m m}_{1212} δ δ + + {m m}_{1313})) dδ dδ = = \frac{4949}{1212} {m m}_{1111} - - {22 m m}_{1212} + + {m m}_{1313} - - - - - - ((44))$

${f f}_{c c 11} = = {&Integral; &Integral;}_{- - 1.5 1.5}^{- - 0.5 0.5} (({m m}_{1111} {δ δ}^{22} + + {m m}_{1212} δ δ + + {m m}_{1313})) dδ dδ = = \frac{1313}{1212} {m m}_{1111} - - {m m}_{1212} + + {m m}_{1313} - - - - - - ((55))$

${f f}_{p p 11} = = {&Integral; &Integral;}_{- - 0.5 0.5}^{0.5 0.5} (({m m}_{1111} {δ δ}^{22} + + {m m}_{1212} δ δ + + {m m}_{1313})) dδ dδ = = \frac{11}{1212} {m m}_{1111} + + {m m}_{1313} - - - - - - ((66))$

${f f}_{a a 11} = = {&Integral; &Integral;}_{0.5 0.5}^{1.5 1.5} (({m m}_{1111} {δ δ}^{22} + + {m m}_{1212} δ δ + + {m m}_{1313})) dδ dδ = = \frac{1313}{1212} {m m}_{1111} + + {m m}_{1212} + + {m m}_{1313} - - - - - - ((77))$

${f f}_{b b 11} = = {&Integral; &Integral;}_{1.5 1.5}^{2.5 2.5} (({m m}_{1111} {δ δ}^{22} + + {m m}_{1212} δ δ + + {m m}_{1313})) dδ dδ = = \frac{4949}{1212} {m m}_{1111} + + {22 m m}_{1212} + + {m m}_{1313} - - - - - - ((88))$

根据式(4)～(8)，联立方程组，用最小二乘法可求得m₁₁，m₁₂，m₁₃关于f_d1、f_c1、f_p1、f_a1和f_b1的表达式，将其代入抛物线顶点坐标值δ＝-m₁₂/2m₁₁，得：According to the equations (4)~(8) and the simultaneous equations, the expression of m ₁₁ , m ₁₂ , m ₁₃ about f _d1 , f _c1 , f _p1 , f _a1 and f _b1 can be obtained by using the least square method. Substituting it into the parabola vertex coordinate value δ=-m ₁₂ /2m ₁₁ , we get:

$δ δ = = - - \frac{- - 0.2 0.2 {f f}_{d d 11} - - 0.1 0.1 {f f}_{c c 11} + + 0.1 0.1 {f f}_{a a 11} - - 0.2 0.2 {f f}_{b b 11}}{22 (({0.1429 0.1429 f f}_{d d 11} - - 0.0714 0.0714 {f f}_{c c 11} - - 0.1429 0.1429 {f f}_{p p 11} - - 0.0714 0.0714 {f f}_{a a 11} + + 0.1429 0.1429 {f f}_{b b 11}}$

注意式(1)到式(2)的过程中，

经过了取自然对数的操作，因此，像素灰度差值也应取对数，所以P点与P′点的距离差为：Note that in the process from formula (1) to formula (2),

After the operation of taking natural logarithm, the pixel gray level difference should also take logarithm, so the distance difference between point P and point P′ is:

$δ δ = = \frac{0.1 0.1 ln ln {f f}_{} d d 11 + + 0.05 0.05 ln ln {f f}_{c c 11} - - 0.05 0.05 ln ln {f f}_{a a 11} + + 00 . . 11 ln ln {f f}_{b b 11}}{0.1429 0.1429 ln ln {f f}_{d d 11} - - 0.0714 0.0714 ln ln {f f}_{c c 11} - - 0.1429 0.1429 ln ln {f f}_{p p 11} - - 0.714 0.714 ln ln {f f}_{a a 11} + + 0.1429 0.1429 ln ln {f f}_{} b b 11}$

若P在第2区域内，P点对应的灰度差为f_p2，其邻近点D₂，C₂，A₂，B₂序号分别表示为-2、-1、1和2，对应的灰度差为f_d2、f_c2、f_a2和f_b2，根据式(3)，有：If P is in the second area, the gray level difference corresponding to point P is f _p2 , and the serial numbers of its adjacent points D ₂ , C ₂ , A ₂ , and B ₂ are respectively expressed as -2, -1, 1 and 2, and the corresponding gray Degree difference is f _d2 , f _c2 , f _a2 and f _b2 , according to formula (3), we have:

${f f}_{d d 22} = = {&Integral; &Integral;}_{- - 2.5 2.5}^{- - 1.5 1.5} (({m m}_{1111} {δ δ}^{22} + + {m m}_{1212} δ δ + + {m m}_{1313})) dδ dδ = = \frac{4949}{1212} {m m}_{1111} - - {22 m m}_{1212} + + {m m}_{1313} - - - - - - ((99))$

${f f}_{c c 22} = = {&Integral; &Integral;}_{- - 1.5 1.5}^{- - 0.5 0.5} (({m m}_{1111} {δ δ}^{22} + + {m m}_{1212} δ δ + + {m m}_{1313})) dδ dδ = = \frac{1313}{1212} {m m}_{1111} - - {m m}_{1212} + + {m m}_{1313} - - - - - - ((1010))$

${f f}_{p p 22} = = {&Integral; &Integral;}_{- - 0.5 0.5}^{00 . . 55} (({m m}_{1111} {δ δ}^{22} + + {m m}_{1212} δ δ + + {m m}_{1313})) dδ dδ = = \frac{11}{1212} {m m}_{1111} + + {m m}_{1313} - - - - - - ((1111))$

${f f}_{a a 22} = = {&Integral; &Integral;}_{00 . . 55}^{1.5 1.5} (({m m}_{1111} {δ δ}^{22} + + {m m}_{1212} δ δ + + {m m}_{1313})) dδ dδ = = \frac{1313}{1212} {m m}_{1111} + + {m m}_{1212} + + {m m}_{1313} - - - - - - ((1212))$

${f f}_{b b 22} = = {&Integral; &Integral;}_{1.5 1.5}^{2.5 2.5} (({m m}_{1111} {δ δ}^{22} + + {m m}_{1212} δ δ + + {m m}_{1313})) dδ dδ = = \frac{4949}{1212} {m m}_{1111} + + {22 m m}_{1212} + + {m m}_{1313} - - - - - - ((1313))$

根据式(9)～(13)，联立方程组，用最小二乘法可求得m₁₁，m₁₂，m₁₃关于f_d2、f_c2、f_p2、f_a2和f_b2的表达式，将其代入抛物线顶点坐标值δ＝-m₁₂/2m₁₁，得：According to equations (9)~(13), the simultaneous equations, the expression of m ₁₁ , m ₁₂ , m ₁₃ with respect to f _d2 , f _c2 , f _p2 , f _a2 and f _b2 can be obtained by using the least square method. Substituting it into the parabola vertex coordinate value δ=-m ₁₂ /2m ₁₁ , we get:

$δ δ = = - - \frac{- - 0.2 0.2 {f f}_{d d 22} - - 0.1 0.1 {f f}_{c c 22} + + 0.1 0.1 {f f}_{a a 22} - - 0.2 0.2 {f f}_{b b 22}}{22 ((0.1429 0.1429 {f f}_{d d 22} - - 0.0714 0.0714 {f f}_{c c 22} - - 0.1429 0.1429 {f f}_{p p 22} - - 0.0714 0.0714 {f f}_{a a 22} + + 0.1429 0.1429 {f f}_{b b 22}}$

注意式(1)到式(2)的过程中，

$δ δ = = \frac{0.1 0.1 ln ln {f f}_{} d d 22 + + 0.05 0.05 ln ln {f f}_{c c 22} - - 0.05 0.05 ln ln {f f}_{a a 22} + + 0.1 0.1 ln ln {f f}_{} b b 22}{0.1429 0.1429 ln ln {f f}_{d d 22} - - 0.0714 0.0714 ln ln {f f}_{c c 22} - - 0.1429 0.1429 ln ln {f f}_{p p 22} - - 0.0174 0.0174 ln ln {f f}_{a a 22} + + 0.142 0.142 {99 ln ln f f}_{b b 22}}$

综上所述，根据第3.4步中所得边缘点及其沿梯度方向的邻近点的对应的灰度差，求取亚像素边缘点P′与像素边缘点P的距离差δ为To sum up, according to the corresponding gray level difference between the edge point and its adjacent points along the gradient direction obtained in step 3.4, the distance difference δ between the sub-pixel edge point P′ and the pixel edge point P is calculated as

若边缘点及其沿梯度方向的邻近点对应的灰度差中有一个为0，则令δ＝0，并存储δ＝0出现的次数Count_zero。If one of the grayscale differences corresponding to the edge point and its adjacent points along the gradient direction is 0, then set δ=0, and store the number of occurrences of δ=0 as Count_zero.

第3.6步：求取亚像素边缘点P′的坐标，具体方法为：根据3.5步所求亚像素边缘点P与像素级边缘点P的距离差δ以及像素级边缘点P的梯度方向直线的斜率k_gradient，分别求出P点与P′点在x方向和y方向上的坐标差δ_x，δ_y，公式如下：Step 3.6: Calculate the coordinates of the sub-pixel edge point P′, the specific method is: according to the distance difference δ between the sub-pixel edge point P and the pixel-level edge point P obtained in step 3.5 and the gradient direction straight line of the pixel-level edge point P Slope k _gradient , respectively calculate the coordinate difference δ _x , δ _y between point P and point P′ in the x direction and y direction, the formula is as follows:

那么，对于像素点P(x_p，y_p)，其对应的亚像素边缘点为P′(x_p+δ_x，y_p+δ_y)，Then, for a pixel point P(x _p , y _p ), its corresponding sub-pixel edge point is P′(x _p +δ _x , y _p +δ _y ),

第3.7步：对所有边缘像素点P(x_p，y_p)，依据第3.2步～第3.6步依次求取对应的亚像素边缘点为P′(x_p+δ_x，y_p+δ_y)，其中，(x_p+δ_x，y_p+δ_y)为亚像素边缘点P′的坐标。Step 3.7: For all edge pixel points P(x _p , y _p ), according to steps 3.2 to 3.6, obtain the corresponding sub-pixel edge points as P′(x _p +δ _x , y _p +δ _y ), where (x _p +δ _x , y _p +δ _y ) is the coordinate of the sub-pixel edge point P′.

第4步：对亚像素边缘点进行滤波，针对亚像素边缘点出现的“孤立点”和噪声，分别用曲率滤波的方法和均值滤波的方法进行处理。Step 4: Filter the sub-pixel edge points, and use the curvature filtering method and the mean value filtering method to process the "isolated points" and noises that appear at the sub-pixel edge points, respectively.

(1)求取所有亚像素边缘点的曲率，方法如下：在顺序存储的亚像素边缘点中，为亚像素边缘点P′的坐标，

为P′点的紧邻前一点的坐标，

为P′点的紧邻后一点的坐标，那么亚像素边缘点P′的曲率为(1) Calculate the curvature of all sub-pixel edge points, the method is as follows: in the sequentially stored sub-pixel edge points, is the coordinate of the sub-pixel edge point P′,

is the coordinates of the immediately preceding point of point P′,

is the coordinate of the next point of point P′, then the curvature of the sub-pixel edge point P′ is

$a a = = (({x x}_{{p p}_{11}^{' '}} + + {x x}_{{p p}_{22}^{' '}})) (({x x}_{{p p}_{22}^{' '}} - - {x x}_{{p p}_{11}^{' '}})) (({y the y}_{{p p}_{33}^{' '}} - - {y the y}_{{p p}_{22}^{' '}}))$

$b b = = (({x x}_{{p p}_{22}^{' '}} + + {x x}_{{p p}_{33}^{' '}})) (({x x}_{{p p}_{33}^{' '}} - - {x x}_{{p p}_{22}^{' '}})) (({y the y}_{{p p}_{22}^{' '}} - - {y the y}_{{p p}_{11}^{' '}}))$

$c c = = (({y the y}_{{p p}_{11}^{' '}} - - {y the y}_{{p p}_{33}^{' '}})) (({y the y}_{{p p}_{22}^{' '}} - - {y the y}_{{p p}_{11}^{' '}})) (({y the y}_{{p p}_{33}^{' '}} - - {y the y}_{{p p}_{22}^{' '}}))$

$d d = = 22 [[(({x x}_{{p p}_{22}^{' '}} - - {x x}_{{p p}_{11}^{' '}})) (({y the y}_{{p p}_{33}^{' '}} - - {y the y}_{{p p}_{22}^{' '}})) - - (({x x}_{{p p}_{33}^{' '}} - - {x x}_{{p p}_{22}^{' '}})) (({y the y}_{{p p}_{22}^{' '}} - - {y the y}_{{p p}_{11}^{' '}}))]]$

$e e = = (({y the y}_{{p p}_{11}^{' '}} + + {y the y}_{{p p}_{22}^{' '}})) (({y the y}_{{p p}_{22}^{' '}} - - {y the y}_{{p p}_{11}^{' '}})) (({x x}_{{p p}_{33}^{' '}} - - {x x}_{{p p}_{22}^{' '}}))$

$f f = = (({y the y}_{{p p}_{22}^{' '}} + + {y the y}_{{p p}_{33}^{' '}})) (({y the y}_{{p p}_{33}^{' '}} - - {y the y}_{{p p}_{22}^{' '}})) (({x x}_{{p p}_{22}^{' '}} - - {x x}_{{p p}_{11}^{' '}}))$

$g g = = (({x x}_{{p p}_{11}^{' '}} - - {x x}_{{p p}_{33}^{' '}})) (({x x}_{{p p}_{22}^{' '}} - - {x x}_{{p p}_{11}^{' '}})) (({x x}_{{p p}_{33}^{' '}} - - {x x}_{{p p}_{22}^{' '}}))$

(2)根据所求亚像素边缘点的曲率对所有亚像素边缘点进行滤波，滤波准则如下：(2) Filter all sub-pixel edge points according to the curvature of the sub-pixel edge points. The filtering criteria are as follows:

首先将各亚像素边缘点的曲率进行降序排列，取曲率阈值为该序列曲率值由高到低的第3*n个元素，n为孤立点的个数，而造成“孤立点”的原因是因为在第3.5步中，边缘像素点及其沿梯度方向的邻近点所对应的灰度差有一个或多个为0，因此，孤立点的个数n数值为第3.5步中亚像素边缘定位时δ＝0的出现次数Count_zero；然后用曲率阈值对亚像素边缘点的曲率进行分割，如果亚像素边缘点P的曲率大于曲率阈值，且P的曲率大于其前后紧邻点的曲率，则认为点P是孤立点，给予滤除。First, the curvature of each sub-pixel edge point is arranged in descending order, and the curvature threshold is taken as the 3*nth element of the curvature value of the sequence from high to low, n is the number of isolated points, and the reason for the "isolated point" is Because in step 3.5, one or more grayscale differences corresponding to edge pixels and their adjacent points along the gradient direction are 0, therefore, the number n of isolated points is the sub-pixel edge positioning in step 3.5 When δ=0, the number of occurrences Count_zero; then use the curvature threshold to segment the curvature of the sub-pixel edge point, if the curvature of the sub-pixel edge point P is greater than the curvature threshold, and the curvature of P is greater than the curvature of its immediate front and rear points, then the point is considered P is an isolated point, which is filtered out.

第4.2步：将滤除“孤立点”后的所有亚像素边缘点用均值方法进行滤波，具体方法为：选取亚像素边缘点前后紧邻的两点的横、纵坐标与亚像素边缘点本身的横、纵坐标的均值作为该亚像素边缘点的新的横、纵坐标。Step 4.2: Filter all sub-pixel edge points after filtering out "isolated points" with the mean value method, the specific method is: select the horizontal and vertical coordinates of the two points immediately before and after the sub-pixel edge point and the sub-pixel edge point itself The average value of the horizontal and vertical coordinates is used as the new horizontal and vertical coordinates of the sub-pixel edge point.

第5步：对滤波后的亚像素边缘点用最小二乘法拟合圆，最终得到圆形目标的圆心和半径。最小二乘拟合圆的算法如下：Step 5: Use the least square method to fit the circle to the filtered sub-pixel edge points, and finally get the center and radius of the circular target. The least squares algorithm for fitting circles is as follows:

设滤波后的亚像素边缘点个数为N_sub，二次曲线的一般方程为Let the number of filtered sub-pixel edge points be N_sub, and the general equation of the quadratic curve is

x²+2B_exy+C_ey²+2D_ex+2E_ey+F_e＝0x ² +2B _e xy+C _e y ² +2D _e x+2E _e y+F _e ＝0

利用N_sub个边缘点进行曲线拟合时，其均方差和为When using N_sub edge points for curve fitting, the mean square error sum is

${e e}^{22} = = {Σ Σ}_{i i = = 11}^{N N__sub sub} {(({x x}_{}^{i i} + + 22 {B B}_{e e} {x x}_{i i} {y the y}_{i i} + + {C C}_{e e} {y the y}_{}^{i i} + + 22 {D D.}_{e e} {x x}_{i i} + + 22 {E E.}_{e e} {y the y}_{i i} + + {F f}_{e e}))}^{22}$

对上式关于B_e，C_e，D_e，E_e，F_e分别取偏导数，并令每个式子为零，可以得到一个包含5个方程和5个未知数的静定方程组。用矩阵求逆或高斯列主元消去等方法可以求得最终的中心坐标(x_cf，y_cf)为：Taking partial derivatives of the above formula with respect to Be _e , C _e , De _e , E _e , and F _e respectively, and setting each formula to be zero, a statically determinate equation system containing 5 equations and 5 unknowns can be obtained. The final center coordinates (x _cf , y _cf ) can be obtained by matrix inversion or Gaussian column pivot elimination:

${x x}_{cf cf} = = \frac{{B B}_{e e} {E E.}_{e e} - - {C C}_{e e} {D D.}_{e e}}{{C C}_{e e} - - {B B}_{e e}^{22}}$

${y the y}_{cf cf} = = \frac{{B B}_{e e} {D D.}_{e e} - - {E E.}_{e e}}{{C C}_{e e} - - {B B}_{e e}^{22}} . .$

Claims

1, a kind of circle center positioning method of circular target, it is characterized in that:

Step 1: Perform rough positioning on the center and radius of the circular target in the image, and obtain the rough positioning center (x _oc , y _oc ) and rough positioning circle radius R _c of the circular target, where (x _oc , y _oc ) is the coordinates of the center of the circle in the image coordinate system,

Step 2: Perform pixel-level edge positioning on the circular target in the image, the specific steps are as follows:

Step 2.1: According to the center of the rough positioning circle and the radius of the rough positioning circle, extract a square area from the image, called the "region of interest". The center point of the area, the rough positioning diameter of the circular target plus 2 to 6 pixels as the side length of the square area,

Step 2.2: Use the canny operator to perform pixel-level edge detection on the region of interest extracted in the image to obtain the pixel-level edge of the circular target, and then perform boundary tracking on the pixel-level edge points of the circular target extracted by the canny operator , get the coordinates of each pixel-level edge point, and store the coordinates of each pixel-level edge point in clockwise order,

Step 3: Perform sub-pixel edge positioning on the circular target in the image, the specific steps are as follows:

Step 3.1: Divide the edge pixels of the circular target into regions, and the division method is as follows: the x-axis is the straight line with the rough positioning center of the circular target as the origin, the positive direction horizontally to the right, and the unit length in pixels ; It is called the y-axis with the coarse positioning center of the circular target as the origin, the vertical upward as the positive direction, and the length in units of pixels as the y-axis; the coarse positioning center as the rotation point, and the angle of counterclockwise rotation from the positive direction of the x-axis is θ, according to the size of θ, the edge pixels of the circular target are divided into two parts, one part is θ∈[0°, 45°]∪[135°, 225°]∪[315°, 360°] corresponding to The edge pixels are called the edge pixels of the first area, and the other part is the edge pixels corresponding to θ∈(45°, 135°)∪(225°, 315°), which are called the edge pixels of the second area. The line connecting the center of the coarse positioning circle and the edge pixel point is the gradient direction of the edge pixel point,

Step 3.2: Calculate the adjacent points of the edge pixel point along the gradient direction, (x _p , y _p ) is the coordinate value of the edge point P in units of pixels, the specific method is as follows: if the edge pixel point P(x _p , y _p ) In the first area, take the intersection point A ₁ (x _a1 , y _a1 ) of the gradient direction line of the pixel point and the line x=x _p +1 as the first adjacent point on the right side of the edge pixel point along the gradient direction; Take the intersection point B ₁ (x _b1 , y _b1 ) of the gradient direction straight line of the pixel point and the straight line x=x _p +2 as the second adjacent point on the right side of the edge pixel point along the gradient direction; take the gradient direction straight line of the pixel point The point of intersection F(x _f , y _f ) with the straight line x=x _p +3 is the third adjacent point on the right side of the edge pixel along the gradient direction; take the gradient direction straight line of the pixel point and the straight line x=x _p -1 Intersection point C ₁ (x _c1 , y _c1 ) is the first adjacent point on the left side of the edge pixel along the gradient direction; take the intersection point D ₁ (x _d1 , _y _d1 ) is the second adjacent point on the left side of the edge pixel along the gradient direction, and the intersection point E ₁ (x _e1 , y _e1 ) of the gradient direction line of the pixel point and the line x=x _p -3 is the edge of the edge pixel point The third adjacent point on the left side of the gradient direction; if the edge pixel point P(x _p , y _p ) is in the second area, _take the intersection point A ₂ (x _a2 , y _a2 ) is the first adjacent point above the edge pixel along the gradient direction; take the intersection point B ₂ (x _b2 , y _b2 ) of the gradient direction line of the pixel point and the line y=y _p +2 as the edge pixel The second point next to the upper side of the point along the gradient direction; take the intersection point F ₂ (x _f2 , y _f2 ) of the gradient direction line of the pixel point and the straight line y=y _p +3 as the third point of the upper side of the edge pixel point along the gradient direction Close to the point; take the intersection point C ₂ (x _c2 , y _c2 ) of the gradient direction straight line of the pixel point and the straight line y=y _p -1 as the first adjacent point of the edge pixel point along the gradient direction; take the pixel point The intersection point D ₂ (x _d2 , y _d2 ) of the gradient direction straight line and the straight line y=y _p -2 is the second adjacent point on the upper edge of the edge pixel along the gradient direction; take the gradient direction straight line of the pixel point and the straight line y=y The intersection point E ₂ (x _e2 , y _e2 ) of _p -3 is the third immediate point above the edge pixel along the gradient direction,

Step 3.3: Use the method of linear grayscale interpolation to obtain the grayscale value of the adjacent points of the edge pixel point P(x _p , y _p ) along the gradient direction in units of pixels, where (x _p , y _p ) is The coordinate value of the edge pixel point P in pixels, use f(x, y) to represent the gray value of the pixel point whose coordinates are (x, y), use the symbol [] to represent the integer part, the above is used for The linear gray interpolation method to obtain the gray value of adjacent points is as follows:

(1) If the edge pixel point P(x _p , y _p ) is in the first area,

Gray value of adjacent point A ₁ (x _a1 , y _a1 )

f(x _a1 , y _a1 )=(1-λ)*f(x _a1 , [y _a1 ])+λ*f(x _a1 , [y _a1 ]+1)

λ=y _a1 -[y _a1 ],

Gray value of adjacent point B ₁ (x _b1 , y _b1 )

f(x _b1 , y _b1 )=(1-λ)*f(x _b1 ,[y _b1 ])+λ*f(x _b1 ,[y _b1 ]+1)

λ=y _b1 -[y _b1 ]

Gray value of adjacent point F ₁ (x _f1 , y _f1 )

f(x _f1 , y _f1 )=(1-λ)*f(x _f1 , [y _f1 ])+λ*f(x _f1 , [y _f1 ]+1)

λ=y _f1 -[y _f1 ]

Gray value of adjacent point C ₁ (x _c1 , y _c1 )

f(x _c1 , y _c1 )=(1-λ)*f(x _c1 ,[y _c1 ])+λ*f(x _c1 ,[y _c1 ]+1)

λ=y _c1 -[y _c1 ]

Gray value of adjacent point D ₁ (x _d , y _d )

f(x _d1 , y _d1 )=(1-λ)*f(x _d1 , [y _d1 ])+λ*f(x _d1 , [y _d1 ]+1)

λ=y _d1 -[y _d1 ]

Gray value of adjacent point E ₁ (x _e1 , y _e1 )

f(x _e1 , y _e1 )=(1-λ)*f(x _e1 ,[y _e1 ])+λ*f(x _e1 ,[y _e1 ]+1)

λ=y _e1 -[y _e1 ]

(2) If the edge pixel point P(x _p , y _p ) is in the second area,

Gray value of adjacent point A ₂ (x _a2 , y _a2 )

f(x _a2 , y _a2 )=(1-λ)*f([x _a2 ], y _a2 )+λ*f([x _a2 ]+1, y _a2 )

λ=x _a2 -[x _a2 ]

Gray value of adjacent point B ₂ (x _b2 , y _b2 )

f(x _b2 , y _b2 )=(1-λ)*f([x _b2 ], y _b2 )+λ*f([x _b2 ]+1, y _b2 )

λ=x _b2 -[x _b2 ]

Gray value of adjacent point F ₂ (x _f2 , y _f2 )

f(x _f2 , y _f2 )=(1-λ)*f([x _f2 ], y _f2 )+λ*f([x _f2 ]+1, y _f2 )

λ=x _f2 -[x _f2 ]

Gray value of adjacent point C ₂ (x _c2 , y _c2 )

f(x _c2 , y _c2 )=(1-λ)*f([x _c2 ], y _c2 )+λ*f([x _c2 ]+1, y _c2 )

λ=x _c2 -[x _c2 ]

Gray value of adjacent point D ₂ (x _d2 , y _d2 )

f(x _d2 , y _d2 )=(1-λ)*f([x _d2 ], y _d2 )+λ*f([x _d2 ]+1, y _d2 )

λ=x _d2 -[x _d2 ]

Gray value of adjacent point E ₂ (x _e2 , y _e2 )

f(x _e2 , y _e2 )=(1-λ)*f([x _e2 ], y _e2 )+λ*f([x _e2 ]+1, y _e2 )

λ=x _e2 -[x _e2 ]

Step 3.4: According to the gray value of the adjacent points of the edge point along the gradient direction obtained in step 3.3 and the gray value of the edge point P(x _p , y _p ) itself, select the forward difference and backward direction of the gray value The average value of the difference is taken as the corresponding gray level difference between the edge point and its adjacent points along the gradient direction, where (x _p , y _p ) is the coordinate value of the edge pixel point P in pixels, using f(x, y) represents the gray value of the pixel whose coordinates are (x, y), and the method for obtaining the gray difference is as follows:

(1) If the edge pixel point P(x _p , y _p ) is in the first area,

The gray level difference corresponding to the adjacent point A ₁ (x _a1 , y _a1 )

{f f}_{a a 11} = = \frac{| | f f (({x x}_{p p},, {y the y}_{p p})) - - f f (({x x}_{a a 11},, {y the y}_{a a 11})) | |}{22} + + \frac{| | f f (({x x}_{a a 11},, {y the y}_{a a 11})) - - f f (({x x}_{b b 11},, {y the y}_{b b 11})) | |}{22}

The gray level difference corresponding to the adjacent point B ₁ (x _b1 , y _b1 )

{f f}_{b b 11} = = \frac{| | f f (({x x}_{a a 11},, {y the y}_{a a 11})) - - f f (({x x}_{b b 11},, {y the y}_{b b 11})) | |}{22} + + \frac{| | f f (({x x}_{b b 11},, {y the y}_{b b 11})) - - f f (({x x}_{f f 11},, {y the y}_{f f 11})) | |}{22}

The gray level difference corresponding to the edge point P(x _p , y _p )

{f f}_{p p 11} = = \frac{| | f f (({x x}_{c c 11},, {y the y}_{c c 11})) - - f f (({x x}_{p p},, {y the y}_{p p})) | |}{22} + + \frac{| | f f (({x x}_{p p},, {y the y}_{p p})) - - f f (({x x}_{a a 11},, {y the y}_{a a 11})) | |}{22}

The gray level difference corresponding to the adjacent point C ₁ (x _c1 , y _c1 )

{f f}_{c c 11} = = \frac{| | f f (({x x}_{d d 11},, {y the y}_{d d 11})) - - f f (({x x}_{c c 11},, {y the y}_{c c 11})) | |}{22} + + \frac{| | f f (({c c}_{c c 11},, {y the y}_{c c 11})) - - f f (({x x}_{p p},, {y the y}_{p p})) | |}{22}

The gray level difference corresponding to the adjacent point D ₁ (x _d , y _d )

{f f}_{d d 11} = = \frac{| | f f (({x x}_{e e 11},, {y the y}_{e e 11})) - - f f (({x x}_{d d 11},, {y the y}_{d d 11})) | |}{22} + + \frac{| | f f (({x x}_{d d 11},, {y the y}_{d d 11})) - - f f (({x x}_{c c 11},, {y the y}_{c c 11})) | |}{22}

(2) If the edge pixel point P(x _p , y _p ) is in the second area,

The gray level difference corresponding to the adjacent point A ₂ (x _a2 , y _a2 )

{f f}_{a a 22} = = \frac{| | f f (({x x}_{p p},, {y the y}_{p p})) - - f f (({x x}_{a a 22},, {y the y}_{a a 22})) | |}{22} + + \frac{| | f f (({x x}_{a a 22},, {y the y}_{a a 22})) - - f f (({x x}_{b b 22},, {y the y}_{b b 22})) | |}{22}

The gray level difference corresponding to the adjacent point B ₂ (x _b2 , y _b2 )

{f f}_{b b 22} = = \frac{| | f f (({x x}_{a a 22},, {y the y}_{a a 22})) - - f f (({x x}_{b b 22},, {y the y}_{b b 22})) | |}{22} + + \frac{| | f f (({x x}_{b b 22},, {y the y}_{b b 22})) - - f f (({x x}_{f f 22},, {y the y}_{f f 22})) | |}{22}

The gray level difference corresponding to the edge point P(x _p , y _p )

{f f}_{p p 22} = = \frac{| | f f (({x x}_{c c 22},, {y the y}_{c c 22})) - - f f (({x x}_{p p},, {y the y}_{p p})) | |}{22} + + \frac{| | f f (({x x}_{p p},, {y the y}_{p p})) - - f f (({x x}_{a a 22},, {y the y}_{a a 22})) | |}{22}

The gray level difference corresponding to the adjacent point C ₂ (x _c2 , y _c2 )

{f f}_{c c 22} = = \frac{| | f f (({x x}_{d d 22},, {y the y}_{d d 22})) - - f f (({x x}_{c c 22},, {y the y}_{c c 22})) | |}{22} + + \frac{| | f f (({x x}_{c c 22},, {y the y}_{c c 22})) - - f f (({x x}_{p p},, {y the y}_{p p})) | |}{22}

The gray level difference corresponding to the adjacent point D ₂ (x _d2 , y _d2 )

{f f}_{d d 22} = = \frac{| | f f (({x x}_{e e 22},, {y the y}_{e e 22})) - - f f (({x x}_{d d 22},, {y the y}_{d d 22})) | |}{22} + + \frac{| | f f (({x x}_{d d 22},, {y the y}_{d d 22})) - - f f (({x x}_{c c 22},, {y the y}_{c c 22})) | |}{22}

Step 3.5: Calculate the distance difference δ between the sub-pixel edge point P′ and the pixel edge point P according to the corresponding gray level difference between the edge point and its adjacent points along the gradient direction obtained in step 3.4.

If one of the grayscale differences corresponding to the edge point and its adjacent points along the gradient direction is 0, then let δ=0, and store the count_zero of the number of occurrences of δ=0,

Step 3.6: Find the coordinates of the sub-pixel edge point P′, the specific method is: according to the distance difference δ between the sub-pixel edge point P′ and the pixel-level edge point P obtained in step 3.5 and the gradient direction straight line of the pixel-level edge point P The slope k _gradient , respectively calculate the coordinate difference δ _x , δ _y between point P and point P′ in the x direction and y direction, the formula is as follows:

{δ δ}_{x x} = = \frac{δ δ}{\sqrt{{k k}_{grad grad ient ient}^{22} + + 11}}

{δ δ}_{y the y} = = \frac{{k k}_{grad grad ient ient} * * δ δ}{\sqrt{{k k}_{grad grad ient ient}^{22} + + 11}}

Among them, the slope k _gradient of the gradient direction line of the pixel-level edge point P can be obtained from the coordinates (x _p , y _p ) of the pixel edge point and the coordinates (x _oc , y _oc ) of the coarse positioning circle, namely:

{k k}_{grad grad ient ient} = = \frac{{y the y}_{oc oc} - - {y the y}_{p p}}{{x x}_{oc oc} - - {x x}_{p p}}

Then, for a pixel point P(x _p , y _p ), its corresponding sub-pixel edge point is P′(x _p +δ _x , y _p +δ _y ),

Step 3.7: For all edge pixel points P(x _p , y _p ), according to steps 3.2 to 3.6, obtain the corresponding sub-pixel edge points as P′(x _p +δ _x , y _p +δ _y ), where (x _p +δ _x , y _p +δ _y ) is the coordinate of the sub-pixel edge point P′,

Step 4: Filter the sub-pixel edge points, and use the curvature filtering method and the mean value filtering method to process the "isolated points" and noises that appear at the sub-pixel edge points, respectively.

Step 4.1: For the "isolated point" that appears on the edge of the sub-pixel, use the method of curvature filtering to process it. The specific steps are as follows:

(1) Calculate the curvature of all sub-pixel edge points, the method is as follows: in the sequentially stored sub-pixel edge points,

is the coordinate of the sub-pixel edge point P′, is the coordinates of the immediately preceding point of point P′,

k k = = \frac{11}{r r} = = \frac{11}{\sqrt{(({(({x x}_{00} - - {x x}_{{p p}_{22}^{' '}}))}^{22} + + {(({y the y}_{00} - - {y the y}_{{p p}_{22}^{' '}}))}^{22}))}}

Among them, (x ₀ , y ₀ ) is the center of the circle formed by P', the point immediately before P', and the point immediately after P',

{x x}_{00} = = \frac{a a - - b b + + c c}{d d},, {y the y}_{00} = = \frac{e e - - f f + + g g}{- - d d},,

a a = = (({x x}_{{p p}_{11}^{' '}} + + {x x}_{{p p}_{22}^{' '}})) (({x x}_{{p p}_{22}^{' '}} - - {x x}_{{p p}_{11}^{' '}})) (({y the y}_{{p p}_{33}^{' '}} - - {y the y}_{{p p}_{22}^{' '}}))

b b = = (({x x}_{{p p}_{22}^{' '}} + + {x x}_{{p p}_{33}^{' '}})) (({x x}_{{p p}_{33}^{' '}} - - {x x}_{{p p}_{22}^{' '}})) (({y the y}_{{p p}_{22}^{' '}} - - {y the y}_{{p p}_{11}^{' '}}))

c c = = (({y the y}_{{p p}_{11}^{' '}} - - {y the y}_{{p p}_{33}^{' '}})) (({y the y}_{{p p}_{22}^{' '}} - - {y the y}_{{p p}_{11}^{' '}})) (({y the y}_{{p p}_{33}^{' '}} - - {y the y}_{{p p}_{22}^{' '}}))

d d = = 22 [[(({x x}_{{p p}_{22}^{' '}} - - {x x}_{{p p}_{11}^{' '}})) (({y the y}_{{p p}_{33}^{' '}} - - {y the y}_{{p p}_{22}^{' '}})) - - (({x x}_{{p p}_{33}^{' '}} - - {x x}_{{p p}_{22}^{' '}})) (({y the y}_{{p p}_{22}^{' '}} - - {y the y}_{{p p}_{11}^{' '}}))]]

e e = = (({y the y}_{{p p}_{11}^{' '}} + + {y the y}_{{p p}_{22}^{' '}})) (({y the y}_{{p p}_{22}^{' '}} - - {y the y}_{{p p}_{11}^{' '}})) (({x x}_{{p p}_{33}^{' '}} - - {x x}_{{p p}_{22}^{' '}}))

f f = = (({y the y}_{{p p}_{22}^{' '}} + + {y the y}_{{p p}_{33}^{' '}})) (({y the y}_{{p p}_{33}^{' '}} - - {y the y}_{{p p}_{22}^{' '}})) (({x x}_{{p p}_{22}^{' '}} - - {x x}_{{p p}_{11}^{' '}}))

g g = = (({x x}_{{p p}_{11}^{' '}} - - {x x}_{{p p}_{33}^{' '}})) (({x x}_{{p p}_{22}^{' '}} - - {x x}_{{p p}_{11}^{' '}})) (({x x}_{{p p}_{33}^{' '}} - - {x x}_{{p p}_{22}^{' '}}))

(2) Filter all sub-pixel edge points according to the curvature of the sub-pixel edge points, and the filtering criteria are as follows:

First, arrange the curvature of each sub-pixel edge point in descending order, and take the curvature threshold as the 3*n element of the curvature value of the sequence from high to low, where n is the number of isolated points, and its value is the value in step 3.5 Count_zero the number of occurrences of δ=0 during sub-pixel edge positioning; then use the curvature threshold to segment the curvature of the sub-pixel edge point, if the curvature of the sub-pixel edge point P is greater than the curvature threshold, and the curvature of P is greater than the curvature of its immediate front and rear points , the point P is considered to be an isolated point, and it is filtered out,

Step 4.2: Filter all sub-pixel edge points after filtering out "isolated points" with the mean value method, the specific method is: select the horizontal and vertical coordinates of the two points immediately before and after the sub-pixel edge point and the sub-pixel edge point itself The mean value of the horizontal and vertical coordinates is used as the new horizontal and vertical coordinates of the sub-pixel edge point,

Step 5: Use the least square method to fit the circle to the filtered sub-pixel edge points, and finally get the center and radius of the circular target.

2. The method for locating the center of a circular target according to claim 1, characterized in that:

The method for roughly positioning the center and radius of the circular target in the image as described in step 1 is:

Step 1.1: Carry out noise removal and threshold segmentation on the image to obtain a binary image with a gray value of 255 or 0 for each pixel in the image.

Step 1.2: Perform boundary extraction and boundary tracking on the circular target in the binary image,

Step 1.3: Fit the circle with the centroid method of the tracked boundary point, obtain the rough positioning circle center and the rough positioning circle radius of the circular target, and store the coarse positioning circle center coordinates (x _oc , y _oc ) and the rough positioning circle Radius R _c .