Xref: utzoo comp.graphics:6167 sci.math:7040 Path: utzoo!utgpu!jarvis.csri.toronto.edu!rutgers!cs.utexas.edu!milano!cadillac!sunshine!finn From: finn@sunshine.cad.mcc.com (Chris Finn) Newsgroups: comp.graphics,sci.math Subject: Re: Need to fit a circle to some points Keywords: least squares Message-ID: <1241@cadillac.CAD.MCC.COM> Date: 15 Jun 89 19:49:57 GMT References: <573@lehi3b15.csee.Lehigh.EDU> <221@obs.unige.ch> Sender: news@cadillac.CAD.MCC.COM Reply-To: finn@MCC.COM (Chris Finn) Organization: MCC CAD Program, Austin, TX Lines: 190 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