CN101334263A - Circular target circular center positioning method - Google Patents

Abstract

The invention discloses a circle center locating method of a circular target, mainly relating to a great deal of target identification and target location. Firstly, rough circle center location is carried out to an image by using a simple contour centroid method, a key square region is extracted as a region of interest according to the information of a rough location circle center and a rough location radius, and pixel-level edge location is carried out to the circular target in the region of interest by a canny operator; then sub-pixel location is carried out to the circular target according to the geometric features and the gray information of the circle, therefore, a precise coordinate of sub-pixel edge points is obtained; after that, a curvature filtering method and an average filtering method are respectively used for filtering 'isolated points' and noise occurring in the sub-pixel edge points; finally, a least squares method is utilized to fit a circle to the filtered sub-pixel edge points so as to obtain the final circular center and radius. The method not only effectively improves the precision of circle center location, but also improves the robustness thereof, thus further improving the measurement precision of a measurement system and perfecting the stability of the measurement system.

Description

The circle center locating method of circular target

Technical field

The present invention relates to a kind of sub-pixel edge location of circular target and the method for sub-pixel edge point filtering, relate in particular to a kind of circle center locating method of circular target.

Background technology

In numerous graphical analyses and disposal system, obtain the center of circle quickly and accurately and be widely used on Target Recognition and the location, as in the system of three-dimensionalreconstruction, obtaining the scaling board center of circle accurately has fundamental influence to whole three-dimensionalreconstruction precision.The common algorithm of asking for the center of circle at present has: centroid method, Hough converter technique, Gauss's surface fitting method, least square curve fitting method.It is more even that centroid method requires gradation of image to distribute, otherwise can produce than mistake; The Hough converter technique need be to each frontier point pointwise ballot, record, and computing time is longer, and it is bigger to take calculator memory, is restricted in actual applications; Gauss's surface fitting is that the principle of utilizing the intensity profile of circular light spot to be similar to Gauss model is carried out surface fitting, its bearing accuracy height, if bigger but handle the area of a circle, then operand is very big, and computing time is long, is not suitable in the real-time optical detection system; The least square curve fitting method is the most common method, because it realizes simple, calculated amount is little, be used widely in practice, but often all be to utilize the pixel edge point of asking for to carry out the least square fitting center of circle at present in the least square curve fitting method setting circle, make its fitting precision be subjected to certain restriction, thus sub-pixel edge circular in the accurate image can be obtained, for the fitting precision important influence that improves the center of circle.

Common sub-pixel edge location algorithm has at present: parameter fitting method, the method for square, method of interpolation and the method for utilizing other geometric properties of image.Document " Subpixel edge localization and the interpolation of still images " (Jensen K, Anastassiou, IEEE Trans on PAMI, 1995,4 (3): 285～295) be method of interpolation, be that the pixel edge location algorithm is used for the interpolation enlarged image, under real image intensity profile rule and the corresponding to prerequisite of interpolation rule, make the location, edge reach sub-pixel precision; Document " utilize tangential direction information detect sub-pixel edge " (Zhang Yujin, Fu is notable, pattern-recognition and artificial intelligence, 1997,10 (1): 69～74) mainly utilized the information of tangential direction, the target of known form is carried out the sub-pixel edge location, and this method is too dependent on the center of circle and the radius information of coarse positioning, does not make full use of the gray value information of image; Document " 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) middle method realization edge sub-pixel location with square.Method application at present based on square is more extensive, but owing to will use a plurality of 7 * 7 templates, therefore is not suitable for the real-time processing of large-scale bitmap; Document " the fast sub-picture element edge detection method of image " (Liu Lishuan, open small pot with a handle and a spout for boiling water or herbal medicine, photoelectron laser, 2005,16 (8): 993～996) introduced the principle of sub-pixel edge parameter fitting, it adopts is that the method for Sobel operator and parameter fitting realizes the sub-pixel edge location, wherein there are two problems, the first, the Sobel operator all can not be realized whole pixel edge location to most images, causes sub-pixel positioning inaccurate; The second, the method that the document provides can not be guaranteed the edge precision of quafric curve type.

Summary of the invention

At existing in prior technology shortcoming and restriction, the object of the present invention is to provide a kind of circle center locating method that can improve the circular target of measuring system precision.

The present invention designs and a kind ofly at first the circular target in the image is carried out center of circle coarse positioning, after obtaining the pixel edge point of circular target with the canny operator, take all factors into consideration round geometric properties and the edge gray distribution features is asked for the sub-pixel edge point coordinate, then accurate sub-pixel edge point is carried out filtering, try to achieve the method in the center of circle of circular target at last with the filtered sub-pixel edge point of least square fitting.The present invention adopts following technical scheme:

A kind of circle center locating method of circular target, its key step is:

The 1st step: the center of circle and radius to circular target in the image carry out coarse positioning, obtain the coarse positioning center of circle (x of circular target _Oc, y _Oc) and coarse positioning radius of circle R _c, wherein, (x _Oc, y _Oc) be the coordinate in the center of circle under the image coordinate system, concrete grammar is as follows:

The 1.1st step: image removed make an uproar, Threshold Segmentation, obtain each pixel gray-scale value in the image and be 255 or 0 bianry image;

The 1.2nd step: the circular target in the bianry image is carried out Boundary Extraction and border tracking;

The 1.3rd step: round to the frontier point of following the tracks of with the match of contour centroid method, obtain the coarse positioning center of circle and the coarse positioning radius of circle of this circular target, and store coarse positioning central coordinate of circle (x _Oc, y _Oc) and coarse positioning radius of circle R _c

The 2nd step: the circular target in the image is carried out the pixel edge location, and concrete steps are as follows:

The 2.1st step: according to the coarse positioning center of circle and coarse positioning radius of circle, from image, extract a square area, be called " area-of-interest ", the extracting method of this square area is: with the coarse positioning center of circle of circular target in the image central point as square area, add the length of side of 2～6 pixels as square area with the coarse positioning diameter of circular target;

The 2.2nd step: with the canny operator area-of-interest that is extracted in the image is carried out pixel edge and detect, obtain the pixel edge of circular target, again the pixel edge point of the circular target of canny operator extraction being carried out the border follows the tracks of, obtain the coordinate of each pixel edge point, and press the coordinate of clockwise each pixel edge point of sequential storage;

The 3rd step: the circular target in the image is carried out the sub-pixel edge location, and concrete steps are as follows:

The 3.1st step: the edge pixel point of circular target is carried out area dividing, and division methods is as follows: claim the coarse positioning center of circle with circular target to be initial point, to be positive dirction, to be that the straight line of unit length is the x axle with the pixel with level to the right; Title with the coarse positioning center of circle of circular target be initial point, to be positive dirction vertically upward, to be that the straight line of unit length is the y axle with the pixel; With the coarse positioning center of circle is the point of rotation, the angle that is rotated counterclockwise by x axle positive dirction is θ, edge pixel point according to big young pathbreaker's circular target of θ is divided into two parts, a part is [0 ° of θ ∈, 45 °] [135 ° of ∪, 225 °] [315 ° of ∪, 360 °] pairing edge pixel point, the edge pixel point that is called the 1st zone, another part is the corresponding edge pixel point of θ ∈ (45 °, 135 °) ∪ (225 °, 315 °), be called the edge pixel point in the 2nd zone, the line of the coarse positioning center of circle and edge pixel point is the gradient direction of this edge pixel point;

The 3.2nd step: ask for the neighbor point of edge pixel point, (x along gradient direction _p, y _p) for being the coordinate figure of the marginal point P of unit with the pixel, concrete grammar is as follows: if edge pixel point P (x _p, y _p) in the 1st zone, get the gradient direction straight line and the straight line x=x of this pixel _p+ 1 intersection point A ₁(x _A1, y _A1) be the right first neighbor point of this edge pixel point along gradient direction; Get the gradient direction straight line and the straight line x=x of this pixel _p+ 2 intersection points B ₁(x _B1, y _B1) be the right second neighbor point of this edge pixel point along gradient direction; Get the gradient direction straight line and the straight line x=x of this pixel _p+ 3 intersection point F (x _f, y _f) be the right three neighbor point of this edge pixel point along gradient direction; Get the gradient direction straight line and the straight line x=x of this pixel _p-1 intersection point C ₁(x _C1, y _C1) be the left side first neighbor point of this edge pixel point along gradient direction; Get the gradient direction straight line and the straight line x=x of this pixel _p-2 intersection point D ₁(x _D1, y _D1) be the left side second neighbor point of this edge pixel point along gradient direction, get the gradient direction straight line and the straight line x=x of this pixel _p-3 intersection point E ₁(x _E1, y _E1) be the left side three neighbor point of this edge pixel point along gradient direction; If edge pixel point P (x _p, y _p) in the 2nd zone, get the gradient direction straight line and the straight line y=y of this pixel _p+ 1 intersection point A ₂(x _A2, y _A2) be the top first next-door neighbour point of this edge pixel point along gradient direction; Get the gradient direction straight line and the straight line y=y of this pixel _p+ 2 intersection points B ₂(x _B2, y _B2) be the top second next-door neighbour point of this edge pixel point along gradient direction; Get the gradient direction straight line and the straight line y=y of this pixel _p+ 3 intersection point F ₂(x _F2, y _F2) be the top three next-door neighbour point of this edge pixel point along gradient direction; Get the gradient direction straight line and the straight line y=y of this pixel _p-1 intersection point C ₂(x _C2, y _C2) be the bottom first next-door neighbour point of this edge pixel point along gradient direction; Get the gradient direction straight line and the straight line y=y of this pixel _p-2 intersection point D ₂(x _D2, y _D2) be the top second next-door neighbour point of this edge pixel point along gradient direction; Get the gradient direction straight line and the straight line y=y of this pixel _p-3 intersection point E ₂(x _E2, y _E2) be the top three next-door neighbour point of this edge pixel point along gradient direction;

The 3.3rd step: adopting the method for linear gray-level interpolation to ask for the pixel is the edge pixel point P (x of unit _p, y _p) along the gray-scale value of the neighbor point of gradient direction, wherein, (x _p, y _p) for being the coordinate figure of the edge pixel point P of unit with the pixel, use f (x, y) denotation coordination be (use symbol [] is represented the round numbers part for x, gray values of pixel points y), and the described linear gray-level interpolation method that is used to obtain the neighbor point gray-scale value is as follows:

(1) if edge pixel point P (x _p, y _p) in the 1st zone,

Neighbor point A ₁(x _A1, y _A1) gray-scale value

f(x _a1，y _a1)＝(1-λ)*f(x _a1，[y _a1])+λ*f(x _a1，[y _a1]+1)

λ＝y _a1-[y _a1]，

Neighbor point B ₁(x _B1, y _B1) gray-scale value

f(x _b1，y _b1)＝(1-λ)*f(x _b1，[y _b1])+λ*f(x _b1，[y _b1]+1)

λ＝y _b1-[y _b1]

Neighbor point F ₁(x _F1, y _F1) gray-scale value

f(x _f1，y _f1)＝(1-λ)*f(x _f1，[y _f1])+λ*f(x _f1，[y _f1]+1)

λ＝y _f1-[y _f1]

Neighbor point C ₁(x _C1, y _C1) gray-scale value

f(x _c1，y _c1)＝(1-λ)*f(x _c1，[y _c1])+λ*f(x _c1，[y _c1]+1)

λ＝y _c1-[y _c1]

Neighbor point D ₁(x _d, y _d) gray-scale value

f(x _d1，y _d1)＝(1-λ)*f(x _d1，[y _d1])+λ*f(x _d1，[y _d1]+1)

λ＝y _d1-[y _d1]

Neighbor point E ₁(x _E1, y _E1) gray-scale value

f(x _e1，y _e1)＝(1-λ)*f(x _e1，[y _e1])+λ*f(x _e1，[y _e1]+1)

λ＝y _e1-[y _e1]

(2) if edge pixel point P (x _p, y _p) in the 2nd zone,

Neighbor point A ₂(x _A2, y _A2) gray-scale value

f(x _a2，y _a2)＝(1-λ)*f([x _a2]，y _a2)+λ*f([x _a2]+1，y _a2)

λ＝x _a2-[x _a2]

Neighbor point B ₂(x _B2, y _B2) gray-scale value

f(x _b2，y _b2)＝(1-λ)*f([x _b2]，y _b2)+λ*f([x _b2]+1，y _b2)

λ＝x _b2-[x _b2]

Neighbor point F ₂(x _F2, y _F2) gray-scale value

f(x _f2，y _f2)＝(1-λ)*f([x _f2]，y _f2)+λ*f([x _f2]+1，y _f2)

λ＝x _f2-[x _f2]

Neighbor point C ₂(x _C2, y _C2) gray-scale value

f(x _c2，y _c2)＝(1-λ)*f([x _c2]，y _c2)+λ*f([x _c2]+1，y _c2)

λ＝x _c2-[x _c2]

Neighbor point D ₂(x _D2, y _D2) gray-scale value

f(x _d2，y _d2)＝(1-λ)*f([x _d2]，y _d2)+λ*f([x _d2]+1，y _d2)

λ＝x _d2-[x _d2]

Neighbor point E ₂(x _E2, y _E2) gray-scale value

f(x _e2，y _e2)＝(1-λ)*f([x _e2]，y _e2)+λ*f([x _e2]+1，y _e2)

λ＝x _e2-[x _e2]

The 3.4th step: according to the 3.3rd step the marginal point of asking along the gray-scale value and the marginal point P (x of the neighbor point of gradient direction _p, y _p) gray-scale value itself, the mean value of choosing the forward difference of gray-scale value and backward difference is as marginal point and along the gray scale difference of the correspondence of the neighbor point of gradient direction, wherein, and (x _p, y _p) for being the coordinate figure of the edge pixel point P of unit with the pixel, use f (x, y) denotation coordination be (the described method of asking for gray scale difference is as follows for x, gray values of pixel points y):

(1) if edge pixel point P (x _p, y _p) in the 1st zone,

Neighbor point A ₁(x _A1, y _A1) corresponding gray scale difference

$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}$

Neighbor point B ₁(x _B1, y _B1) corresponding gray scale difference

$f_{b1}=\frac{|f(x_{a1},y_{a1})-f(x_{b1},y_{b1})|}{2}+\frac{|f(x_{b1},y_{b1})-f(x_{\mathrm{<figure-callout\; id="1"\; label="f"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">f</figure-callout>}1},y_{f1})|}{2}$

Marginal point P (x _p, y _p) corresponding gray scale difference

${\mathrm{<figure-callout\; id="1"\; label="f\; p"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">f</figure-callout>}}_p=\frac{|f(x_{\mathrm{<figure-callout\; id="1"\; label="c"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">c</figure-callout>}1},y_{c1})-f(x_p,y_p)|}{2}+\frac{|f(x_p,y_p)-f(x_{a1},y_{a1})|}{2}$

Neighbor point C ₁(x _C1, y _C1) corresponding gray scale difference

${\mathrm{<figure-callout\; id="1"\; label="f\; c"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">f</figure-callout>}}_c=\frac{|f({\mathrm{<figure-callout\; id="1"\; label="x\; d"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">x</figure-callout>}}_d,y_{d1})-f(x_{\mathrm{<figure-callout\; id="1"\; label="c"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">c</figure-callout>}1},y_{c1})|}{2}+\frac{|f(x_{\mathrm{<figure-callout\; id="1"\; label="c"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">c</figure-callout>}1},y_{c1})-f(x_p,y_p)|}{2}$

Neighbor point D ₁(x _d, y _d) corresponding gray scale difference

${\mathrm{<figure-callout\; id="1"\; label="f\; d"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">f</figure-callout>}}_d=\frac{|f(x_{\mathrm{<figure-callout\; id="1"\; label="e"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">e</figure-callout>}1},y_{e1})-f({\mathrm{<figure-callout\; id="1"\; label="x\; d"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">x</figure-callout>}}_d,y_{d1})|}{2}+\frac{|f({\mathrm{<figure-callout\; id="1"\; label="x\; d"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">x</figure-callout>}}_d,y_{d1})-f(x_{\mathrm{<figure-callout\; id="1"\; label="c"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">c</figure-callout>}1},y_{c1})|}{2}$

(2) if edge pixel point P (x _p, y _p) in the 2nd zone,

Neighbor point A ₂(x _A2, y _A2) corresponding gray scale difference

$f_{a2}=\frac{|f(x_p,y_p)-f(x_{a2},y_{a2})|}{2}+\frac{|f(x_{a2},y_{a2})-f(x_{\mathrm{<figure-callout\; id="2"\; label="b"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">b</figure-callout>}2},y_{b2})|}{2}$

Neighbor point B ₂(x _B2, y _B2) corresponding gray scale difference

${\mathrm{<figure-callout\; id="2"\; label="f\; b"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">f</figure-callout>}}_b=\frac{|f(x_{a2},y_{a2})-f(x_{\mathrm{<figure-callout\; id="2"\; label="b"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">b</figure-callout>}2},y_{b2})|}{2}+\frac{|f(x_{\mathrm{<figure-callout\; id="2"\; label="b"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">b</figure-callout>}2},y_{b2})-f({\mathrm{<figure-callout\; id="2"\; label="x\; f"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">x</figure-callout>}}_f,y_{f2})|}{2}$

Marginal point P (x _p, y _p) corresponding gray scale difference

${\mathrm{<figure-callout\; id="2"\; label="f\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">f</figure-callout>}}_p=\frac{|f(x_{\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">c</figure-callout>}2},y_{c2})-f(x_p,y_p)|}{2}+\frac{|f(x_p,y_p)-f(x_{a2},y_{a2})|}{2}$

Neighbor point C ₂(x _C2, y _C2) corresponding gray scale difference

${\mathrm{<figure-callout\; id="2"\; label="f\; c"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">f</figure-callout>}}_c=\frac{|f({\mathrm{<figure-callout\; id="2"\; label="x\; d"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">x</figure-callout>}}_d,y_{d2})-f(x_{\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">c</figure-callout>}2},y_{c2})|}{2}+\frac{|f(x_{\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">c</figure-callout>}2},y_{c2})-f(x_p,y_p)|}{2}$

Neighbor point D ₂(x _D2, y _D2) corresponding gray scale difference

${\mathrm{<figure-callout\; id="2"\; label="f\; d"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">f</figure-callout>}}_d=\frac{|f(x_{\mathrm{<figure-callout\; id="2"\; label="e"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">e</figure-callout>}2},y_{e2})-f({\mathrm{<figure-callout\; id="2"\; label="x\; d"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">x</figure-callout>}}_d,y_{d2})|}{2}+\frac{|f({\mathrm{<figure-callout\; id="2"\; label="x\; d"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">x</figure-callout>}}_d,y_{d2})-f(x_{\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">c</figure-callout>}2},y_{c2})|}{2}$

The 3.5th step: according to gained marginal point in the 3.4th step and along the gray scale difference of the correspondence of the neighbor point of gradient direction, the range difference δ that asks for sub-pixel edge point P ' and pixel edge point P is

If marginal point and along having to be 0 in the gray scale difference of the neighbor point correspondence of gradient direction then makes δ=0, and storage δ=0 number of times Count_zero that occurs;

The 3.6th step: ask for the coordinate of sub-pixel edge point P ', concrete grammar is: according to 3.5 go on foot the slope k of gradient direction straight line of the range difference δ of the sub-pixel edge point P ' that asks and pixel edge point P and pixel edge point P _Gradient, obtain P point and P ' some coordinate difference δ on x direction and y direction respectively _x, δ _y, formula is as follows:

${\mathrm{\&delta;}}_x=\frac{\mathrm{\&delta;}}{\sqrt{k_{\mathrm{grad}\mathrm{ient}}^2+1}}$

${\mathrm{\&delta;}}_y=\frac{k_{\mathrm{grad}\mathrm{ient}}*\mathrm{\&delta;}}{\sqrt{k_{\mathrm{grad}\mathrm{ient}}^2+1}}$

Wherein, the slope k of the gradient direction straight line of pixel edge point P _GradientCan be by the coordinate (x of this pixel edge point _p, y _p) and coarse positioning central coordinate of circle (x _Oc, y _Oc) obtain, that is:

$k_{\mathrm{grad}\mathrm{ient}}=\frac{y_{\mathrm{oc}}-y_p}{x_{\mathrm{oc}}-x_p}$

So, for pixel P (x _p, y _p), its corresponding sub-pixel edge point is P ' (x _p+ δ _x, y _p+ δ _y);

The 3.7th step: to all edge pixel point P (x _p, y _p), asking for corresponding sub-pixel edge point successively according to～the 3.6 step of the 3.2nd step is P ' (x _p+ δ _x, y _p+ δ _y), wherein, (x _p+ δ _x, y _p+ δ _y) be the coordinate of sub-pixel edge point P ';

The 4th step: the sub-pixel edge point is carried out filtering,, handle with the method for curvature filtering and the method for mean filter respectively at " isolated point " and the noise that sub-pixel edge point occurs;

The 4.1st step: at " isolated point " that sub-pixel edge occurs, handle with the method for curvature filtering, concrete steps are as follows:

(1) ask for the curvature of all sub-pixel edge points, method is as follows: in the sub-pixel edge point of sequential storage, (x ' _P2, y ' _P2) be the coordinate of sub-pixel edge point P ', (x ' _P1, y ' _P1) be any coordinate before P ' the next-door neighbour, (x ' _P3, y ' _P3) be any coordinate behind P ' the next-door neighbour, the curvature of sub-pixel edge point P ' is so

$k=\frac{1}{r}=\frac{1}{\sqrt{({(x_0-{\mathrm{<figure-callout\; id="2"\; label="x\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">x</figure-callout>}}_{p_{}^{}})}^2+{(y_0-{\mathrm{<figure-callout\; id="2"\; label="y\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">y</figure-callout>}}_{p_{}^{}})}^2)}}$

Wherein, (x ₀, y ₀) for before P ', P ' the next-door neighbour a bit, any forms the round center of circle behind the next-door neighbour of P ',

$x_0=\frac{a-b+c}{d},$ ${\mathrm{<figure-callout\; id="0"\; label="y"\; filenames="A20081002279300322.png,A20081002279300341.png"\; state="\{\{state\}\}">y</figure-callout>}}_0=\frac{e-f+g}{-d},$

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)

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)]

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)

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

(2) according to the curvature of sub-pixel edge point all sub-pixel edge points are carried out filtering, filter criteria is as follows:

At first the curvature of each sub-pixel edge point is carried out descending sort, get curvature threshold and be this sequence curvature value 3*n element from high to low, wherein, n is the number of isolated point, and its numerical value is the occurrence number Count_zero of δ=0 when sub-pixel edge is located in the 3.5th step; With curvature threshold the curvature of sub-pixel edge point is cut apart then, if the curvature of sub-pixel edge point P greater than curvature threshold, and the curvature of P thinks then that greater than the curvature of next-door neighbour's point before and after it P is an isolated point, gives filtering;

The 4.2nd step: all the sub-pixel edge points after the filtering " isolated point " are carried out filtering with Mean Method, and concrete grammar is: the average of horizontal stroke, ordinate of choosing 2 horizontal stroke, ordinate and the sub-pixel edge point itself that sub-pixel edge point front and back are close to is as new horizontal stroke, the ordinate of this sub-pixel edge point;

The 5th step: filtered sub-pixel edge point is justified with least square fitting, finally obtained the center of circle and the radius of circular target.

Compared with prior art, the present invention has following advantage:

(1) compares with traditional two dimensional image edge localization method, the method of the sub-pixel edge location that the present invention proposes is carried out based on the determined area-of-interest in the coarse positioning center of circle, and when match edge gray scale difference, Gauss's surface fitting is reduced to gaussian curve approximation along gradient direction, simplify calculating, improved the rapidity of algorithm.

(2) method of the sub-pixel edge location of the present invention's proposition is when considering the gradient direction of each edge pixel point, be not limited to the several special direction that traditional edge detection operator is determined, but considered the concrete gradient direction of each edge pixel point according to the feature of circle, and utilized linear gray-level interpolation to obtain the gray-scale value of edge pixel point, and then made the location of marginal point more accurate along the neighbor point of gradient direction.

(3) the present invention is directed in the sub-pixel edge point " isolated point " and the noise that occurs, proposed a kind of simple and practical filtering method, this method not only principle and algorithm is simple, and filter effect is fine.

(4) circle center locating method of the circular target of the present invention's proposition is compared with existent method, higher precision and better stable is arranged, and have very strong practicality.

Description of drawings

Fig. 1 is the process flow diagram of the circle center locating method concrete steps of circular target.

Fig. 2 is area-of-interest extraction figure.

Fig. 3 is that non-maximum value suppresses synoptic diagram, and wherein Fig. 3 (a) is four sector charts dividing according to gradient direction; 3 (b) are 3*3 neighborhood figure.

Fig. 4 (a) is the schematic diagram of the definition pixel of edge pixel point P when the 1st zone along the gradient direction neighbor point; Fig. 4 (b) is the schematic diagram of the definition pixel of edge pixel point P when the 2nd zone along the gradient direction neighbor point.

Fig. 5 (a) be when edge pixel point P in the 1st when zone, edge pixel point P schemes along the neighbor point definition of gradient direction; Fig. 5 (b) be when edge pixel point P in the 2nd when zone, edge pixel point P schemes along the neighbor point definition of gradient direction.

Fig. 6 is that the sub-pixel edge point is asked for schematic diagram.

Embodiment

Below in conjunction with accompanying drawing the specific embodiment of the present invention is further described.According to said method, in Windows operating system, realized the center of circle positioning action of circular target with the C++ programming by the VC++6.0 platform.

Using this method carries out the location, the center of circle of circular target and comprises that mainly center of circle coarse positioning, pixel edge location, sub-pixel edge location, the filtering of sub-pixel edge point, least square fitting try to achieve five operation stepss in the center of circle, the process flow diagram of concrete steps is used the concrete steps following (not specified unit all is unit with the pixel in the step) that this method is carried out the location, the center of circle of circular target as shown in Figure 1:

The 1st step: the center of circle and radius to the circular target in the image carry out coarse positioning, obtain the coarse positioning center of circle (x of circular target _Oc, y _Oc) and coarse positioning radius of circle R _c, wherein, (x _Oc, y _Oc) be the coordinate in the center of circle under the image coordinate system, concrete steps are as follows:

The 1.1st step: image removed make an uproar, Threshold Segmentation, obtain each pixel gray-scale value in the image and be 255 or 0 bianry image.Wherein, Threshold Segmentation adopts the iteration threshold method.The iteration threshold method can calculate proper segmentation threshold automatically, and its calculation procedure is as follows:

(1) select threshold value T, the average gray value of selecting image usually is as initial threshold;

(2), the average gray value of image is divided into two groups of R by initial threshold T ₁And R ₂

(3) calculate two groups of average gray value μ ₁And μ ₂

(4) reselect threshold value T, new T is defined as: T=(μ ₁+ μ ₂)/2;

Circulation (2) to (4) is up to two groups average gray value μ ₁And μ ₂No longer change, therefore obtain needed threshold value, according to resulting threshold value image is carried out binaryzation again.

The 1.2nd step: the circular target in the bianry image is carried out Boundary Extraction and border tracking.In bianry image, gray-scale value is that 255 pixel is circular target, and the method for circular target being carried out Boundary Extraction is as follows: if having among the former figure a bit for white, and 8 adjacent pixels all are to deceive, and then this gray values of pixel points are made as 0.Through after the Boundary Extraction, the target wheel profile among this figure is followed the tracks of and stored, can obtain the frontier point and the coordinate thereof of circular target.

The 1.3rd step: round to the frontier point of following the tracks of with the match of contour centroid method, obtain the coarse positioning center of circle and the coarse positioning radius of circle of this circular target, and store coarse positioning central coordinate of circle (x _Oc, y _Oc) and coarse positioning radius of circle R _cThe contour centroid method is a kind of the simplest sub-pixel positioning algorithm, and algorithm is as follows: the frontier point coordinate of establishing tracking for (i, j), the lock-on boundary point adds up to N, according to the contour centroid method, the central coordinate of circle of circular target is:

$\left\{\begin{array}{cc}{x}_{\mathrm{oc}}=\frac{\mathrm{\&Sigma;i}}{N}& \\ y_{\mathrm{oc}}=\frac{\mathrm{\&Sigma;j}}{N}& \end{array}\right.$

The radius of circular target for each frontier point to distance of center circle from mean value, that is:

$R_c=\frac{\mathrm{\&Sigma;}\sqrt{{(i-x_{\mathrm{oc}})}^2+{(j-y_{\mathrm{oc}})}^2}}{N}$

The 2nd step: the circular target in the image is carried out the pixel edge location, and concrete steps are as follows:

The 2.1st step: according to the coarse positioning center of circle and coarse positioning radius of circle, from image, extract a square area, be called " area-of-interest ", the extracting method of this square area is: with the coarse positioning center of circle of circular target in the image central point as square area, add the length of side of 2～6 pixels as square area with the coarse positioning diameter of circular target.As shown in Figure 2, the circular target among the figure in the circle representative image, square area is area-of-interest, and wherein, O is the coarse positioning center of circle, R _cBe the coarse positioning radius, d generally gets 2～6 pixels.Subsequent step is only handled at area-of-interest.

The 2.2nd step: with the canny operator area-of-interest that is extracted in the image is carried out pixel edge and detect, obtain the pixel edge of circular target, again the pixel edge point of the circular target of canny operator extraction being carried out the border follows the tracks of, obtain the coordinate of each pixel edge point, and press the coordinate of clockwise each pixel edge point of sequential storage.Wherein, the Canny rim detection is estimated signal to noise ratio (S/N ratio) and location product, belong to earlier level and smooth after the method for differentiate again, its key step and ultimate principle are as follows:

(1) use the Gaussian filter smoothed image, Gauss's smooth function is as follows:

$H(x,y)=e^{-\frac{x^2+{\mathrm{<figure-callout\; id="2"\; label="y"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">y</figure-callout>}}^2}{{2\mathrm{\&sigma;}}^2}}$

Wherein, σ is the standard deviation of Gaussian function, and x, y are respectively horizontal stroke, the ordinates of pixel in the image.Image is carried out gaussian filtering, promptly to image f _I(x, y) with Gaussian function H (x y) carries out convolution, obtain filtered smoothed image G (x, y):

G(x，y)＝f _I(x，y)*H(x，y)

(2) with the finite difference of single order partial derivative the assign to amplitude and the direction of compute gradient.First order difference convolution template is as follows:

${\mathrm{<figure-callout\; id="1"\; label="H"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">H</figure-callout>}}_1=\left|\begin{array}{cc}-1& -1\\ 1& 1\end{array}\right|$ ${\mathrm{<figure-callout\; id="2"\; label="H"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">H</figure-callout>}}_2=\left|\begin{array}{cc}1& -1\\ 1& -1\end{array}\right|$

(x y) carries out convolution with above two templates, and the gradient image that obtains on level and the vertical direction is respectively with the image G after level and smooth And then by

Can obtain the amplitude of gradient

And direction

The reckoning formula is as follows:

(3) (non-maxima suppression, NMS): Fig. 3 is that non-maximum value suppresses synoptic diagram gradient magnitude to be carried out non-maximum value inhibition.According to the direction of gradient, the label of four sectors is 0 to 3, and shown in Fig. 3 (a), four kinds of corresponding 3*3 neighborhood may be made up, and the 3*3 neighborhood is shown in Fig. 3 (b).On each pixel, with the center pixel P of neighborhood _CpCompare with two pixels, if P along gradient line _CpGrad big unlike two neighbor Grad along gradient line, then this is not a marginal point.

(4) usefulness dual threshold algorithm detects and is connected the edge.The dual threshold algorithm suppresses two threshold tau 1 of image effect and τ 2 to non-maximum value, and 2 τ, 1 ≈ τ 2, thereby can obtain two threshold value edge image N ₁[i, j] and N ₂[i, j].Because N ₂[i, j] uses high threshold to obtain, thereby contains false edge seldom, but interruption is arranged.The dual threshold method will be at N ₂In [i, j] edge is connected into profile, when arriving the end points of profile, this algorithm is just at N ₁Seek the edge that can be connected on the profile in the 8 neighborhood points of [i, j], like this, algorithm is constantly at N ₁Collect edge in [i, j], up to N ₂Till [i, j] couples together.

The 3rd step: the circular target in the image is carried out the sub-pixel edge location, and concrete steps are as follows:

The 3.1st step: the edge pixel point of circular target is carried out area dividing, and division methods is as follows: claim the coarse positioning center of circle with circular target to be initial point, to be positive dirction, to be that the straight line of unit length is the x axle with the pixel with level to the right; Title with the coarse positioning center of circle of circular target be initial point, to be positive dirction vertically upward, to be that the straight line of unit length is the y axle with the pixel; With the coarse positioning center of circle is the point of rotation, the angle that is rotated counterclockwise by x axle positive dirction is θ, edge pixel point according to big young pathbreaker's circular target of θ is divided into two parts, a part is [0 ° of θ ∈, 45 °] [135 ° of ∪, 225 °] [315 ° of ∪, 360 °] pairing edge pixel point, the edge pixel point that is called the 1st zone, another part is the corresponding edge pixel point of θ ∈ (45 °, 135 °) ∪ (225 °, 315 °), be called the edge pixel point in the 2nd zone, the line of the coarse positioning center of circle and edge pixel point is the gradient direction of this edge pixel point.

The 3.2nd step: ask for the neighbor point of edge pixel point along gradient direction, concrete principle is as follows: if edge pixel point P (x _p, y _p) in the 1st zone, Fig. 4 (a) is the schematic diagram of definition pixel along the gradient direction neighbor point, (x _p, y _p) for being the coordinate figure of the marginal point P of unit with the pixel, L _GBe the gradient direction straight line that P is ordered, straight line y ₁, y ₂, y ₃, y ₄Is four horizontal lines of unit for the P point is neighbouring with the pixel, respectively with straight line L _GIntersect at a B ' ₁, A ' ₁, C ' ₁, D ' ₁, straight line x ₁, x ₂, x ₃, x ₄Is four ordinates of unit for the P point is neighbouring with the pixel, respectively with straight line L _GIntersect at a D ₁, C ₁, A ₁, B ₁, choose B ' ₁, A ' ₁, C ' ₁, D ' ₁And { B ₁, A ₁, C ₁, D ₁Two groups of points are as the candidate neighbor point of marginal point along gradient direction since P o'clock in the 1st zone, gradient direction is based on the x direction, i.e. L _GThe absolute value of slope less than 1, therefore, candidate's neighbor point B ' ₁, A ' ₁, C ' ₁, D ' ₁The distance of ordering to P obviously is greater than candidate's neighbor point { B ₁, A ₁, C ₁, D ₁Arrive the distance that P is ordered, if select for use B ' ₁, A ' ₁, C ' ₁, D ' ₁As the neighbor point of P point, can cause the location, edge accurate inadequately along gradient direction, therefore, select { B for use ₁, A ₁, C ₁, D ₁As its neighbor point; If edge pixel point P (x _p, y _p) in the 2nd zone, Fig. 4 (b) is the schematic diagram of definition pixel along the gradient direction neighbor point, (x _p, y _p) for being the coordinate figure of the marginal point P of unit with the pixel, L _GBe the gradient direction straight line that P is ordered, straight line y ₁, y ₂, y ₃, y ₄Is four horizontal lines of unit for the P point is neighbouring with the pixel, respectively with straight line I _GIntersect at a B ₂, A ₂, C ₂, D ₂, straight line x ₁, x ₂, x ₃, x ₄Is four ordinates of unit for the P point is neighbouring with the pixel, respectively with straight line L _GIntersect at a D ' ₂, C ' ₂, A ' ₂, B ' ₂, choose B ' ₂, A ' ₂, C ' ₂, D ' ₂And { B ₂, A ₂, C ₂, D ₂Two groups of points are as the candidate neighbor point of marginal point along gradient direction.Because gradient direction was based on y direction, i.e. L in the 2nd zone in P o'clock _GThe absolute value of slope greater than 1, therefore, candidate's neighbor point B ' ₂, A ' ₂, C ' ₂, D ' ₂The distance of ordering to P obviously is greater than candidate's neighbor point { B ₂, A ₂, C ₂, D ₂Arrive the distance that P is ordered, if select for use B ' ₂, A ' ₂, C ' ₂, D ' ₂As the neighbor point of P point, can cause the location, edge accurate inadequately along gradient direction, therefore, select { B for use ₂, A ₂, C ₂, D ₂As its neighbor point.Therefore, the pixel of two zoness of different of dividing at the 3.1st step, its neighbor point along gradient direction is defined as follows: (x _p, y _p) for being the coordinate figure of the marginal point P of unit, if edge pixel point P (x with the pixel _p, y _p) in the 1st zone, Fig. 5 (a) is the neighbor point definition figure of pixel along gradient direction, L _GBe edge pixel point P (x _p, y _p) the gradient direction straight line, L ₁Represent straight line x=x _p-3, L ₂Represent straight line x=x _p-2, L ₃Represent straight line x=x _p-1, L ₄Represent straight line x=x _p+ 1, L ₅Represent straight line x=x _p+ 2, L ₆Represent straight line x=x _p+ 3.Get the gradient direction straight line and the straight line x=x of this pixel _p+ 1 intersection point A ₁(x _A1, y _A1) be the right first neighbor point of this edge pixel point along gradient direction; Get the gradient direction straight line and the straight line x=x of this pixel _p+ 2 intersection points B ₁(x _B1, y _B1) be the right second neighbor point of this edge pixel point along gradient direction; Get the gradient direction straight line and the straight line x=x of this pixel _p+ 3 intersection point F (x _f, y _f) be the right three neighbor point of this edge pixel point along gradient direction; Get the gradient direction straight line and the straight line x=x of this pixel _p-1 intersection point C ₁(x _C1, y _C1) be the left side first neighbor point of this edge pixel point along gradient direction; Get the gradient direction straight line and the straight line x=x of this pixel _p-2 intersection point D ₁(x _D1, y _D1) be the left side second neighbor point of this edge pixel point along gradient direction, get the gradient direction straight line and the straight line x=x of this pixel _p-3 intersection point E ₁(x _E1, y _E1) be the left side three neighbor point of this edge pixel point along gradient direction; If edge pixel point P (x _p, y _p) in the 2nd zone, Fig. 5 (b) is the neighbor point definition figure of pixel along gradient direction, L _GBe edge pixel point P (x _p, y _p) the gradient direction straight line, get the gradient direction straight line and the straight line y=y of this pixel _p+ 1 intersection point A ₂(x _A2, y _A2) be the top first next-door neighbour point of this edge pixel point along gradient direction; Get the gradient direction straight line and the straight line y=y of this pixel _p+ 2 intersection points B ₂(x _B2, y _B2) be the top second next-door neighbour point of this edge pixel point along gradient direction; Get the gradient direction straight line and the straight line y=y of this pixel _p+ 3 intersection point F ₂(x _F2, y _F2) be the top three next-door neighbour point of this edge pixel point along gradient direction; Get the gradient direction straight line and the straight line y=y of this pixel _p-1 intersection point C ₂(x _C2, y _C2) be the bottom first next-door neighbour point of this edge pixel point along gradient direction; Get the gradient direction straight line and the straight line y=y of this pixel _p-2 intersection point D ₂(x _D2, y _D2) be the top second next-door neighbour point of this edge pixel point along gradient direction; Get the gradient direction straight line and the straight line y=y of this pixel _p-3 intersection point E ₂(x _E2, y _E2) be the top three next-door neighbour point of this edge pixel point along gradient direction.

The 3.3rd step: adopting the method for linear gray-level interpolation to ask for the pixel is the edge pixel point P (x of unit _p, y _p) along the gray-scale value of the neighbor point of gradient direction, concrete grammar is: at first find out two the whole pixels nearest with the neighbor point along gradient direction of this pixel, then these two whole pixels are carried out linear gray-level interpolation, interpolation result is the gray-scale value along the neighbor point of gradient direction of this edge pixel point, wherein, (x _p, y _p) for being the coordinate figure of the edge pixel point P of unit with the pixel, use f (x, y) denotation coordination be (use symbol [] is represented the round numbers part for x, gray values of pixel points y), and the described linear gray-level interpolation method that is used to obtain the neighbor point gray-scale value is as follows:

(1) if edge pixel point P (x _p, y _p) in the 1st zone,

Neighbor point A ₁(x _A1, y _A1) gray-scale value

f(x _a1，y _a1)＝(1-λ)*f(x _a1，[y _a1])+λ*f(x _a1，[y _a1]+1)

λ＝y _a1-[y _a1]，

Neighbor point B ₁(x _B1, y _B1) gray-scale value

f(x _b1，y _b1)＝(1-λ)*f(x _b1，[y _b1])+λ*f(x _b1，[y _b1]+1)

λ＝y _b1-[y _b1]

Neighbor point F ₁(x _F1, y _F1) gray-scale value

f(x _f1，y _f1)＝(1-λ)*f(x _f1，[y _f1])+λ*f(x _f1，[y _f1]+1)

λ＝y _f1-[y _f1]

Neighbor point C ₁(x _C1, y _C1) gray-scale value

f(x _c1，y _c1)＝(1-λ)*f(x _c1，[y _c1])+λ*f(x _c1，[y _c1]+1)

λ＝y _c1-[y _c1]

Neighbor point D ₁(x _d, y _d) gray-scale value

f(x _d1，y _d1)＝(1-λ)*f(x _d1，[y _d1])+λ*f(x _d1，[y _d1]+1)

λ＝y _d1-[y _d1]

Neighbor point E ₁(x _E1, y _E1) gray-scale value

f(x _e1，y _e1)＝(1-λ)*f(x _e1，[y _e1])+λ*f(x _e1，[y _e1]+1)

λ＝y _e1-[y _e1]

(2) if edge pixel point P (x _p, y _p) in the 2nd zone,

Neighbor point A ₂(x _A2, y _A2) gray-scale value

f(x _a2，y _a2)＝(1-λ)*f([x _a2]，y _a2)+λ*f([x _a2]+1，y _a2)

λ＝x _a2-[x _a2]

Neighbor point B ₂(x _B2, y _B2) gray-scale value

f(x _b2，y _b2)＝(1-λ)*f([x _b2]，y _b2)+λ*f([x _b2]+1，y _b2)

λ＝x _b2-[x _b2]

Neighbor point F ₂(x _F2, y _F2) gray-scale value

f(x _f2，y _f2)＝(1-λ)*f([x _f2]，y _f2)+λ*f([x _f2]+1，y _f2)

λ＝x _f2-[x _f2]

Neighbor point C ₂(x _C2, y _C2) gray-scale value

f(x _c2，y _c2)＝(1-λ)*f([x _c2]，y _c2)+λ*f([x _c2]+1，y _c2)

λ＝x _c2-[x _c2]

Neighbor point D ₂(x _D2, y _D2) gray-scale value

f(x _d2，y _d2)＝(1-λ)*f([x _d2]，y _d2)+λ*f([x _d2]+1，y _d2)

λ＝x _d2-[x _d2]

Neighbor point E ₂(x _E2, y _E2) gray-scale value

f(x _e2，y _e2)＝(1-λ)*f([ _e2]，y _e2)+λ*f([ _e2]+1，y _e2)

λ＝x _e2-[x _e2]

The 3.4th step: according to the 3.3rd step the marginal point of asking along the gray-scale value and the marginal point P (x of the neighbor point of gradient direction _p, y _p) gray-scale value itself, the mean value of choosing the forward difference of gray-scale value and backward difference is as marginal point and along the gray scale difference of the correspondence of the neighbor point of gradient direction, wherein, and (x _p, y _p) for being the coordinate figure of the edge pixel point P of unit with the pixel, use f (x, y) denotation coordination be (the described method of asking for gray scale difference is as follows for x, gray values of pixel points y):

(1) if edge pixel point P (x _p, y _p) in the 1st zone,

Neighbor point A ₁(x _A1, y _A1) corresponding gray scale difference

$f_{a1}=\frac{|f(x_p,y_p)-f(x_{a1},y_{a1})|}{2}+\frac{|f(x_{a1},y_{a1})-f(x_{\mathrm{<figure-callout\; id="1"\; label="b"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">b</figure-callout>}1},y_{b1})|}{2}$

Neighbor point B ₁(x _B1, y _B1) corresponding gray scale difference

${\mathrm{<figure-callout\; id="1"\; label="f\; b"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">f</figure-callout>}}_b=\frac{|f(x_{a1},y_{a1})-f(x_{\mathrm{<figure-callout\; id="1"\; label="b"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">b</figure-callout>}1},y_{b1})|}{2}+\frac{|f(x_{\mathrm{<figure-callout\; id="1"\; label="b"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">b</figure-callout>}1},y_{b1})-f({\mathrm{<figure-callout\; id="1"\; label="x\; f"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">x</figure-callout>}}_f,y_{f1})|}{2}$

Marginal point P (x _p, y _p) corresponding gray scale difference

${\mathrm{<figure-callout\; id="1"\; label="f\; p"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">f</figure-callout>}}_p=\frac{|f(x_{\mathrm{<figure-callout\; id="1"\; label="c"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">c</figure-callout>}1},y_{c1})-f(x_p,y_p)|}{2}+\frac{|f(x_p,y_p)-f(x_{a1},y_{a1})|}{2}$

Neighbor point C ₁(x _C1, y _C1) corresponding gray scale difference

${\mathrm{<figure-callout\; id="1"\; label="f\; c"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">f</figure-callout>}}_c=\frac{|f({\mathrm{<figure-callout\; id="1"\; label="x\; d"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">x</figure-callout>}}_d,y_{d1})-f(x_{\mathrm{<figure-callout\; id="1"\; label="c"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">c</figure-callout>}1},y_{c1})|}{2}+\frac{|f(x_{\mathrm{<figure-callout\; id="1"\; label="c"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">c</figure-callout>}1},y_{c1})-f(x_p,y_p)|}{2}$

Neighbor point D ₁(x _d, y _d) corresponding gray scale difference

${\mathrm{<figure-callout\; id="1"\; label="f\; d"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">f</figure-callout>}}_d=\frac{|f(x_{\mathrm{<figure-callout\; id="1"\; label="e"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">e</figure-callout>}1},y_{e1})-f({\mathrm{<figure-callout\; id="1"\; label="x\; d"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">x</figure-callout>}}_d,y_{d1})|}{2}+\frac{|f({\mathrm{<figure-callout\; id="1"\; label="x\; d"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">x</figure-callout>}}_d,y_{d1})-f(x_{\mathrm{<figure-callout\; id="1"\; label="c"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">c</figure-callout>}1},y_{c1})|}{2}$

(2) if edge pixel point P (x _p, y _p) in the 2nd zone,

Neighbor point A ₂(x _A2, y _A2) corresponding gray scale difference

$f_{a2}=\frac{|f(x_p,y_p)-f(x_{a2},y_{a2})|}{2}+\frac{|f(x_{a2},y_{a2})-f(x_{\mathrm{<figure-callout\; id="2"\; label="b"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">b</figure-callout>}2},y_{b2})|}{2}$

Neighbor point B ₂(x _B2, y _B2) corresponding gray scale difference

${\mathrm{<figure-callout\; id="2"\; label="f\; b"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">f</figure-callout>}}_b=\frac{|f(x_{a2},y_{a2})-f(x_{\mathrm{<figure-callout\; id="2"\; label="b"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">b</figure-callout>}2},y_{b2})|}{2}+\frac{|f(x_{\mathrm{<figure-callout\; id="2"\; label="b"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">b</figure-callout>}2},y_{b2})-f({\mathrm{<figure-callout\; id="2"\; label="x\; f"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">x</figure-callout>}}_f,y_{f2})|}{2}$

Marginal point P (x _p, y _p) corresponding gray scale difference

${\mathrm{<figure-callout\; id="2"\; label="f\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">f</figure-callout>}}_p=\frac{|f(x_{\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">c</figure-callout>}2},y_{c2})-f(x_p,y_p)|}{2}+\frac{|f(x_p,y_p)-f(x_{a2},y_{a2})|}{2}$

Neighbor point C ₂(x _C2, y _C2) corresponding gray scale difference

${\mathrm{<figure-callout\; id="2"\; label="f\; c"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">f</figure-callout>}}_c=\frac{|f({\mathrm{<figure-callout\; id="2"\; label="x\; d"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">x</figure-callout>}}_d,y_{d2})-f(x_{\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">c</figure-callout>}2},y_{c2})|}{2}+\frac{|f(x_{\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">c</figure-callout>}2},y_{c2})-f(x_p,y_p)|}{2}$

Neighbor point D ₂(x _D2, y _D2) corresponding gray scale difference

${\mathrm{<figure-callout\; id="2"\; label="f\; d"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">f</figure-callout>}}_d=\frac{|f(x_{\mathrm{<figure-callout\; id="2"\; label="e"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">e</figure-callout>}2},y_{e2})-f({\mathrm{<figure-callout\; id="2"\; label="x\; d"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">x</figure-callout>}}_d,y_{d2})|}{2}+\frac{|f({\mathrm{<figure-callout\; id="2"\; label="x\; d"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">x</figure-callout>}}_d,y_{d2})-f(x_{\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">c</figure-callout>}2},y_{c2})|}{2}$

The 3.5th step: ask for the range difference of sub-pix point edge point and pixel edge point, acquiring method and being described as follows: according to the square aperture principle, the gray-scale value of a certain pixel is exported and can be expressed as

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

In the formula, (x y) is the light distribution of consecutive image to g, and (i j) is the result of each several part light intensity combined action on the pixel light-sensitive surface to f, and sampled result is to be the discrete matrix of numerical value with the gray-scale value.Because the convolution effect and the optical diffraction effect of optical component, become the form of gradual change through optical imagery in the gray scale of object space drastic change, the edge is characterized by a kind of intensity profile in image, the image border gray-value variation should be a Gaussian distribution, the position of Gaussian curve summit correspondence is the position of true edge point, and for circular image, the gradient direction of pixel edge point is the rectilinear direction that this point is arrived in the center of circle, only need to carry out the position that gaussian curve approximation can be obtained sub-pixel edge point, Gauss's surface fitting of two dimension can be converted into the gaussian curve approximation of one dimension like this to pixel edge point and along the neighbor point of gradient direction.As shown in Figure 6, P (x _p, y _p) point is for the pixel edge point, L _GBe edge pixel point P (x _p, y _p) the gradient direction straight line, if P is in the 1st zone, B ₁, A ₁, C ₁, D ₁Be the neighbor point of P point along gradient direction, if P is in the 2nd zone, B ₂, A ₂, C ₂, D ₂Be the neighbor point of P point along gradient direction, pairing P ' of Gaussian curve summit M should be its true edge point position, and P point and P ' 's range difference is δ.The expression formula of one dimension Gaussian curve is:

$\tilde{y}=\frac{1}{\sqrt{2\mathrm{\&pi;}}{\mathrm{\&sigma;}}_G}\exp \left(\frac{-{(\mathrm{\&delta;}-\mathrm{\&mu;})}^2}{2{{\mathrm{\&sigma;}}_G}^2}\right)---\left(1\right)$

In the formula, μ is the average of Gaussian function, σ _GStandard deviation for Gaussian function.

Calculate for convenient, taken the logarithm in the following formula both sides, and order ${\tilde{y}}^{*}=\ln \tilde{y},$ Following formula can be converted into:

$\stackrel{~*}{y}=m_{11}{\mathrm{\&delta;}}^2+{\mathrm{<figure-callout\; id="12"\; label="m"\; filenames="A20081002279300331.png"\; state="\{\{state\}\}">m</figure-callout>}}_{12}\mathrm{\&delta;}+m_{13}---\left(2\right)$

According to the square aperture sampling thheorem, the pixel grey scale difference is

$\tilde{y}*\left(n\right)={\&Integral;}_{n-0.5}^{n+0.5}(m_{11}{\mathrm{\&delta;}}^2+{\mathrm{<figure-callout\; id="12"\; label="m"\; filenames="A20081002279300331.png"\; state="\{\{state\}\}">m</figure-callout>}}_{12}\mathrm{\&delta;}+m_{13})\mathrm{d\&delta;}---\left(3\right)$

The sequence number that makes marginal point P is 0, and in the 1st zone, the gray scale difference that the P point is corresponding is f as if P _P1, its neighbor point D ₁, C ₁, A ₁, B ₁Sequence number is expressed as-2 ,-1,1 and 2 respectively, and corresponding gray scale difference is f _D1, f _C1, f _A1And f _B1, according to formula (3), have:

${\mathrm{<figure-callout\; id="1"\; label="f\; d"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">f</figure-callout>}}_d={\&Integral;}_{-2.5}^{-1.5}(m_{11}{\mathrm{\&delta;}}^2+{\mathrm{<figure-callout\; id="12"\; label="m"\; filenames="A20081002279300331.png"\; state="\{\{state\}\}">m</figure-callout>}}_{12}\mathrm{\&delta;}+m_{13})\mathrm{d\&delta;}=\frac{49}{12}{m}_{11}-{2\mathrm{<figure-callout\; id="12"\; label="m"\; filenames="A20081002279300331.png"\; state="\{\{state\}\}">m</figure-callout>}}_{12}+m_{13}---\left(4\right)$

${\mathrm{<figure-callout\; id="1"\; label="f\; c"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">f</figure-callout>}}_c={\&Integral;}_{-1.5}^{-0.5}(m_{11}{\mathrm{\&delta;}}^2+{\mathrm{<figure-callout\; id="12"\; label="m"\; filenames="A20081002279300331.png"\; state="\{\{state\}\}">m</figure-callout>}}_{12}\mathrm{\&delta;}+m_{13})\mathrm{d\&delta;}=\frac{13}{12}{m}_{11}-{\mathrm{<figure-callout\; id="12"\; label="m"\; filenames="A20081002279300331.png"\; state="\{\{state\}\}">m</figure-callout>}}_{12}+m_{13}---\left(5\right)$

${\mathrm{<figure-callout\; id="1"\; label="f\; p"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">f</figure-callout>}}_p={\&Integral;}_{-0.5}^{0.5}(m_{11}{\mathrm{\&delta;}}^2+{\mathrm{<figure-callout\; id="12"\; label="m"\; filenames="A20081002279300331.png"\; state="\{\{state\}\}">m</figure-callout>}}_{12}\mathrm{\&delta;}+m_{13})\mathrm{d\&delta;}=\frac{1}{12}{m}_{11}+m_{13}---\left(6\right)$

$f_{a1}={\&Integral;}_{0.5}^{1.5}(m_{11}{\mathrm{\&delta;}}^2+{\mathrm{<figure-callout\; id="12"\; label="m"\; filenames="A20081002279300331.png"\; state="\{\{state\}\}">m</figure-callout>}}_{12}\mathrm{\&delta;}+m_{13})\mathrm{d\&delta;}=\frac{13}{12}{m}_{11}+m_{12}+m_{13}---\left(7\right)$

${\mathrm{<figure-callout\; id="1"\; label="f\; b"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">f</figure-callout>}}_b={\&Integral;}_{1.5}^{2.5}(m_{11}{\mathrm{\&delta;}}^2+{\mathrm{<figure-callout\; id="12"\; label="m"\; filenames="A20081002279300331.png"\; state="\{\{state\}\}">m</figure-callout>}}_{12}\mathrm{\&delta;}+m_{13})\mathrm{d\&delta;}=\frac{49}{12}{m}_{11}+{2\mathrm{<figure-callout\; id="12"\; label="m"\; filenames="A20081002279300331.png"\; state="\{\{state\}\}">m</figure-callout>}}_{12}+m_{13}---\left(8\right)$

According to formula (4)～(8), Simultaneous Equations can be tried to achieve m with least square method ₁₁, m ₁₂, m ₁₃About f _D1, f _C1, f _P1, f _A1And f _B1Expression formula, with its substitution para-curve apex coordinate value δ=-m ₁₂/ 2m ₁₁:

$\mathrm{\&delta;}=-\frac{-0.2f_{d1}-0.1{\mathrm{<figure-callout\; id="1"\; label="f\; c"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">f</figure-callout>}}_c+0.1f_{a1}-0.2{\mathrm{<figure-callout\; id="1"\; label="f\; b"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">f</figure-callout>}}_b}{2({0.1429f}_{d1}-0.0714f_{c1}-0.1429f_{p1}-0.0714f_{a1}+0.1429{\mathrm{<figure-callout\; id="1"\; label="f\; b"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">f</figure-callout>}}_b}$

Attention formula (1) arrives in the process of formula (2),

Passed through the operation of taking from right logarithm, therefore, the pixel grey scale difference also should be taken the logarithm, so P point and P ' 's range difference is:

$\mathrm{\&delta;}=\frac{0.1\mathrm{<figure-callout\; id="1"\; label="ln\; f\; d"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">ln</figure-callout>}{f}_d{1}_{}+0.05\ln f_{c1}-0.05\ln f_{a1}+0.1\ln {\mathrm{<figure-callout\; id="1"\; label="f\; b"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">f</figure-callout>}}_b}{0.1429\ln f_{d1}-0.0714\ln f_{c1}-0.1429\ln f_{p1}-0.714\ln f_{a1}+0.1429\mathrm{<figure-callout\; id="1"\; label="ln\; f\; b"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">ln</figure-callout>}{f}_b{1}_{}}$

If P is in the 2nd zone, the gray scale difference that the P point is corresponding is f _P2, its neighbor point D ₂, C ₂, A ₂, B ₂Sequence number is expressed as-2 ,-1,1 and 2 respectively, and corresponding gray scale difference is f _D2, f _C2, f _A2And f _B2, according to formula (3), have:

${\mathrm{<figure-callout\; id="2"\; label="f\; d"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">f</figure-callout>}}_d={\&Integral;}_{-2.5}^{-1.5}(m_{11}{\mathrm{\&delta;}}^2+{\mathrm{<figure-callout\; id="12"\; label="m"\; filenames="A20081002279300331.png"\; state="\{\{state\}\}">m</figure-callout>}}_{12}\mathrm{\&delta;}+m_{13})\mathrm{d\&delta;}=\frac{49}{12}{m}_{11}-{2\mathrm{<figure-callout\; id="12"\; label="m"\; filenames="A20081002279300331.png"\; state="\{\{state\}\}">m</figure-callout>}}_{12}+m_{13}---\left(9\right)$

${\mathrm{<figure-callout\; id="2"\; label="f\; c"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">f</figure-callout>}}_c={\&Integral;}_{-1.5}^{-0.5}(m_{11}{\mathrm{\&delta;}}^2+{\mathrm{<figure-callout\; id="12"\; label="m"\; filenames="A20081002279300331.png"\; state="\{\{state\}\}">m</figure-callout>}}_{12}\mathrm{\&delta;}+m_{13})\mathrm{d\&delta;}=\frac{13}{12}{m}_{11}-{\mathrm{<figure-callout\; id="12"\; label="m"\; filenames="A20081002279300331.png"\; state="\{\{state\}\}">m</figure-callout>}}_{12}+m_{13}---\left(10\right)$

${\mathrm{<figure-callout\; id="2"\; label="f\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">f</figure-callout>}}_p={\&Integral;}_{-0.5}^{0.5}(m_{11}{\mathrm{\&delta;}}^2+{\mathrm{<figure-callout\; id="12"\; label="m"\; filenames="A20081002279300331.png"\; state="\{\{state\}\}">m</figure-callout>}}_{12}\mathrm{\&delta;}+m_{13})\mathrm{d\&delta;}=\frac{1}{12}{m}_{11}+m_{13}---\left(11\right)$

$f_{a2}={\&Integral;}_{0.5}^{1.5}(m_{11}{\mathrm{\&delta;}}^2+{\mathrm{<figure-callout\; id="12"\; label="m"\; filenames="A20081002279300331.png"\; state="\{\{state\}\}">m</figure-callout>}}_{12}\mathrm{\&delta;}+m_{13})\mathrm{d\&delta;}=\frac{13}{12}{m}_{11}+m_{12}+m_{13}---\left(12\right)$

${\mathrm{<figure-callout\; id="2"\; label="f\; b"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">f</figure-callout>}}_b={\&Integral;}_{1.5}^{2.5}(m_{11}{\mathrm{\&delta;}}^2+{\mathrm{<figure-callout\; id="12"\; label="m"\; filenames="A20081002279300331.png"\; state="\{\{state\}\}">m</figure-callout>}}_{12}\mathrm{\&delta;}+m_{13})\mathrm{d\&delta;}=\frac{49}{12}{m}_{11}+{2\mathrm{<figure-callout\; id="12"\; label="m"\; filenames="A20081002279300331.png"\; state="\{\{state\}\}">m</figure-callout>}}_{12}+m_{13}---\left(13\right)$

According to formula (9)～(13), Simultaneous Equations can be tried to achieve m with least square method ₁₁, m ₁₂, m ₁₃About f _D2, f _C2, f _P2, f _A2And f _B2Expression formula, with its substitution para-curve apex coordinate value δ=-m ₁₂/ 2m ₁₁:

$\mathrm{\&delta;}=-\frac{-0.2f_{d2}-0.1{\mathrm{<figure-callout\; id="2"\; label="f\; c"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">f</figure-callout>}}_c+0.1f_{a2}-0.2{\mathrm{<figure-callout\; id="2"\; label="f\; b"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">f</figure-callout>}}_b}{2(0.1429f_{d2}-0.0714f_{c2}-0.1429f_{p2}-0.0714f_{a2}+0.1429{\mathrm{<figure-callout\; id="2"\; label="f\; b"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">f</figure-callout>}}_b}$

Attention formula (1) arrives in the process of formula (2),

Passed through the operation of taking from right logarithm, therefore, the pixel grey scale difference also should be taken the logarithm, so P point and P ' 's range difference is:

$\mathrm{\&delta;}=\frac{0.1\mathrm{<figure-callout\; id="2"\; label="ln\; f\; d"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">ln</figure-callout>}{f}_d{2}_{}+0.05\ln f_{c2}-0.05\ln f_{a2}+0.1\mathrm{<figure-callout\; id="2"\; label="ln\; f\; b"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">ln</figure-callout>}{f}_b{2}_{}}{0.1429\ln f_{d2}-0.0714\ln f_{c2}-0.1429\ln f_{p2}-0.0174\ln f_{a2}+0.142{9\mathrm{<figure-callout\; id="2"\; label="ln\; f\; b"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">ln</figure-callout>}f}_b}$

In sum, according to gained marginal point in the 3.4th step and along the gray scale difference of the correspondence of the neighbor point of gradient direction, the range difference δ that asks for sub-pixel edge point P ' and pixel edge point P is

If marginal point and along having to be 0 in the gray scale difference of the neighbor point correspondence of gradient direction then makes δ=0, and storage δ=0 number of times Count_zero that occurs.

The 3.6th step: ask for the coordinate of sub-pixel edge point P ', concrete grammar is: according to 3.5 go on foot the slope k of gradient direction straight line of the range difference δ of the sub-pixel edge point P that asks and pixel edge point P and pixel edge point P _Gradient, obtain P point and P ' some coordinate difference δ on x direction and y direction respectively _x, δ _y, formula is as follows:

${\mathrm{\&delta;}}_x=\frac{\mathrm{\&delta;}}{\sqrt{k_{\mathrm{grad}\mathrm{ient}}^2+1}}$

${\mathrm{\&delta;}}_y=\frac{k_{\mathrm{grad}\mathrm{ient}}*\mathrm{\&delta;}}{\sqrt{k_{\mathrm{grad}\mathrm{ient}}^2+1}}$

Wherein, the slope k of the gradient direction straight line of pixel edge point P _GradientCan be by the coordinate (x of this pixel edge point _p, y _p) and coarse positioning central coordinate of circle (x _Oc, y _Oc) obtain, that is:

$k_{\mathrm{grad}\mathrm{ient}}=\frac{y_{\mathrm{oc}}-y_p}{x_{\mathrm{oc}}-x_p}$

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

The 3.7th step: to all edge pixel point P (x _p, y _p), asking for corresponding sub-pixel edge point successively according to～the 3.6 step of the 3.2nd step is P ' (x _p+ δ _x, y _p+ δ _y), wherein, (x _p+ δ _x, y _p+ δ _y) be the coordinate of sub-pixel edge point P '.

The 4th step: the sub-pixel edge point is carried out filtering,, handle with the method for curvature filtering and the method for mean filter respectively at " isolated point " and the noise that sub-pixel edge point occurs.

The 4.1st step: at " isolated point " that sub-pixel edge occurs, handle with the method for curvature filtering, concrete steps are as follows:

(1) ask for the curvature of all sub-pixel edge points, method is as follows: in the sub-pixel edge point of sequential storage, Be the coordinate of sub-pixel edge point P ',

Be any coordinate before P ' the next-door neighbour,

Be any coordinate behind P ' the next-door neighbour, the curvature of sub-pixel edge point P ' is so

$k=\frac{1}{r}=\frac{1}{\sqrt{({(x_0-{\mathrm{<figure-callout\; id="2"\; label="x\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">x</figure-callout>}}_{p_{}^{}})}^2+{(y_0-{\mathrm{<figure-callout\; id="2"\; label="y\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">y</figure-callout>}}_{p_{}^{}})}^2)}}$

Wherein, (x ₀, y ₀) for before P ', P ' the next-door neighbour a bit, any forms the round center of circle behind the next-door neighbour of P ',

$x_0=\frac{a-b+c}{d},$ ${\mathrm{<figure-callout\; id="0"\; label="y"\; filenames="A20081002279300322.png,A20081002279300341.png"\; state="\{\{state\}\}">y</figure-callout>}}_0=\frac{e-f+g}{-d},$

$a=({\mathrm{<figure-callout\; id="1"\; label="x\; p"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">x</figure-callout>}}_{p_{}^{}}+{\mathrm{<figure-callout\; id="2"\; label="x\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">x</figure-callout>}}_{p_{}^{}})({\mathrm{<figure-callout\; id="2"\; label="x\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">x</figure-callout>}}_{p_{}^{}}-{\mathrm{<figure-callout\; id="1"\; label="x\; p"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">x</figure-callout>}}_{p_{}^{}})({\mathrm{<figure-callout\; id="3"\; label="y\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">y</figure-callout>}}_{p_{}^{}}-{\mathrm{<figure-callout\; id="2"\; label="y\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">y</figure-callout>}}_{p_{}^{}})$

$b=({\mathrm{<figure-callout\; id="2"\; label="x\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">x</figure-callout>}}_{p_{}^{}}+{\mathrm{<figure-callout\; id="3"\; label="x\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">x</figure-callout>}}_{p_{}^{}})({\mathrm{<figure-callout\; id="3"\; label="x\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">x</figure-callout>}}_{p_{}^{}}-{\mathrm{<figure-callout\; id="2"\; label="x\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">x</figure-callout>}}_{p_{}^{}})({\mathrm{<figure-callout\; id="2"\; label="y\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">y</figure-callout>}}_{p_{}^{}}-{\mathrm{<figure-callout\; id="1"\; label="y\; p"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">y</figure-callout>}}_{p_{}^{}})$

$c=({\mathrm{<figure-callout\; id="1"\; label="y\; p"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">y</figure-callout>}}_{p_{}^{}}-{\mathrm{<figure-callout\; id="3"\; label="y\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">y</figure-callout>}}_{p_{}^{}})({\mathrm{<figure-callout\; id="2"\; label="y\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">y</figure-callout>}}_{p_{}^{}}-{\mathrm{<figure-callout\; id="1"\; label="y\; p"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">y</figure-callout>}}_{p_{}^{}})({\mathrm{<figure-callout\; id="3"\; label="y\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">y</figure-callout>}}_{p_{}^{}}-{\mathrm{<figure-callout\; id="2"\; label="y\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">y</figure-callout>}}_{p_{}^{}})$

$d=2[({\mathrm{<figure-callout\; id="2"\; label="x\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">x</figure-callout>}}_{p_{}^{}}-{\mathrm{<figure-callout\; id="1"\; label="x\; p"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">x</figure-callout>}}_{p_{}^{}})({\mathrm{<figure-callout\; id="3"\; label="y\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">y</figure-callout>}}_{p_{}^{}}-{\mathrm{<figure-callout\; id="2"\; label="y\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">y</figure-callout>}}_{p_{}^{}})-({\mathrm{<figure-callout\; id="3"\; label="x\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">x</figure-callout>}}_{p_{}^{}}-{\mathrm{<figure-callout\; id="2"\; label="x\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">x</figure-callout>}}_{p_{}^{}})({\mathrm{<figure-callout\; id="2"\; label="y\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">y</figure-callout>}}_{p_{}^{}}-{\mathrm{<figure-callout\; id="1"\; label="y\; p"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">y</figure-callout>}}_{p_{}^{}})]$

$e=({\mathrm{<figure-callout\; id="1"\; label="y\; p"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">y</figure-callout>}}_{p_{}^{}}+{\mathrm{<figure-callout\; id="2"\; label="y\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">y</figure-callout>}}_{p_{}^{}})({\mathrm{<figure-callout\; id="2"\; label="y\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">y</figure-callout>}}_{p_{}^{}}-{\mathrm{<figure-callout\; id="1"\; label="y\; p"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">y</figure-callout>}}_{p_{}^{}})({\mathrm{<figure-callout\; id="3"\; label="x\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">x</figure-callout>}}_{p_{}^{}}-{\mathrm{<figure-callout\; id="2"\; label="x\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">x</figure-callout>}}_{p_{}^{}})$

$f=({\mathrm{<figure-callout\; id="2"\; label="y\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">y</figure-callout>}}_{p_{}^{}}+{\mathrm{<figure-callout\; id="3"\; label="y\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">y</figure-callout>}}_{p_{}^{}})({\mathrm{<figure-callout\; id="3"\; label="y\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">y</figure-callout>}}_{p_{}^{}}-{\mathrm{<figure-callout\; id="2"\; label="y\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">y</figure-callout>}}_{p_{}^{}})({\mathrm{<figure-callout\; id="2"\; label="x\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">x</figure-callout>}}_{p_{}^{}}-{\mathrm{<figure-callout\; id="1"\; label="x\; p"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">x</figure-callout>}}_{p_{}^{}})$

$g=({\mathrm{<figure-callout\; id="1"\; label="x\; p"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">x</figure-callout>}}_{p_{}^{}}-{\mathrm{<figure-callout\; id="3"\; label="x\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">x</figure-callout>}}_{p_{}^{}})({\mathrm{<figure-callout\; id="2"\; label="x\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">x</figure-callout>}}_{p_{}^{}}-{\mathrm{<figure-callout\; id="1"\; label="x\; p"\; filenames="A20081002279300311.png,A20081002279300322.png"\; state="\{\{state\}\}">x</figure-callout>}}_{p_{}^{}})({\mathrm{<figure-callout\; id="3"\; label="x\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">x</figure-callout>}}_{p_{}^{}}-{\mathrm{<figure-callout\; id="2"\; label="x\; p"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">x</figure-callout>}}_{p_{}^{}})$

(2) according to the curvature of asking sub-pixel edge point all sub-pixel edge points are carried out filtering, filter criteria is as follows:

At first the curvature of each sub-pixel edge point is carried out descending sort, get curvature threshold and be this sequence curvature value 3*n element from high to low, n is the number of isolated point, and the reason that causes " isolated point " is because in the 3.5th step, edge pixel point and to have one or more along the neighbor point corresponding gray scale difference of gradient direction be 0, therefore, the number n numerical value of isolated point is the occurrence number Count_zero of δ=0 when sub-pixel edge is located in the 3.5th step; With curvature threshold the curvature of sub-pixel edge point is cut apart then, if the curvature of sub-pixel edge point P greater than curvature threshold, and the curvature of P thinks then that greater than the curvature of next-door neighbour's point before and after it P is an isolated point, gives filtering.

The 4.2nd step: all the sub-pixel edge points after the filtering " isolated point " are carried out filtering with Mean Method, and concrete grammar is: the average of horizontal stroke, ordinate of choosing 2 horizontal stroke, ordinate and the sub-pixel edge point itself that sub-pixel edge point front and back are close to is as new horizontal stroke, the ordinate of this sub-pixel edge point.

The 5th step: filtered sub-pixel edge point is justified with least square fitting, finally obtained the center of circle and the radius of circular target.The algorithm of least square fitting circle is as follows:

If filtered sub-pixel edge point number is N_sub, the general equation of quafric curve is

x ²+2B _exy+C _ey ²+2D _ex+2E _ey+F _e＝0

When utilizing N_sub marginal point to carry out curve fitting, its mean square deviation and be

$e^2=\underset{i=1}{\overset{N\_\mathrm{sub}}{\mathrm{\&Sigma;}}}{({\mathrm{<figure-callout\; id="2"\; label="x\; i"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">x</figure-callout>}}_i^{}+2B_e{x}_i{y}_i+C_e{y}_i^2+2D_e{x}_i+2E_e{y}_i+F_e)}^2$

To following formula about B _e, C _e, D _e, E _e, F _eGet partial derivative respectively, and make that each formula is zero, can obtain a static determinacy system of equations that comprises 5 equations and 5 unknown numbers.Can be with methods such as matrix inversion or the cancellations of Gauss's pivot in a column in the hope of final centre coordinate (x _Cf, y _Cf) be:

$x_{\mathrm{cf}}=\frac{B_e{E}_e-C_e{D}_e}{C_e-{{\mathrm{<figure-callout\; id="2"\; label="B\; e"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">B</figure-callout>}}_e}^2}$

$y_{\mathrm{cf}}=\frac{B_e{D}_e-E_e}{C_e-{{\mathrm{<figure-callout\; id="2"\; label="B\; e"\; filenames="A20081002279300322.png,A20081002279300333.png"\; state="\{\{state\}\}">B</figure-callout>}}_e}^2}.$

Claims (2)

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

The 1st step: the center of circle and radius to circular target in the image carry out coarse positioning, obtain the coarse positioning center of circle (x of circular target _Oc, y _Oc) and coarse positioning radius of circle R _c, wherein, (x _Oc, y _Oc) be the coordinate in the center of circle under the image coordinate system,

The 2nd step: the circular target in the image is carried out the pixel edge location, and concrete steps are as follows:

The 2.1st step: according to the coarse positioning center of circle and coarse positioning radius of circle, from image, extract a square area, be called " area-of-interest ", the extracting method of this square area is: with the coarse positioning center of circle of circular target in the image central point as square area, coarse positioning diameter with circular target adds the length of side of 2～6 pixels as square area

The 2.2nd step: with the canny operator area-of-interest that is extracted in the image is carried out pixel edge and detect, obtain the pixel edge of circular target, again the pixel edge point of the circular target of canny operator extraction being carried out the border follows the tracks of, obtain the coordinate of each pixel edge point, and coordinate by clockwise each pixel edge point of sequential storage

The 3rd step: the circular target in the image is carried out the sub-pixel edge location, and concrete steps are as follows:

The 3.1st step: the edge pixel point of circular target is carried out area dividing, and division methods is as follows: claim the coarse positioning center of circle with circular target to be initial point, to be positive dirction, to be that the straight line of unit length is the x axle with the pixel with level to the right; Title with the coarse positioning center of circle of circular target be initial point, to be positive dirction vertically upward, to be that the straight line of unit length is the y axle with the pixel; With the coarse positioning center of circle is the point of rotation, the angle that is rotated counterclockwise by x axle positive dirction is θ, edge pixel point according to big young pathbreaker's circular target of θ is divided into two parts, a part is [0 ° of θ ∈, 45 °] [135 ° of ∪, 225 °] [315 ° of ∪, 360 °] pairing edge pixel point, the edge pixel point that is called the 1st zone, another part are (225 ° of θ ∈ (45 °, 135 °) ∪, 315 °) corresponding edge pixel point, be called the edge pixel point in the 2nd zone, the line of the coarse positioning center of circle and edge pixel point is the gradient direction of this edge pixel point

The 3.2nd step: ask for the neighbor point of edge pixel point, (x along gradient direction _p, y _p) for being the coordinate figure of the marginal point P of unit with the pixel, concrete grammar is as follows: if edge pixel point P (x _p, y _p) in the 1st zone, get the gradient direction straight line and the straight line x=x of this pixel _p+ 1 intersection point A ₁(x _A1, y _A1) be the right first neighbor point of this edge pixel point along gradient direction; Get the gradient direction straight line and the straight line x=x of this pixel _p+ 2 intersection points B ₁(x _B1, y _B1) be the right second neighbor point of this edge pixel point along gradient direction; Get the gradient direction straight line and the straight line x=x of this pixel _p+ 3 intersection point F (x _f, y _f) be the right three neighbor point of this edge pixel point along gradient direction; Get the gradient direction straight line and the straight line x=x of this pixel _p-1 intersection point C ₁(x _C1, y _C1) be the left side first neighbor point of this edge pixel point along gradient direction; Get the gradient direction straight line and the straight line x=x of this pixel _p-2 intersection point D ₁(x _D1, y _D1) be the left side second neighbor point of this edge pixel point along gradient direction, get the gradient direction straight line and the straight line x=x of this pixel _p-3 intersection point E ₁(x _E1, y _E1) be the left side three neighbor point of this edge pixel point along gradient direction; If edge pixel point P (x _p, y _p) in the 2nd zone, get the gradient direction straight line and the straight line y=y of this pixel _p+ 1 intersection point A ₂(x _A2, y _A2) be the top first next-door neighbour point of this edge pixel point along gradient direction; Get the gradient direction straight line and the straight line y=y of this pixel _p+ 2 intersection points B ₂(x _B2, y _B2) be the top second next-door neighbour point of this edge pixel point along gradient direction; Get the gradient direction straight line and the straight line y=y of this pixel _p+ 3 intersection point F ₂(x _F2, y _F2) be the top three next-door neighbour point of this edge pixel point along gradient direction; Get the gradient direction straight line and the straight line y=y of this pixel _p-1 intersection point C ₂(x _C2, y _C2) be the bottom first next-door neighbour point of this edge pixel point along gradient direction; Get the gradient direction straight line and the straight line y=y of this pixel _p-2 intersection point D ₂(x _D2, y _D2) be the top second next-door neighbour point of this edge pixel point along gradient direction; Get the gradient direction straight line and the straight line y=y of this pixel _p-3 intersection point E ₂(x _E2, y _E2) be the top three next-door neighbour point of this edge pixel point along gradient direction,

The 3.3rd step: adopting the method for linear gray-level interpolation to ask for the pixel is the edge pixel point P (x of unit _p, y _p) along the gray-scale value of the neighbor point of gradient direction, wherein, (x _p, y _p) for being the coordinate figure of the edge pixel point P of unit with the pixel, use f (x, y) denotation coordination be (use symbol [] is represented the round numbers part for x, gray values of pixel points y), and the described linear gray-level interpolation method that is used to obtain the neighbor point gray-scale value is as follows:

(1) if edge pixel point P (x _p, y _p) in the 1st zone,

Neighbor point A ₁(x _A1, y _A1) gray-scale value

f(x _a1，y _a1)＝(1-λ)*f(x _a1，[y _a1])+λ*f(x _a1，[y _a1]+1)

λ＝y _a1-[y _a1]，

Neighbor point B ₁(x _B1, y _B1) gray-scale value

f(x _b1，y _b1)＝(1-λ)*f(x _b1，[y _b1])+λ*f(x _b1，[y _b1]+1)

λ＝y _b1-[y _b1]

Neighbor point F ₁(x _F1, y _F1) gray-scale value

f(x _f1，y _f1)＝(1-λ)*f(x _f1，[y _f1])+λ*f(x _f1，[y _f1]+1)

λ＝y _f1-[y _f1]

Neighbor point C ₁(x _C1, y _C1) gray-scale value

f(x _c1，y _c1)＝(1-λ)*f(x _c1，[y _c1])+λ*f(x _c1，[y _c1]+1)

λ＝y _c1-[y _c1]

Neighbor point D ₁(x _d, y _d) gray-scale value

f(x _d1，y _d1)＝(1-λ)*f(x _d1，[y _d1])+λ*f(x _d1，[y _d1]+1)

λ＝y _d1-[y _d1]

Neighbor point E ₁(x _E1, y _E1) gray-scale value

f(x _e1，y _e1)＝(1-λ)*f(x _e1，[y _e1])+λ*f(x _e1，[y _e1]+1)

λ＝y _e1-[y _e1]

(2) if edge pixel point P (x _p, y _p) in the 2nd zone,

Neighbor point A ₂(x _A2, y _A2) gray-scale value

f(x _a2，y _a2)＝(1-λ)*f([x _a2]，y _a2)+λ*f([x _a2]+1，y _a2)

λ＝x _a2-[x _a2]

Neighbor point B ₂(x _B2, y _B2) gray-scale value

f(x _b2，y _b2)＝(1-λ)*f([x _b2]，y _b2)+λ*f([x _b2]+1，y _b2)

λ＝x _b2-[x _b2]

Neighbor point F ₂(x _F2, y _F2) gray-scale value

f(x _f2，y _f2)＝(1-λ)*f([x _f2]，y _f2)+λ*f([x _f2]+1，y _f2)

λ＝x _f2-[x _f2]

Neighbor point C ₂(x _C2, y _C2) gray-scale value

f(x _c2，y _c2)＝(1-λ)*f([x _c2]，y _c2)+λ*f([x _c2]+1，y _c2)

λ＝x _c2-[x _c2]

Neighbor point D ₂(x _D2, y _D2) gray-scale value

f(x _d2，y _d2)＝(1-λ)*f([x _d2]，y _d2)+λ*f([x _d2]+1，y _d2)

λ＝x _d2-[x _d2]

Neighbor point E ₂(x _E2, y _E2) gray-scale value

f(x _e2，y _e2)＝(1-λ)*f([x _e2]，y _e2)+λ*f([x _e2]+1，y _e2)

λ＝x _e2-[x _e2]

The 3.4th step: according to the 3.3rd step the marginal point of asking along the gray-scale value and the marginal point P (x of the neighbor point of gradient direction _p, y _p) gray-scale value itself, the mean value of choosing the forward difference of gray-scale value and backward difference is as marginal point and along the gray scale difference of the correspondence of the neighbor point of gradient direction, wherein, and (x _p, y _p) for being the coordinate figure of the edge pixel point P of unit with the pixel, use f (x, y) denotation coordination be (the described method of asking for gray scale difference is as follows for x, gray values of pixel points y):

(1) if edge pixel point P (x _p, y _p) in the 1st zone,

Neighbor point A ₁(x _A1, y _A1) corresponding gray scale difference

$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}$

Neighbor point B ₁(x _B1, y _B1) corresponding gray scale difference

$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}$

Marginal point P (x _p, y _p) corresponding gray scale difference

$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}$

Neighbor point C ₁(x _C1, y _C1) corresponding gray scale difference

$f_{c1}=\frac{|f(x_{d1},y_{d1})-f(x_{c1},y_{c1})|}{2}+\frac{|f(c_{c1},y_{c1})-f(x_p,y_p)|}{2}$

Neighbor point D ₁(x _d, y _d) corresponding gray scale difference

$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) if edge pixel point P (x _p, y _p) in the 2nd zone,

Neighbor point A ₂(x _A2, y _A2) corresponding gray scale difference

$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}$

Neighbor point B ₂(x _B2, y _B2) corresponding gray scale difference

$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}$

Marginal point P (x _p, y _p) corresponding gray scale difference

$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}$

Neighbor point C ₂(x _C2, y _C2) corresponding gray scale difference

$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}$

Neighbor point D ₂(x _D2, y _D2) corresponding gray scale difference

$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}$

The 3.5th step: according to gained marginal point in the 3.4th step and along the gray scale difference of the correspondence of the neighbor point of gradient direction, the range difference δ that asks for sub-pixel edge point P ' and pixel edge point P is

If marginal point and along having to be 0 in the gray scale difference of the neighbor point correspondence of gradient direction then makes δ=0, and storage δ=0 number of times Count_zero that occurs,

The 3.6th step: ask for the coordinate of sub-pixel edge point P ', concrete grammar is: according to 3.5 go on foot the slope k of gradient direction straight line of the range difference δ of the sub-pixel edge point P ' that asks and pixel edge point P and pixel edge point P _Gradient, obtain P point and P ' some coordinate difference δ on x direction and y direction respectively _x, δ _y, formula is as follows:

${\mathrm{\&delta;}}_x=\frac{\mathrm{\&delta;}}{\sqrt{k_{\mathrm{grad}\mathrm{ient}}^2+1}}$

${\mathrm{\&delta;}}_y=\frac{k_{\mathrm{grad}\mathrm{ient}}*\mathrm{\&delta;}}{\sqrt{k_{\mathrm{grad}\mathrm{ient}}^2+1}}$

Wherein, the slope k of the gradient direction straight line of pixel edge point P _GradientCan be by the coordinate (x of this pixel edge point _p, y _p) and coarse positioning central coordinate of circle (x _Oc, y _Oc) obtain, that is:

$k_{\mathrm{grad}\mathrm{ient}}=\frac{y_{\mathrm{oc}}-y_p}{x_{\mathrm{oc}}-x_p}$

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

The 3.7th step: to all edge pixel point P (x _p, y _p), asking for corresponding sub-pixel edge point successively according to～the 3.6 step of the 3.2nd step is P ' (x _p+ δ _x, y _p+ δ _y), wherein, (x _p+ δ _x, y _p+ δ _y) be the coordinate of sub-pixel edge point P ',

The 4th step: the sub-pixel edge point is carried out filtering,, handle with the method for curvature filtering and the method for mean filter respectively at " isolated point " and the noise that sub-pixel edge point occurs,

The 4.1st step: at " isolated point " that sub-pixel edge occurs, handle with the method for curvature filtering, concrete steps are as follows:

(1) ask for the curvature of all sub-pixel edge points, method is as follows: in the sub-pixel edge point of sequential storage,

Be the coordinate of sub-pixel edge point P ', Be any coordinate before P ' the next-door neighbour,

Be any coordinate behind P ' the next-door neighbour, the curvature of sub-pixel edge point P ' is so

$k=\frac{1}{r}=\frac{1}{\sqrt{({(x_0-x_{p_2^{\&prime;}})}^2+{(y_0-y_{p_2^{\&prime;}})}^2)}}$

Wherein, (x ₀, y ₀) for before P ', P ' the next-door neighbour a bit, any forms the round center of circle behind the next-door neighbour of 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) according to the curvature of sub-pixel edge point all sub-pixel edge points are carried out filtering, filter criteria is as follows:

At first the curvature of each sub-pixel edge point is carried out descending sort, get curvature threshold and be this sequence curvature value 3*n element from high to low, wherein, n is the number of isolated point, and its numerical value is the occurrence number Count_zero of δ=0 when sub-pixel edge is located in the 3.5th step; With curvature threshold the curvature of sub-pixel edge point is cut apart then, if the curvature of sub-pixel edge point P greater than curvature threshold, and the curvature of P thinks then that greater than the curvature of next-door neighbour's point before and after it P is an isolated point, gives filtering,

The 4.2nd step: all the sub-pixel edge points after the filtering " isolated point " are carried out filtering with Mean Method, concrete grammar is: the average of horizontal stroke, ordinate of choosing next-door neighbour's 2 horizontal stroke, ordinate and sub-pixel edge point itself before and after the sub-pixel edge point is as new horizontal stroke, the ordinate of this sub-pixel edge point

The 5th step: filtered sub-pixel edge point is justified with least square fitting, finally obtained the center of circle and the radius of circular target.

2, the circle center locating method of circular target according to claim 1 is characterized in that:

The 1st step, the described method that the center of circle and the radius of circular target in the image are carried out coarse positioning was:

The 1.1st step: image removed make an uproar, Threshold Segmentation, obtain each pixel gray-scale value in the image and be 255 or 0 bianry image,

The 1.2nd step: the circular target in the bianry image is carried out Boundary Extraction and border tracking,

The 1.3rd step: round to the frontier point of following the tracks of with the match of contour centroid method, obtain the coarse positioning center of circle and the coarse positioning radius of circle of this circular target, and store coarse positioning central coordinate of circle (x _Oc, y _Oc) and coarse positioning radius of circle R _c

Images