- Email: [email protected]
Effective algorithms for circle fitting
Computer Physics Communications 33(1984) 329—333 North-Holland, Amsterdam
329
EFFECTIVE ALGORITHMS FOR CIRCLE FITTING N.!. CHERNOV, G.A. OSOSKOV Joint institute for Nuclear research, Head Post Office, P.O. Box 79, 101000 Moscow, USSR
Received 20 February 1984
After an investigation of known methods for circle fitting two new effective algorithms are proposed which meet all the requirements of mass data handling.
1. In many data handling problems there is a need to approximate data points by a circle. Some cases (pictures of lunar craters [11, eye cornea surface [2], etc.) are related to circle-like images, but more frequently a circular arc is used as an approximating curve that possesses some advantages when compared with other curves even allowing for their simplicity in fitting. One can consider as an essential example the problem of automatic recognition of elementary particle trajectories on images of multiple track events produced by some physical device that detects particle collisions. This problem is quite typical in high energy physics. The data are fed to a computer in the form of a set of point coordinates. The problem is then to identify points lying on the same track and to find all tracks belonging to the given event. To do this one has to follow these tracks through all the regions of low contrast or crossings with other tracks, noise spots and other ambiguities. Since the projection of a charged particle moving in the homogeneous magnetic field is very close to some circle, it is natural to use the circle equation in track-following proximate thethe points already found.process to apThis approximation is important for two reasons: (i) to check whether all the points belong to the same track, (ii) to extrapolate the track in order to predict a region in which new points can be found to be added to the track-bank.
The circle is better for such a prediction than the parabola, for instance, because the circle retains its curvature. It is necessary to point out the extremely rigid requirements on the efficiency of the fitting algorithm in tracking problems. This follows from the very high rate at which the fitting procedure is called (~iO~times for each frame, with 105_106 frames for one experiment). A satisfactory solution of this problem, however, was not yet available, mainly because of its nonlinearity. In the present paper, after presenting the known algorithms, two new algorithms are proposed for full circle-fitting and also for the circular arc. The comparative characteristics (accuracy, numbers of operations, computer time) of the proposed algorithms show their high efficiency and suitability not only for off-line computations but also for on-line applications on minicomputers. 2. Let us formulate the problem. For n points given on the plane by their coordinates (x1, y,), i =
~IT~ it is necessary 2 + (y
b)2
to draw the circle R2
(1) (x a) minimizing the mean square distance from the given points to this circle, i.e. to minimize the functional —
=
—
L(a, b, R)
OO1O-4655/84/$03.OO © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)
~ p~,
=
I
(2)
330
N.J. Chernov, GA. Ososkov
/
where p,
=
f(x,
Effective algorithms for circlefitting
serious disadvantage of PM is its in-principal map—
a)2 +(y,
—
b)2
—
R.
Due to the nonlinear dependence of p, on the circle parameters the application of the least square method (LSM) leads to cumbersome computations. Usually, the LSM is carried over as an iterative procedure, applying linearisation at each step and demanding on elaborate choice of the first approximation (a 0, b0, R0) (see, for example, refs. [3—5]). The corresponding subroutine is very time-consuming and therefore in practical applications LSM is not used or is applied at the final stages of data handling [2,5]. A number of simplified fast algorithms were proposed [2,5,7], but they are worse than LSM from the point of view of parameter accuracy. We shall describe two of them since this turns out to be useful later on. 3. The first method is valid only for small circular arcs (<300) and is based on the least square approximation to the set (x1, y,), i 1,n by a parabola =
y=px2+qx+r,
(4)
plicability in the along frequently when points are placed a full occurring circle or oncase some big arc of it. The second method is based on introducing the new variable Z in the circle equation to transform it to x2 the+y2 linear regression equation(1) =
ax + /3y + y, where a 2a, $
z
=
=
2b, y
=
=
=
—
a2
—
b2 are treated
as the new unknown parameters. Evaluating them by LSM one can find the initial parameters a, b, R [8].Let us derive the corresponding formulae which will be needed later. To simplify the notation les Ut denote by Gauss brackets expressions like
X P~q
=
[x Py ~]
and suppose that the origin of the coordinate system has been transferred to the centre of grayity of the point set (x1, y,); hence [x] [y] 0. The normal equations for the parameters a, ~8,y are 2]+$[xy]= [x(x2+y2)1, a[x a[xy]+$[y21 [y(x2+y2)j, (5) =
=
after rotating the coordinate system in order to have the axis OX parallel to the straight line going “along” the arc. This line can be drawn either via two end-points (x 1, y1), (xe, y~)or by the LSM, using the full point set [2,5]. The approximating circle in question is the circle tangential to the parabola (4) at its top point and having the radius R l/2p. For the sake of brevity we shall call this method “parabolic” and denote it by PM. In ref. [5] the circle parameters obtained by PM are used as the first approximation for LSM iterations. PM is clearly faster than LSM but the accuracy of circle parameter estimates is definitely worse (see table 1 below). Even in the case when all points lie on the same circle PM cannot give an exact result. Therefore some improvements in the PM parameters were suggested [7], which inevitably demanded additional computer resources. The
R2
=
[x2
+
y2
=
1.
We shall refer to the second method of the circle fitting as the linear regression method (LRM). LRM is fast enough and in addition has the advantage of giving the exact result in the case when all points (x,, y,) lie on the same circle. Nevertheless there is a very essential disadvantage of LRM which should be explained in detail. It is clear that LRM consists in minimizing the functional 2+(y M(a,b,R)=~((xe_a) I = I
2_R2)(6) 1_b)
that is, according to (3) n
M(a, b, R)
=
~
((R
+ p,)2
—
R2)
I =
=
4R2
~ ,=1
p~(1+ p,/2R)2.
N.J. Chernov, GA. Ososkov
/
The usual assumption in any practical application is that p1 << R, otherwise it will be impossible
331
Effective algorithms for circle fitting
Equating to zero the derivatives of this expression
to get a reasonable result. So using (2) we can obtain 2 ~n p~= 4R2L(a, b, R). (7)
by a, b, R we obtain the normal equations, which form a second-order nonlinear system: Fa + Hb ay = Ha + Gb by = Q, (9)
M(a, b, R)
2Pa + 2Qb + .~2= T,
4R
Thus, for LRM it is more “profitable” to minimize not the mean squares L(a, b, R), but R2, whose contribution to the functional (7) is more valuable. This leads to the above-mentioned disadvantage of LRM, i.e. its low stability a superfluous sensitivity to any small errors in measurements. However, it has to be pointed out that the accuracy of the LRM estimations increases rapidly with increasing arc length, and if the points (x,, ~) lie along the full circle the LRM estimations get very close to the true LSM. The method can be generalized to the case of fitting any algebraic curve PN(x, y) = 0, where ~N is a polynomial of the Nth degree. For example, for N = 2 we have the curve equation —
—
—
where as above y = R2 a2 Gauss brackets we denote —
F=
-~-
n
H=~[xy], ~
b2. Using again the
~ [x~ + 3y2]
n
P=~~~[x(x2+y2)J,
[y(x2
!
G=
+ y2]
[3x2
—
+y2)],
~-
!
[(x2
In the system (9) one can exclude any pair of unknown variables to get an equation of the 4th degree. The most suitable choice is to exclude a, b and thereby obtain the equation (10)
y4+Ay3+By2+C-y+D=0 =
a
+y2)2].
2+a 1x
2xy + a3x + a4y + a5.
2=z,x2=u, Theintroductionofnewvariablesy zy = t allows us to obtain the linear regression equation. Then the parameters a 1—a5 are evaluated by solving the corresponding linear normal equation system. Some versions of this approach, leading to a very fast algorithm, will be applied below in no. 5. 4. After the above discussion we consider our own algorithms for circle fitting. The first one is universal, i.e. it is well-suited for fitting any size of circular arc including the full circle. The second is faster, but works only for small size arcs. It will be described in the next section. The idea of the first algorithm is close to the linear regression approach, but it improves the estimation accuracy essentially. Using (7) we shall minimize the functional 2 K(a, b, R)
=
M(a, b, R)R =
+
~
R
—
a
1 Yo
b
)
(8)
1).
=
4+A
3+B 0x
2+C 0x
B0
=
C 0
=
(11) y/y~with the coeffi-
0x+D0=0
related to the variable x
cients A 0 = A/y0,
—
2
2
([x + [~2 It is easy to obtain ~y1~ as a by-product when calculating the coefficiei~itsof (10). Since y~is the first good approximation we can minimize roundoff errors dividing all coefficients of (10) by y~. The new equation is
2-~x, 2—y,
2+b2—R2 R
In spite of the apparent complexity of this method in comparison with PM or LRM we succeeded to find a very fast numerical solution of (10), which now will be briefly described. The eq. (10) has 4 roots; in order to obtain the first approximation to the root we look for, we can use the solution of the normal equations (5)
x a
____
=
H2, C= coefficients T(F+ G)—2(P2 + Q2), T(H2 — with the A = —F— G, BD= = FG — T— FG) + 2(P2G + Q2F) — 4PQH.
C/y~, D0 = D/y.
=
N.J. Chernov, GA. Ososkov
332
/
It is easy to obtain A0 = —4. The evaluations of other coefficients are rather cumbersome, so we shall result limit ourseles to giving here only the final 3
—
Effective algorithms for circle fitting
be rewritten as 2 K(a, b, R)
=
4L(a, b, R).(13)
M(a, b, R)b~
Thus, in order to find the arc parameters in
n~ B
0 ~ 3, 4 — 2n —4n~D0~<2.5n—3.
~
C0
~ 4n,
question,
it
the functional
is sufficient to minimize
~ (x1 b
k(a, b, R)= ~
—
a 2~x,
—
)
For solving (11) any iterative procedure with 1theis value suitable. We root use the of the by 2—5 iterations with the maximal accuracy 10_i3 admissible on the CDC 6500. Calculating y = y 0x we obtain a and b from (9) and 2 + y. (12)
We have been really happy to find that the normal equations for (14) were linear (!) relative to the unknown parameters a, bandy = R2 a2 h2.
R
2aP+2bQ+y([x2]+[y2])=
i~12 + b2
the starting valuegiving x0 = Newton method
=
~/a2 + b
Thus, our first algorithm consists of 5 stages: (i) Transfer the origin of the coordinate system to the centre of gravity of the set (x,,y,); (ii) Calculate the Gauss-bracket values for [x2], [y2], [xy], [x(x2 +y2)], [y(x2 +y2)], [(x2 +
—
a
—
b
R2 2
(14)
—
2 a [x2] + 2 b [ xy] ny = [x2] + [y2]
=
—
T,
(15)
P,
Here we use again the previous notations P,
Q,
T.
(iii) Calculate coefficients for eq. (11) (iv) Solve (11) by the Newton method with the starting value x 0 = 1; (v) Calculate the parameters y, a, b, R. The Fortran-IV program implementing these 5 stages consists of 90 operations. The only external call to the SQRT routine (for calculating R by (12)) is not obligatory because the decision in 2 insteadrules of R, the majority of applications use R so that the SQRT call can be omitted. The iterational method just described will be referred to as the ILRM Iterational Linear Regression Method.
Obtaining a, b, R from (15) is a trivial task, so that we omit the corresponding formulae.
5. The second algorithm for circle fitting is valued only in the case of a small circular arc
point distances from the circle during simulation wasThe taken as 1. comparative—estimate results of parameters obtained by PM, LRM, ILRM and MLRM are
—
(<400).
Let us suppose from the outset that the (x,, y.) points have their centre of gravity located at the origin of the coordinate system and that the OX axis coincides with the straight line drawn by the LSM our points (when these are more or lessthrough equidistant) via theortwo end-points. In such conditions the arc is almost horizontal and clearly Ibi R. Therefore the relation (7) can
We shall call this method the Modified Linear Regression method MLRM. 6. An investigation was made on the CDC-6500 computer in order to compare the characteristics of all the methods As a was model for the experimental data, described. a test of points simulated according to the normal distribution along the circular arc with a = b = 0, R = 1000. The length of the arc was about 400 (from the point (0, 1000) to (400, 100V~),i.e. about 25°).The number of points varied from n = 4 to n = 50. The variance of
summarized in table 1 in the form of their ratios to LSM estimates which were assumed to be the most exact. The variance estimate was calculated as
a2
=
~
—
a)2 +(y,
—
b)2
—
R). (16)
I =
A good approximation for (16) can be derived
N.!. Chernov, G.A. Ososkov
/
Table 1 Comparative—estimate results of parameters obtained by PM. LRM, ILRM and MLRM
Table 3 Number of operations required for execution of each algorithm +,—
a, b, R estimates a2 estimate _____________________________________________ PM 1% 2% LRM 2% 3% IRLM 0.01 + 0.02% 0.001% MLRM 0.05% 0.005%
from (3), (6) and (7) ((x,
“
-2
=
+(y a)24~2(n 1
—
2
— —
2 —
R2)
b) 3)
(17)
I
where a, 1 and ~kare the parameter estimates from MLRM (see also ref. [6]). The computer-time costs for each method were found by the CDC 6500 time-counting routine SECOND. Results in 10-i s are given in table 2. The resulting comparison shows that ILRM is obviously not much slower than PM, so that the effort in improving on PM (see no. 3) is almost useless. One can also note that in the conditions of our test example, the MLRM accuracy decreases sharply when the arc angle increases; MLRM is practically as fast as the fastest LRM and as exact as ILRM. The numbers of operations required for executions of each algorithm are shown in table 3. These numbers were found according to the common rule: addition and subtraction are treated as one operation, multiplication and division as 3, square root as 25 operations. The symbol i in the ILRM row is the iteration number of the Newton-method-solving equation (11). The mean value of i in our test example was 2.5. However, in order to evaluate it in the general case and also to check the stability of ILRM with respect to different point configurations a numeriTable 2 Computer time in
uo—~
n
4
10
50
LSM PM LRM ILRM MLRM
8 1.1 0.7 1.2 0.8
17 2.0 1.5 2.0 1.6
75 8.1 5.9 6.6 6.0
333
Effective algorithms for circle fitting
LSM PM LRM ILRM MLRM
X,/
97n+169 149n+245 13n+29 9n+49 un +8 6n +20 un + 16i + 21 6n + 7i + 42 lln+14 6n+39
[
total
7n+1 I 1 1 I
719n+537 40n+201 29n +93 29n + 37i + 172 29n+156
cal experiment was made. 10000 x 10 simulated points were uniformly distributed in the square (0< x <1,0
.
[1] W.A. Astarev et al., Thesisy I seminara p0 automatlzatzyi nauchnykh issledovanii v jadernoi fizike i smezhnykh oblastjakh (Donish. Dushanbe. 1980) p. 267 (in Russian). [2] G.A. Ososkov, JINR-83-187 Dubna (1983) p. 40 (in Russian). [31B.L. van der Warden, Matematische Statistik (SpringerVerlag, Berlin, 1957). [4] SN. Sokolov and IN. Sum, JINR-D-810. Dubna (1960) (in Russian). [5] I.M. Lipkin et al., Preprint ITEP no. 47. Moscow (1975) (in Russian). [6] D. Hall, UCRL-118544, ~erkeIey,California (1968). [7] N.N. Govorun et al., JINR-1101, Dubna (1962) (in Russian). [8] N.N. Govorun et al., JINR-2036, Dubna (1965) (in Russian).
329
EFFECTIVE ALGORITHMS FOR CIRCLE FITTING N.!. CHERNOV, G.A. OSOSKOV Joint institute for Nuclear research, Head Post Office, P.O. Box 79, 101000 Moscow, USSR
Received 20 February 1984
After an investigation of known methods for circle fitting two new effective algorithms are proposed which meet all the requirements of mass data handling.
1. In many data handling problems there is a need to approximate data points by a circle. Some cases (pictures of lunar craters [11, eye cornea surface [2], etc.) are related to circle-like images, but more frequently a circular arc is used as an approximating curve that possesses some advantages when compared with other curves even allowing for their simplicity in fitting. One can consider as an essential example the problem of automatic recognition of elementary particle trajectories on images of multiple track events produced by some physical device that detects particle collisions. This problem is quite typical in high energy physics. The data are fed to a computer in the form of a set of point coordinates. The problem is then to identify points lying on the same track and to find all tracks belonging to the given event. To do this one has to follow these tracks through all the regions of low contrast or crossings with other tracks, noise spots and other ambiguities. Since the projection of a charged particle moving in the homogeneous magnetic field is very close to some circle, it is natural to use the circle equation in track-following proximate thethe points already found.process to apThis approximation is important for two reasons: (i) to check whether all the points belong to the same track, (ii) to extrapolate the track in order to predict a region in which new points can be found to be added to the track-bank.
The circle is better for such a prediction than the parabola, for instance, because the circle retains its curvature. It is necessary to point out the extremely rigid requirements on the efficiency of the fitting algorithm in tracking problems. This follows from the very high rate at which the fitting procedure is called (~iO~times for each frame, with 105_106 frames for one experiment). A satisfactory solution of this problem, however, was not yet available, mainly because of its nonlinearity. In the present paper, after presenting the known algorithms, two new algorithms are proposed for full circle-fitting and also for the circular arc. The comparative characteristics (accuracy, numbers of operations, computer time) of the proposed algorithms show their high efficiency and suitability not only for off-line computations but also for on-line applications on minicomputers. 2. Let us formulate the problem. For n points given on the plane by their coordinates (x1, y,), i =
~IT~ it is necessary 2 + (y
b)2
to draw the circle R2
(1) (x a) minimizing the mean square distance from the given points to this circle, i.e. to minimize the functional —
=
—
L(a, b, R)
OO1O-4655/84/$03.OO © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)
~ p~,
=
I
(2)
330
N.J. Chernov, GA. Ososkov
/
where p,
=
f(x,
Effective algorithms for circlefitting
serious disadvantage of PM is its in-principal map—
a)2 +(y,
—
b)2
—
R.
Due to the nonlinear dependence of p, on the circle parameters the application of the least square method (LSM) leads to cumbersome computations. Usually, the LSM is carried over as an iterative procedure, applying linearisation at each step and demanding on elaborate choice of the first approximation (a 0, b0, R0) (see, for example, refs. [3—5]). The corresponding subroutine is very time-consuming and therefore in practical applications LSM is not used or is applied at the final stages of data handling [2,5]. A number of simplified fast algorithms were proposed [2,5,7], but they are worse than LSM from the point of view of parameter accuracy. We shall describe two of them since this turns out to be useful later on. 3. The first method is valid only for small circular arcs (<300) and is based on the least square approximation to the set (x1, y,), i 1,n by a parabola =
y=px2+qx+r,
(4)
plicability in the along frequently when points are placed a full occurring circle or oncase some big arc of it. The second method is based on introducing the new variable Z in the circle equation to transform it to x2 the+y2 linear regression equation(1) =
ax + /3y + y, where a 2a, $
z
=
=
2b, y
=
=
=
—
a2
—
b2 are treated
as the new unknown parameters. Evaluating them by LSM one can find the initial parameters a, b, R [8].Let us derive the corresponding formulae which will be needed later. To simplify the notation les Ut denote by Gauss brackets expressions like
X P~q
=
[x Py ~]
and suppose that the origin of the coordinate system has been transferred to the centre of grayity of the point set (x1, y,); hence [x] [y] 0. The normal equations for the parameters a, ~8,y are 2]+$[xy]= [x(x2+y2)1, a[x a[xy]+$[y21 [y(x2+y2)j, (5) =
=
after rotating the coordinate system in order to have the axis OX parallel to the straight line going “along” the arc. This line can be drawn either via two end-points (x 1, y1), (xe, y~)or by the LSM, using the full point set [2,5]. The approximating circle in question is the circle tangential to the parabola (4) at its top point and having the radius R l/2p. For the sake of brevity we shall call this method “parabolic” and denote it by PM. In ref. [5] the circle parameters obtained by PM are used as the first approximation for LSM iterations. PM is clearly faster than LSM but the accuracy of circle parameter estimates is definitely worse (see table 1 below). Even in the case when all points lie on the same circle PM cannot give an exact result. Therefore some improvements in the PM parameters were suggested [7], which inevitably demanded additional computer resources. The
R2
=
[x2
+
y2
=
1.
We shall refer to the second method of the circle fitting as the linear regression method (LRM). LRM is fast enough and in addition has the advantage of giving the exact result in the case when all points (x,, y,) lie on the same circle. Nevertheless there is a very essential disadvantage of LRM which should be explained in detail. It is clear that LRM consists in minimizing the functional 2+(y M(a,b,R)=~((xe_a) I = I
2_R2)(6) 1_b)
that is, according to (3) n
M(a, b, R)
=
~
((R
+ p,)2
—
R2)
I =
=
4R2
~ ,=1
p~(1+ p,/2R)2.
N.J. Chernov, GA. Ososkov
/
The usual assumption in any practical application is that p1 << R, otherwise it will be impossible
331
Effective algorithms for circle fitting
Equating to zero the derivatives of this expression
to get a reasonable result. So using (2) we can obtain 2 ~n p~= 4R2L(a, b, R). (7)
by a, b, R we obtain the normal equations, which form a second-order nonlinear system: Fa + Hb ay = Ha + Gb by = Q, (9)
M(a, b, R)
2Pa + 2Qb + .~2= T,
4R
Thus, for LRM it is more “profitable” to minimize not the mean squares L(a, b, R), but R2, whose contribution to the functional (7) is more valuable. This leads to the above-mentioned disadvantage of LRM, i.e. its low stability a superfluous sensitivity to any small errors in measurements. However, it has to be pointed out that the accuracy of the LRM estimations increases rapidly with increasing arc length, and if the points (x,, ~) lie along the full circle the LRM estimations get very close to the true LSM. The method can be generalized to the case of fitting any algebraic curve PN(x, y) = 0, where ~N is a polynomial of the Nth degree. For example, for N = 2 we have the curve equation —
—
—
where as above y = R2 a2 Gauss brackets we denote —
F=
-~-
n
H=~[xy], ~
b2. Using again the
~ [x~ + 3y2]
n
P=~~~[x(x2+y2)J,
[y(x2
!
G=
+ y2]
[3x2
—
+y2)],
~-
!
[(x2
In the system (9) one can exclude any pair of unknown variables to get an equation of the 4th degree. The most suitable choice is to exclude a, b and thereby obtain the equation (10)
y4+Ay3+By2+C-y+D=0 =
a
+y2)2].
2+a 1x
2xy + a3x + a4y + a5.
2=z,x2=u, Theintroductionofnewvariablesy zy = t allows us to obtain the linear regression equation. Then the parameters a 1—a5 are evaluated by solving the corresponding linear normal equation system. Some versions of this approach, leading to a very fast algorithm, will be applied below in no. 5. 4. After the above discussion we consider our own algorithms for circle fitting. The first one is universal, i.e. it is well-suited for fitting any size of circular arc including the full circle. The second is faster, but works only for small size arcs. It will be described in the next section. The idea of the first algorithm is close to the linear regression approach, but it improves the estimation accuracy essentially. Using (7) we shall minimize the functional 2 K(a, b, R)
=
M(a, b, R)R =
+
~
R
—
a
1 Yo
b
)
(8)
1).
=
4+A
3+B 0x
2+C 0x
B0
=
C 0
=
(11) y/y~with the coeffi-
0x+D0=0
related to the variable x
cients A 0 = A/y0,
—
2
2
([x + [~2 It is easy to obtain ~y1~ as a by-product when calculating the coefficiei~itsof (10). Since y~is the first good approximation we can minimize roundoff errors dividing all coefficients of (10) by y~. The new equation is
2-~x, 2—y,
2+b2—R2 R
In spite of the apparent complexity of this method in comparison with PM or LRM we succeeded to find a very fast numerical solution of (10), which now will be briefly described. The eq. (10) has 4 roots; in order to obtain the first approximation to the root we look for, we can use the solution of the normal equations (5)
x a
____
=
H2, C= coefficients T(F+ G)—2(P2 + Q2), T(H2 — with the A = —F— G, BD= = FG — T— FG) + 2(P2G + Q2F) — 4PQH.
C/y~, D0 = D/y.
=
N.J. Chernov, GA. Ososkov
332
/
It is easy to obtain A0 = —4. The evaluations of other coefficients are rather cumbersome, so we shall result limit ourseles to giving here only the final 3
—
Effective algorithms for circle fitting
be rewritten as 2 K(a, b, R)
=
4L(a, b, R).(13)
M(a, b, R)b~
Thus, in order to find the arc parameters in
n~ B
0 ~ 3, 4 — 2n —4n~D0~<2.5n—3.
~
C0
~ 4n,
question,
it
the functional
is sufficient to minimize
~ (x1 b
k(a, b, R)= ~
—
a 2~x,
—
)
For solving (11) any iterative procedure with 1theis value suitable. We root use the of the by 2—5 iterations with the maximal accuracy 10_i3 admissible on the CDC 6500. Calculating y = y 0x we obtain a and b from (9) and 2 + y. (12)
We have been really happy to find that the normal equations for (14) were linear (!) relative to the unknown parameters a, bandy = R2 a2 h2.
R
2aP+2bQ+y([x2]+[y2])=
i~12 + b2
the starting valuegiving x0 = Newton method
=
~/a2 + b
Thus, our first algorithm consists of 5 stages: (i) Transfer the origin of the coordinate system to the centre of gravity of the set (x,,y,); (ii) Calculate the Gauss-bracket values for [x2], [y2], [xy], [x(x2 +y2)], [y(x2 +y2)], [(x2 +
—
a
—
b
R2 2
(14)
—
2 a [x2] + 2 b [ xy] ny = [x2] + [y2]
=
—
T,
(15)
P,
Here we use again the previous notations P,
Q,
T.
(iii) Calculate coefficients for eq. (11) (iv) Solve (11) by the Newton method with the starting value x 0 = 1; (v) Calculate the parameters y, a, b, R. The Fortran-IV program implementing these 5 stages consists of 90 operations. The only external call to the SQRT routine (for calculating R by (12)) is not obligatory because the decision in 2 insteadrules of R, the majority of applications use R so that the SQRT call can be omitted. The iterational method just described will be referred to as the ILRM Iterational Linear Regression Method.
Obtaining a, b, R from (15) is a trivial task, so that we omit the corresponding formulae.
5. The second algorithm for circle fitting is valued only in the case of a small circular arc
point distances from the circle during simulation wasThe taken as 1. comparative—estimate results of parameters obtained by PM, LRM, ILRM and MLRM are
—
(<400).
Let us suppose from the outset that the (x,, y.) points have their centre of gravity located at the origin of the coordinate system and that the OX axis coincides with the straight line drawn by the LSM our points (when these are more or lessthrough equidistant) via theortwo end-points. In such conditions the arc is almost horizontal and clearly Ibi R. Therefore the relation (7) can
We shall call this method the Modified Linear Regression method MLRM. 6. An investigation was made on the CDC-6500 computer in order to compare the characteristics of all the methods As a was model for the experimental data, described. a test of points simulated according to the normal distribution along the circular arc with a = b = 0, R = 1000. The length of the arc was about 400 (from the point (0, 1000) to (400, 100V~),i.e. about 25°).The number of points varied from n = 4 to n = 50. The variance of
summarized in table 1 in the form of their ratios to LSM estimates which were assumed to be the most exact. The variance estimate was calculated as
a2
=
~
—
a)2 +(y,
—
b)2
—
R). (16)
I =
A good approximation for (16) can be derived
N.!. Chernov, G.A. Ososkov
/
Table 1 Comparative—estimate results of parameters obtained by PM. LRM, ILRM and MLRM
Table 3 Number of operations required for execution of each algorithm +,—
a, b, R estimates a2 estimate _____________________________________________ PM 1% 2% LRM 2% 3% IRLM 0.01 + 0.02% 0.001% MLRM 0.05% 0.005%
from (3), (6) and (7) ((x,
“
-2
=
+(y a)24~2(n 1
—
2
— —
2 —
R2)
b) 3)
(17)
I
where a, 1 and ~kare the parameter estimates from MLRM (see also ref. [6]). The computer-time costs for each method were found by the CDC 6500 time-counting routine SECOND. Results in 10-i s are given in table 2. The resulting comparison shows that ILRM is obviously not much slower than PM, so that the effort in improving on PM (see no. 3) is almost useless. One can also note that in the conditions of our test example, the MLRM accuracy decreases sharply when the arc angle increases; MLRM is practically as fast as the fastest LRM and as exact as ILRM. The numbers of operations required for executions of each algorithm are shown in table 3. These numbers were found according to the common rule: addition and subtraction are treated as one operation, multiplication and division as 3, square root as 25 operations. The symbol i in the ILRM row is the iteration number of the Newton-method-solving equation (11). The mean value of i in our test example was 2.5. However, in order to evaluate it in the general case and also to check the stability of ILRM with respect to different point configurations a numeriTable 2 Computer time in
uo—~
n
4
10
50
LSM PM LRM ILRM MLRM
8 1.1 0.7 1.2 0.8
17 2.0 1.5 2.0 1.6
75 8.1 5.9 6.6 6.0
333
Effective algorithms for circle fitting
LSM PM LRM ILRM MLRM
X,/
97n+169 149n+245 13n+29 9n+49 un +8 6n +20 un + 16i + 21 6n + 7i + 42 lln+14 6n+39
[
total
7n+1 I 1 1 I
719n+537 40n+201 29n +93 29n + 37i + 172 29n+156
cal experiment was made. 10000 x 10 simulated points were uniformly distributed in the square (0< x <1,0
.
[1] W.A. Astarev et al., Thesisy I seminara p0 automatlzatzyi nauchnykh issledovanii v jadernoi fizike i smezhnykh oblastjakh (Donish. Dushanbe. 1980) p. 267 (in Russian). [2] G.A. Ososkov, JINR-83-187 Dubna (1983) p. 40 (in Russian). [31B.L. van der Warden, Matematische Statistik (SpringerVerlag, Berlin, 1957). [4] SN. Sokolov and IN. Sum, JINR-D-810. Dubna (1960) (in Russian). [5] I.M. Lipkin et al., Preprint ITEP no. 47. Moscow (1975) (in Russian). [6] D. Hall, UCRL-118544, ~erkeIey,California (1968). [7] N.N. Govorun et al., JINR-1101, Dubna (1962) (in Russian). [8] N.N. Govorun et al., JINR-2036, Dubna (1965) (in Russian).