flash@lehi3b15.csee.Lehigh.EDU (Stephen Corbesero) (06/01/89)
I need an algorithm to fit a set of data points, obtained experimentally, to the best-fit circle. I have searched through many book son least squares anaylsis, but can only find simple cases. I started to do the math, but it became rather involved very quickly. I would appreciate it if i can get code to do this, or perhaps leads as to where I can continue looking. Thanks in advance.... -- Stephen Corbesero flash@lehi3b15.UUCP VLSI Design Automation Lab flash@lehi3b15.CSEE.Lehigh.EDU Computer Science And Electrical Engineering, usgcorb@vax1.CC.Lehigh.EDU Lehigh University, Bethlehem, PA 18015
flowers@osf.OSF.ORG (Ken Flowers) (06/02/89)
In article <573@lehi3b15.csee.Lehigh.EDU> flash@lehi3b15.csee.Lehigh.EDU (Stephen Corbesero) writes: > >I need an algorithm to fit a set of data points, obtained >experimentally, to the best-fit circle. I have searched through many >book son least squares anaylsis, but can only find simple cases. I >started to do the math, but it became rather involved very quickly. > >I would appreciate it if i can get code to do this, or perhaps leads >as to where I can continue looking. > >Thanks in advance.... > >-- >Stephen Corbesero flash@lehi3b15.UUCP >VLSI Design Automation Lab flash@lehi3b15.CSEE.Lehigh.EDU >Computer Science And Electrical Engineering, usgcorb@vax1.CC.Lehigh.EDU >Lehigh University, Bethlehem, PA 18015 Not being a mathmatician, (but I did start a math major at MIT. Finally went on to MechE.) please excuse any level of ignorance. It seems to me that the best fit circle can be simply found by doing the following: 1. Make the center the average point of all the data points. 2. Make the radius the average distance from the center to each data point. This gives also the smallest best fit circle, because as we can see from the two point generalization that an infinate number of circles can fit exactly. Also, if we take a radius significantly larger than the range of the data and then place the data around a small segment of the resulting circle, that seems like a preaty damn good fit also. What it turns out were doing in this case, if we make the radius suffiently large, is creating the best fit line.
flynn@pixel.cps.msu.edu (Patrick J. Flynn) (06/02/89)
In article <1088@osf.OSF.ORG> flowers@osf.org (Ken Flowers) writes: >In article <573@lehi3b15.csee.Lehigh.EDU> flash@lehi3b15.csee.Lehigh.EDU (Stephen Corbesero) writes: >> >>I need an algorithm to fit a set of data points, obtained >>experimentally, to the best-fit circle. I have searched through many >>book son least squares anaylsis, but can only find simple cases. I >>started to do the math, but it became rather involved very quickly. >> >It seems to me that the best fit circle can be simply found by doing >the following: > > 1. Make the center the average point of all the data points. > > 2. Make the radius the average distance from the center to each > data point. Unfortunately, if the points are not distributed uniformly along the *entire* circumference, you'll get a biased estimate of the location of the center. A least-squares error criterion to the problem is nonlinear in the `geometric' parameters of the circle (center coordinates and radius), but if you're willing to use iterative techniques, you can get a solution. Minimize f^2(x,y;x0,y0,r) with respect to (x0,y0) and r, where f(x,y;x0,y0,r) = (x-x0)^2 + (y-y0)^2 - r^2 I've used the Levenberg-Marquardt method to solve related problems in the past (you have to watch out for local minima, though). A second approach which avoids the nonlinear stuff: Start choosing triples of points from your point set. Find the circle passing through those points; save the estimates of the center and radius. Then average the estimates. If the points are not very noisy, that might do; otherwise, you might want to trim those estimates. When choosing triples, be careful to choose points which aren't too close together. ------------------------------------------------------------------------------ Pat Flynn, CS, Mich. State U | "Don't eat with your hands, son; use your flynn@cps.msu.edu | entrenching tool." (517) 353-4638 | -Firesign Theatre
toms@ncifcrf.gov (Tom Schneider) (06/02/89)
In article <1088@osf.OSF.ORG> flowers@osf.org (Ken Flowers) writes: >In article <573@lehi3b15.csee.Lehigh.EDU> flash@lehi3b15.csee.Lehigh.EDU (Stephen Corbesero) writes: >> >>I need an algorithm to fit a set of data points, obtained >>experimentally, to the best-fit circle. I have searched through many >>book son least squares anaylsis, but can only find simple cases. I >>started to do the math, but it became rather involved very quickly. >> >>Stephen Corbesero flash@lehi3b15.UUCP >>VLSI Design Automation Lab flash@lehi3b15.CSEE.Lehigh.EDU >>Computer Science And Electrical Engineering, usgcorb@vax1.CC.Lehigh.EDU >>Lehigh University, Bethlehem, PA 18015 > 1. Make the center the average point of all the data points. > 2. Make the radius the average distance from the center to each data point. That was my first thought also. But suppose that all of the data points lie on one side of the circle, or are confined to a small angular region. Then this method will fail because it will place the center much to close to the circle. If Stephen knows that the points are pretty well spread around then this method should be pretty good, but it wouldn't be as good as a true least squares fit. Seems to me that should be possible, where one is determining the minimum square radial distance from each point to the circle along radii. (Also an MIT grad, but in biology!) Tom Schneider National Cancer Institute Laboratory of Mathematical Biology Frederick, Maryland 21701-1013 toms@ncifcrf.gov
roy@phri.UUCP (Roy Smith) (06/02/89)
flash@lehi3b15.csee.Lehigh.EDU (Stephen Corbesero) writes: > I need an algorithm to fit a set of data points, obtained > experimentally, to the best-fit circle. There was an article about finding the smallest circle which just encloses a set of points on a plane recently on the net which might be of some interest to you. It appeared in rec.humor.funny. Had something to do with sheep. -- Roy Smith, Public Health Research Institute 455 First Avenue, New York, NY 10016 {allegra,philabs,cmcl2,rutgers,hombre}!phri!roy -or- roy@alanine.phri.nyu.edu "The connector is the network"
bill@ut-emx.UUCP (Bill Jefferys) (06/02/89)
In article <573@lehi3b15.csee.Lehigh.EDU> flash@lehi3b15.csee.Lehigh.EDU (Stephen Corbesero) writes:
#
#I need an algorithm to fit a set of data points, obtained
#experimentally, to the best-fit circle. I have searched through many
#book son least squares anaylsis, but can only find simple cases. I
#started to do the math, but it became rather involved very quickly.
Look in Wayne A. Fuller's _Measurement Error Models_, Chapter 3
(esp. Section 3.2). Published by Wiley, 1987.
Warning: This is not easy going.
Bill Jefferys
#
#I would appreciate it if i can get code to do this, or perhaps leads
#as to where I can continue looking.
#
#Thanks in advance....
#
#--
#Stephen Corbesero flash@lehi3b15.UUCP
#VLSI Design Automation Lab flash@lehi3b15.CSEE.Lehigh.EDU
#Computer Science And Electrical Engineering, usgcorb@vax1.CC.Lehigh.EDU
#Lehigh University, Bethlehem, PA 18015ted@nmsu.edu (Ted Dunning) (06/03/89)
interestingly, while the coordinates of the center of the best fit circle involves a non-linear solution, if you are given the center, then the best fit radius is the root mean square of the distances from the center. using this fact will make a non-linear optimization program for finding the best fit circle converge much more rapidly.
pfenniger@obs.unige.ch (06/08/89)
In article <573@lehi3b15.csee.Lehigh.EDU>, flash@lehi3b15.csee.Lehigh.EDU (Stephen Corbesero) writes: > I need an algorithm to fit a set of data points, obtained > experimentally, to the best-fit circle. I have searched through many > book son least squares anaylsis, but can only find simple cases. I > started to do the math, but it became rather involved very quickly. > The problem can be stated as follow : You have a set of points in cartesian coordinates {Xi, Yi} i = 1, ..., n (n>2). You want to fit to these points in a least square sense the quadratic form of a circle of unknown radius R and center {X0, Y0} : f(X,Y) = (X - X0)^2 + (Y - Y0)^2 - R^2 = a X + b Y + c + X^2 + Y^2 , where a = -2 X0, b = -2 Y0, c = X0^2 + Y0^2 - R^2. In matrix form the problem is to find Z which minimizes the euclidian norm : || A Z - B || where [ a ] Z = [ b ] [ c ] and [ X1 Y1 1 ] [ X1^2 + Y1^2 ] A = [ X2 Y2 1 ] , B = -[ X2^2 + Y2^2 ] . ....... ........... [ Xn Yn 1 ] [ Xn^2 + Yn^2 ] The problem is therefore suited to *linear* least squares, since the unknowns a, b, c in Z enter in a linear way. Then use any modern linear least squares algorithm (e.g. HFTI in Lawson & Hanson, 1974, "Solving Least Squares Problems", Prentice-Hall, NJ ; or Linpack etc.) which just finds this minimum euclidian norm by orthogonal decomposition. This approach is superior to the traditional inversion of a square normal matrix. Note that fitting ellipses or more complicated forms would use essentially the same method, since the additional unknown parameters are also linear. Daniel Pfenniger
wjmartiniii@violet.waterloo.edu (Bill Martin) (06/09/89)
In article <1088@osf.OSF.ORG> flowers@osf.org (Ken Flowers) writes: >In article <573@lehi3b15.csee.Lehigh.EDU> flash@lehi3b15.csee.Lehigh.EDU (Stephen Corbesero) writes: >> >>I need an algorithm to fit a set of data points, obtained >>experimentally, to the best-fit circle. I have searched through many >>-- >It seems to me that the best fit circle can be simply found by doing >the following: > > 1. Make the center the average point of all the data points. > > 2. Make the radius the average distance from the center to each > data point. > With regard to the above approach, consider a large number of points lying on a line. The given algorithm would choose a point on that line as the center of the circle. This is certainly not optimal in general.
finn@sunshine.cad.mcc.com (Chris Finn) (06/16/89)
In article <221@obs.unige.ch> pfenniger@obs.unige.ch writes: >In article <573@lehi3b15.csee.Lehigh.EDU>, flash@lehi3b15.csee.Lehigh.EDU (Stephen Corbesero) writes: >> I need an algorithm to fit a set of data points, obtained >> experimentally, to the best-fit circle. I have searched through many >> book son least squares anaylsis, but can only find simple cases. I >> started to do the math, but it became rather involved very quickly. >> > >The problem can be stated as follow : >You have a set of points in cartesian coordinates {Xi, Yi} i = 1, ..., n (n>2). >You want to fit to these points in a least square sense the quadratic form of a >circle of unknown radius R and center {X0, Y0} : > > f(X,Y) = (X - X0)^2 + (Y - Y0)^2 - R^2 > > = a X + b Y + c + X^2 + Y^2 , > >where a = -2 X0, b = -2 Y0, c = X0^2 + Y0^2 - R^2. matrix notation stuff deleted >The problem is therefore suited to *linear* least squares, since the unknowns >a, b, c in Z enter in a linear way. >Then use any modern linear least squares algorithm (e.g. HFTI in Lawson & >Hanson, 1974, "Solving Least Squares Problems", Prentice-Hall, NJ ; or Linpack >etc.) which just finds this minimum euclidian norm by orthogonal decomposition. >This approach is superior to the traditional inversion of a square normal >matrix. > I found myself having to solve this very problem after following the discussion of it in comp.graphics. Since I had not seen anyone post the actually code for doing this I thought I'd send in mine. I followed Mr. Pfenniger's suggestion for linearizing the problem. However, I did not use orthogonal decomposition but rather inverted the square normal matrix ( by Cramer's Rule no less, well it's only 3x3). I think it should be fairly robust if a large number of angles are present in the data. Included is a driver program to generate a circle and test the subroutine which actually performs the fit. Chris Finn MCC CAD Program [512] 338-3637 ARPA: finn@mcc.com UUCP: {uunet,harvard,gatech,pyramid}!cs.utexas.edu!milano!cadillac!finn ----------------------------- Cut Here -------------------------- PROGRAM CIRCLE C PARAMETER (MAXDATA=10000) C REAL X(MAXDATA),Y(MAXDATA) REAL X0,Y0,R,X0E,Y0E,RE INTEGER NDATA C R = 2000.0 X0 = 150.0 Y0 = -300.0 NDATA = 500 C C Fill up data arrays with actual coordinates of a circle C DO 10 I = 1, NDATA/2 X(I) = (X0-R) + (I-1)*2.0*R/FLOAT(NDATA/2) Y(I) = SQRT( R*R - (X(I)-X0)*(X(I)-X0) ) + Y0 10 CONTINUE DO 15 I = NDATA/2 + 1, NDATA J = I - NDATA/2 X(I) = (X0+R) - (J-1)*2.0*R/FLOAT(NDATA/2) Y(I) = -1.0*SQRT( R*R - (X(I)-X0)*(X(I)-X0) ) + Y0 15 CONTINUE C CALL CIRCFIT(X,Y,NDATA,X0E,Y0E,RE) C PRINT *, ' Estimated circle center (x0,y0) =', X0E,Y0E PRINT *, ' Estimated circle radius =', RE C STOP END C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SUBROUTINE CIRCFIT(X,Y,N,X0,Y0,R) C C This subroutine fits a circle to a set of data points C Input : X - array containing x-coordinates of data values C Y - array containing y-coordinates of data values C N - number of X,Y pairs C Output: X0,Y0 - estimated x,y coordinates of circles center C R - estimated circle radius C C Explaination: C The equation of a circle is: C (X-X0)**2 + (Y-Y0)**2 = R**2 (1) C which can be written C X**2 + Y**2 - A*X - B*Y + C = 0 (2) C where C A = 2*X0, B = 2*Y0, and C = X0**2 + Y0**2 - R**2. (3) C A least squares fit of equation (2) to a set of data C (x_i,y_i,i=1,N) is obtained by minimizing C C N C SUM( x_i**2 + y_i**2 - A*x_i - B*y_i + C )**2. (4) C i=1 C C with respect to A,B, and C. This is accomplished by taking C the derivatives of (4) with respect to A, B, and C and setting C them equal to zero. C The resultant set of 3 linear equations in three C unknowns is solved using Cramer's Rule. C REAL X(N),Y(N),X0,Y0,R INTEGER N REAL A,B,C,G(3,3),D(3) REAL XSQ,YSQ,XQ,YQ,DETG C DO 5 I = 1, 3 D(I) = 0.0 DO 6 J = 1, 3 G(I,J) = 0.0 6 CONTINUE 5 CONTINUE C C Construct coefficient matrix G and right hand side C vector D C DO 10 I = 1, N XSQ = X(I)*X(I) XQ = XSQ*X(I) YSQ = Y(I)*Y(I) YQ = YSQ*Y(I) G(1,1) = G(1,1) + XSQ G(1,2) = G(1,2) + X(I)*Y(I) G(1,3) = G(1,3) + X(I) G(2,2) = G(2,2) + YSQ G(2,3) = G(2,3) + Y(I) D(1) = D(1) + XQ + X(I)*YSQ D(2) = D(2) + XSQ*Y(I) + YQ D(3) = D(3) - XSQ - YSQ 10 CONTINUE C G(1,3) = -G(1,3) G(2,3) = -G(2,3) C C Take advantage of the symmetry of the coefficient matrix C G(2,1) = G(1,2) G(3,1) = G(1,3) G(3,2) = G(2,3) G(3,3) = FLOAT(N) C C Compute the determinant of the coefficient matrix C DETG = G(1,1)*( G(2,2)*G(3,3) - G(2,3)*G(3,2) ) C - G(1,2)*( G(2,1)*G(3,3) - G(2,3)*G(3,1) ) C + G(1,3)*( G(2,1)*G(3,2) - G(2,2)*G(3,1) ) C C Return if matrix is singular C IF (DETG .EQ. 0.0 ) THEN PRINT *, ' YOU IN A WORLD A HURT!' RETURN ENDIF C C Solve for the unknown parameters of the circle using C Cramer's Rule C A = D(1) *( G(2,2)*G(3,3) - G(2,3)*G(3,2) ) C - G(1,2)*( D(2) *G(3,3) - G(2,3)*D(3) ) C + G(1,3)*( D(2) *G(3,2) - G(2,2)*D(3) ) C B = G(1,1)*( D(2) *G(3,3) - G(2,3)*D(3) ) C - D(1) *( G(2,1)*G(3,3) - G(2,3)*G(3,1) ) C + G(1,3)*( G(2,1)*D(3) - D(2) *G(3,1) ) C C = G(1,1)*( G(2,2)*D(3) - D(2) *G(3,2) ) C - G(1,2)*( G(2,1)*D(3) - D(2) *G(3,1) ) C + D(1) *( G(2,1)*G(3,2) - G(2,2)*G(3,1) ) C A = A/DETG B = B/DETG C = C/DETG C C Unscramble the coefficients used for the solution C into the actual parameters of a circle C X0 = A/2.0 Y0 = B/2.0 R = SQRT(X0*X0 + Y0*Y0 - C) C RETURN END
finn@sunshine.cad.mcc.com (Chris Finn) (06/20/89)
> Good job Chris! How accurate are the results? The reason I ask is because > of the inversion step. Running the example I posted outputs Estimated circle center (x0,y0) = 149.9992 -299.9999 Estimated circle radius = 2000.001 The values used to generate the circle were (x0,y0) = 150.0 -300.0 and radius = 2000.0. The results are accurate to within the precision of the machine ( approximately six or seven decimal places ). These results are for a set of 500 data points which are an "exact" circle evenly sampled along its circumference. As with any inverse problem the quality of the results are highly dependent on the quality of the data input. Data which lie along a circle but contain noise will be fit in the least squares sense (i.e. the summed squared distance between the predicted circle and the data will be minimized). For reasonable data sets the matrix inversion should be stable. A more fundamental problem is fitting data which do not fully determine the parameters of the circle. An example which comes to mind is to fit a circle to a point. As a more likely example consider a small data set which covers a very small segment of the circles arc ( say five degrees ) and contains some noise. Such a data set could be appear as a straight line requiring a circle of infinite radius with a center at infinity. For those of you out there in netland who solve a lot of inverse problems I highly recommend "Inverse Problem Theory: Methods for Data Fitting and Model Parameter Estiation" by Albert Tarantola, Elsevier Science Publishers. It is very complete and contains lots of examples for special cases. Chris Finn MCC CAD Program, P.O. Box 200195, Austin, TX 78720 [512] 343-0978 ARPA: finn@mcc.com UUCP: {uunet,harvard,gatech,pyramid}!cs.utexas.edu!milano!cadillac!finn
hutch@celerity.uucp (Jim Hutchison) (07/23/89)
finn@MCC.COM (Chris Finn) writes: >pfenniger@obs.unige.ch writes: >>flash@lehi3b15.csee.Lehigh.EDU (Stephen Corbesero) writes: >>> I need an algorithm to fit a set of data points, obtained >>> experimentally, to the best-fit circle. I have searched through many >>> book son least squares anaylsis, but can only find simple cases. ... >>The problem can be stated as follow : >>You have a set of points in cartesian coordinates {Xi, Yi} i = 1, ..., n (n>2). [math problem/?simple case explanation? deleted] >>The problem is therefore suited to *linear* least squares, since the unknowns >>a, b, c in Z enter in a linear way. There is atleast 1 other possibility. It could be reduced to a simpler problem. If the problem where instead fitting a circle around an arbitrary polygon, it'd be much simpler right? 1) Take the first three arbitrary points, construct a triangle. 2) Any points which when added to the polygon make a concave edge, are not added to the polygon. Just drop them on the floor. By concave edge, I mean an edge formed by adding a point inside the convex polygon. 3) Repeat step 2 for all points in your set. 4) Fit a circle to the polygon. This may not be the ideal algorithm. For a large number of points it should atleast reduces the size of the problem with no loss of precision. >[...] I followed Mr. >Pfenniger's suggestion for linearizing the problem. However, I did not use >orthogonal decomposition but rather inverted the square normal matrix ( by >Cramer's Rule no less, well it's only 3x3). I think it should be fairly >robust if a large number of angles are present in the data. Good job Chris! How accurate are the results? The reason I ask is because of the inversion step. /* Jim Hutchison {dcdwest,ucbvax}!ucsd!celerity!hutch */ /* Disclaimor: I am not an official spokesman for FPS computing */