🚀 go-pugleaf

In article <1992Apr9.003631.1350@magnus.acs.ohio-state.edu>,
mjackson@magnus.acs.ohio-state.edu (Michel T Jackson) writes:

|>I have a problem that has turned out to be harder than expected.  The problem
|>is estimating the parameters of the best-fit circle given a set of
|>modestly noisy points along a short arc.


I would be most interested in any references (solutions?) to problems of this
variety.  

In trying to solve a similar problem (fitting a circle to a set of
points), I developed a
solution which worked rather well for my data sets:  One can always find
the center
and radius of a circle which pases through three "good" data points. 
I determined the "best-fitting" circle by first computing the centers of
all such
circles by considering every combination of 3 of my data points, then
by
finding the centroid of this distribution of centers.  This resulted in
circles which
"seemed" to "look like" they fit the data pretty well.

If the data points which lie on your arc do not exhibit enough curvature, then
this technique may not work well for you. 

Recently there was some discussion regarding techniques for
least-squares fitting
of ellipses to data points.  Below I include a listing of a program to
do this. 

Good luck!

Article 3777 of sci.math.num-analysis:


In article <1992Mar16.112316.15434@kulcs.uucp>, stefan@cs.kuleuven.ac.be
(Stefan Vandewalle) writes:
> I am looking for an algorithm/program to fit an ellipse to a number
> of measured (x,y) points. 
> (I need the length of the axes and the rotation angle.)
> 
> I've tried standard linear least squares to determine the parameters 
> a,b,c,d,e in the formula
>             ax^2 + 2bxy + cy^2 + dx + ey = 1.
> This however gives values (axis-lengths and rotation angle) that 
> depend on the location of the coordinate origin.
> 
> I'd very much appreciate any useful hints or references ! Please
e-mail.
> 
> Stefan Vandewalle.
> stefan@cs.kuleuven.ac.be

Try this one.  Good luck, Daniel Pfenniger
------------------------------------------------------------------------------
      PROGRAM FELLIPSE
C
C     Example Fortran program fitting a quadratic form to a random 
C     set of points {X,Y} 
C     Method : linear least square with orthogonal decomposition 
C
C     Quadratic form:
C
C        (1+a) X^2 + (1-a) Y^2 + c X Y + d X + e Y + f = 0
C
C     NMAX is the maximum number of data points
C     A, B : work arrays
C     X, Y : data points coordinates 
C
C     Daniel Pfenniger, Geneva Observatory, 6/1991
C
      PARAMETER (NMAX=10000)
C
      DIMENSION A(NMAX,5), B(NMAX), X(NMAX), Y(NMAX)
C
      PRINT *, 'Generation of a noisy ellipse:'
      PRINT *, 'Seed for random numbers (e.g. 1234567):'
      READ *, ISEED
C
    5 PRINT *, 'Number of data points (  <',NMAX,'):'
      READ *, N
      IF (N.GT.NMAX .OR. N.LT.1) GOTO 5
      PRINT *, 'Major and minor axis lengths, ellipse tilt in degree:'
      READ *, R1, R2, TILT
      PRINT *, 'Center X0,Y0 ?'
      READ *, X0, Y0
      PRINT *, 'Min & max angle along the ellipse:'
      READ *, ANMIN, ANMAX
      PRINT *, 'Noise amplitude:'
      READ *, AMPL
C 
C     Generation of an inclined ellipse with a superposed noise
C     RAN produces uniform random numbers in [0,1[
C
      PRINT *, 'Generating data points...'
      CO = COSD(TILT)
      SI = SIND(TILT)
      DANGLE = ANMAX-ANMIN
      DO 10 I = 1, N
         ANGLE = ANMIN + DANGLE*RAN(ISEED)
         RR1 = R1 + AMPL*(1.-2.*RAN(ISEED))
         RR2 = R2 + AMPL*(1.-2.*RAN(ISEED))
         XX = RR1*COSD(ANGLE)
         YY = RR2*SIND(ANGLE)
         X(I) = CO*XX - SI*YY + X0
         Y(I) = SI*XX + CO*YY + Y0
 10   CONTINUE
C     
      PRINT *, 'Searching the best ellipse...'
      CALL FITQDR(N,NMAX,A,B,X,Y,KRANK,RNORM)
C
      PRINT * 
      PRINT *, 'RESULTS:'
      PRINT * 
      PRINT *, 'Rank, norm', KRANK, RNORM
C
C     Coefficients of the quadratic form 
C
      AA = 1. + B(1)
      BB = 1. - B(1)
      CC = B(2)
      DD = B(3)
      EE = B(4)
      FF = B(5)
      PRINT *, 'Coefficients  a, b, c, d, e, f:'
      PRINT *, AA, BB, CC, DD, EE, FF
C
C     Some useful quantities for polar plot
C
      DET = 4.*AA*BB - CC*CC
      IF (DET .EQ. 0.) STOP 'Not a quadratic form!'
      XX0 = (CC*EE - 2.*BB*DD) / DET
      YY0 = (CC*DD - 2.*AA*EE) / DET
      PHI = 0.5*ATAND(CC/(AA-BB+1E-32))
C
C     XX0, YY0 : Center of the found ellipse
C     PHI : Angle of the ellipse
C
      PRINT *, 'Center X0, Y0: ',XX0, YY0, '  Angle: ', PHI,' deg'
      AAA = 0.5*(AA+BB + (AA-BB)/COSD(2.*PHI))
      BBB = AA+BB - AAA
      FFF = AA*XX0*XX0+BB*YY0*YY0+CC*XX0*YY0 - FF
      IF (FFF/AAA .LT. 0. .OR. FFF/BBB .LT. 0.) STOP 'Not an ellipse!'
C
C     RR1, RR2 : semi-axes of the ellipse
C
      RR1 = SQRT(FFF/AAA)
      RR2 = SQRT(FFF/BBB)
      PRINT *, 'Semi-axes R1, R2 :', RR1, RR2
      CO = COSD(PHI)
      SI = SIND(PHI)
C
C     The data points and points along the found ellipse are written in 
C     files for later use (plot)
C
      DO 15 I = 1, N
 15      WRITE (1,*) X(I), Y(I)
C
      DO 20 ANGLE = 1., 360.1
         XX = RR1*COSD(ANGLE)
         YY = RR2*SIND(ANGLE)
         XXX = CO*XX - SI*YY + XX0
         YYY = SI*XX + CO*YY + YY0
 20      WRITE (2,*) XXX, YYY
C
      PRINT *, 'Data points written on unit #1'
      PRINT *, 'Ellipse written on unit #2'
C
      END
C-----------------------------------------------------------------------
C
      SUBROUTINE FITQDR(N,NMAX,A,B,X,Y,KRANK,RNORM)
C
C     Fortran subroutine using the linear least squares routine HFTI 
C     finding the parameters  a, c, d, e, f  defining the best 
C     quadratic form 
C
C            (1+a) X^2 + (1-a) Y^2 + c X Y + d X + e Y + f = 0
C
C     passing through a set of data points {Xi,Yi} i=1,..N 
C     The result is invariant by rotation and translation.
C     If N<6, this program finds an exact solution passing through 5
C             points
C     If N>5, it finds the quadratic form which minimizes
C
C            SUM(i=1,N)  [ (1+a)X^2 + (1-a)Y^2 + c X Y + d X + e Y + f ]^2
C
C     N     : Actual number of data points
C     NMAX  : Maximum number of data points
C     NPAR  : Number of parameters to find (5)
C     A, B  : Work arrays
C     X, Y  : Data point coordinates
C     KRANK : Rank of the data matrix A (should be 5)
C     RNORM : Norm of the residual
C
C     TAU   : tolerance for a zero in HFTI
C
C     On outpout  a,c,d,e,f  are in B(1),B(2), ... B(5) respectively
C
C     Daniel Pfenniger, Geneva Observatory, 6/1991
C
      PARAMETER (NPAR=5, TAU=1E-7)
C
      DIMENSION A(NMAX,NPAR), B(NMAX), X(NMAX), Y(NMAX)
C
C     H, G, IP : Work arrays for HFTI
C
      DIMENSION H(NPAR), G(NPAR), IP(NPAR)
C
C     Building the matrices A and B 
C
      DO 10 I = 1, N
        A(I,1) = X(I)**2-Y(I)**2
        A(I,2) = X(I)*Y(I)
        A(I,3) = X(I)
        A(I,4) = Y(I)
        A(I,5) = 1. 
   10   B(I) = -(X(I)**2+Y(I)**2)
C
C     Solving the LS problem: ||A Z - B|| minimum
C     Call HFTI, the result Z is in the first NPAR rows of B
C     TAU is the tolerance for zero
C     KRANK is the rank of A (should be NPAR)
C     RNORM is the norm of the residual
C     See Lawson & Hanson (1974) for detail
      CALL HFTI(A,NMAX,N,NPAR,B,NMAX,1,TAU,KRANK,RNORM,H,G,IP)
C
      RETURN
      END
C
C     Some LSQ routines from Lawson & Hanson 
C
      SUBROUTINE HFTI (A,MDA,M,N,B,MDB,NB,TAU,KRANK,RNORM,H,G,IP)
C     C.L.LAWSON & R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12
C     TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL 1974
C     SOLVE LEAST SQUARES PROBLEM USING ALGORITHM HFTI.
C
      DIMENSION A(MDA,N),B(MDB,NB),H(N),G(N),RNORM(NB)
      INTEGER IP(N)
      DOUBLE PRECISION SM,DZERO
      DATA SZERO/ 0./, DZERO / 0.D0 /
      DATA FACTOR / .001 /
C
      K = 0
      LDIAG = MIN0(M,N)
      IF (LDIAG.LE.0) GOTO 270
        DO 80 J=1,LDIAG
          IF (J.EQ.1) GOTO 20
C
C     UPDATE SQUARED COLUMN LENGTHS AND FIND LMAX
C
          LMAX = J
            DO 10 L = J,N
              H(L) = H(L) - A(J-1,L)**2
              IF (H(L).GT.H(LMAX)) LMAX = L
   10       CONTINUE
          IF (DIFF(HMAX+FACTOR*H(LMAX),HMAX)) 20,20,50
C
C     COMPUTE SQUARED COLUMN LENGTHS AND FIND LMAX
C
   20     LMAX = J
            DO 40 L = J,N
              H(L) = SZERO
              DO 30 I = J,M
   30           H(L) = H(L) + A(I,L)**2
              IF (H(L).GT.H(LMAX)) LMAX = L
   40       CONTINUE
          HMAX = H(LMAX)
C
C     LMAX HAS BEEN DETERMINED
C
C     DO COLUMN INTERCHANGES IF NEEDED.
C
   50     CONTINUE
          IP(J) = LMAX
          IF (IP(J).EQ.J) GOTO 70
          DO 60 I = 1,M
            TMP = A(I,J)
            A(I,J) = A(I,LMAX)
   60       A(I,LMAX) = TMP
          H(LMAX) = H(J)
C
C     COMPUTE THE J-TH TRANSFORMATION AND APPLY IT TO A AND B.
C
   70     CALL H12 (1,J,J+1,M,A(1,J),1,H(J),A(1,J+1),1,MDA,N-J)
   80     CALL H12 (2,J,J+1,M,A(1,J),1,H(J),B,1,MDB,NB)
C
C     DETERMINE THE PSEUDORANK, K, USING THE TOLERANCE, TAU.
C
          DO 90 J = 1,LDIAG
            IF (ABS(A(J,J)).LE.TAU) GOTO 100
   90     CONTINUE
      K = LDIAG
      GOTO 110
  100 K = J - 1
  110 KP1 = K + 1
C
C     COMPUTE THE NORMS OF THE RESIDUAL VECTORS.
C
      IF (NB.LE.0) GOTO 140
        DO 130 JB = 1,NB
          TMP = SZERO
          IF (KP1.GT.M) GOTO 130
          DO 120 I = KP1,M
  120       TMP = TMP + B(I,JB)**2
  130     RNORM(JB) = SQRT(TMP)
  140 CONTINUE
C     SPECIAL FOR PSEUDORANK = 0
      IF (K.GT.0) GOTO 160
      IF (NB.LE.0) GOTO 270
        DO 150 JB = 1,NB
          DO 150 I = 1,N
  150       B(I,JB) = SZERO
      GOTO 270
C
C     IF THE PSEUDORANK IS LESS THAN N COMPUTE HOUSHOLDER
C     DECOMPOSITION OF FIRST K ROWS.
C
  160 IF (K.EQ.N) GOTO 180
        DO 170 II = 1,K
          I = KP1 - II
  170     CALL H12 (1,I,KP1,N,A(I,1),MDA,G(I),A,MDA,1,I-1)
  180 CONTINUE
C
C
      IF (NB.LE.0) GOTO 270
        DO 260 JB = 1,NB
C
C     SOLVE THE K BY K TRIANGULAR SYSTEM
C
          DO 210 L = 1,K
            SM = DZERO
            I = KP1 - L
            IF (I.EQ.K) GOTO 200
            IP1 = I + 1
            DO 190 J = IP1,K
  190         SM = SM + A(I,J)*DBLE(B(J,JB))
  200       SM1 = SM
  210       B(I,JB) = (B(I,JB)-SM1)/A(I,I)
C
          IF (K.EQ.N) GOTO 240
            DO 220 J = KP1,N
  220         B(J,JB) = SZERO
            DO 230 I = 1,K
  230         CALL H12 (2,I,KP1,N,A(I,1),MDA,G(I),B(1,JB),1,MDB,1)
C
C     RE-ORDER THE SOLUTION VECTOR TO COMPENSATE FOR THE
C     COLUMN INTERCHANGES.
C
C     COMPLETE COMPUTATION OF SOLUTION VECTOR
C
  240     DO 250 JJ = 1,LDIAG
            J = LDIAG + 1 - JJ
            IF (IP(J).EQ.J) GOTO 250
            L = IP(J)
            TMP = B(L,JB)
            B(L,JB) = B(J,JB)
            B(J,JB) = TMP
  250     CONTINUE
  260   CONTINUE
C
C     THE SOLUTION VECTORS, X, ARE NOW
C     IN THE FIRST N ROWS OF THE ARRAY B(.).
C
  270 KRANK = K
      RETURN
      END
C-----------------------------------------------------------------------
C
C     SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV)
C     C.L.LAWSON & R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12
C     TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974
C
C     CONSTRUCTION AND/OR APPLICATION OF A SINGLE
C     HOUSHOLDER TRANSFORMATION:  Q = I + U*(U**T)/B
C
C     MODE   = 1 OR 2  TO SELECT ALGORITHM H1 OR H2
C     LPIVOT IS THE INDEX OF THE PIVOT ELEMENT.
C     L1,M   IF L1 .LE. M  THE TRANSFORMATION WILL BE CONSTRUCTED TO
C            ZERO ELEMENTS INDEXED FROM L1 THROUGH M.  IF L1 .GT. M
C            THE SUBROUTINE DOES AN IDENTITY TRANSFORMATION.
C     U(),IUE,UP    ON ENTRY TO H1 U() CONTAINS THE PIVOT VECTOR.
C                   IUE IS THE STORAGE INCREMENT BETWEEN ELEMENTS.
C                   ON EXIT FROM H1 U() AND UP
C                   CONTAIN QUANTITES DEFINING THE VECTOR U OF THE
C                   HOUSHOLDER TRANSFORMATION.  ON ENTRY TO H2 U()
C                   AND UP SHOULD CONTAIN QUATITIES PREVIOUSLY COMPUTED
C                   BY H1. THESE WILL NOT BE MODIFIED BY H2.
C     C()    ON ENTRY TO H1 OR H2 C() CONTAINS A MATRIX WHICH WILL BE
C            REGARDED AS A SET OF VECTORS TO WHICH THE HOUSHOLDER 
C            TRANSFORMATION IS TO BE APPLIED. ON EXIT C() CONTAINS THE
C            SET OF TRANSFORMED VECTORS.
C     ICE    STORAGE INCREMENT BETWEEN ELEMENTS OF VECTORS IN C().
C     ICV    STORAGE INCREMENT BETWEEN VECTORS IN C().
C     NCV    NUMBER OF VECTORS IN C() TO BE TRANSFORMED. IF NCV .LE. 0
C            NO OPERATIONS WILL BE DONE ON C().
C
      SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV)
      DIMENSION U(IUE,M),C(1)
      DOUBLE PRECISION SM,B
      DATA ONE / 1. /
C
      IF (0.GE.LPIVOT.OR.LPIVOT.GE.L1.OR.L1.GT.M) RETURN
      CL = ABS(U(1,LPIVOT))
      IF (MODE.EQ.2) GOTO 60
C     *** CONSTRUCT THE TRANSFORMATION ***
      DO 10 J = L1,M
   10   CL = AMAX1(ABS(U(1,J)),CL)
      IF (CL) 130,130,20
   20 CLINV = ONE/CL
      SM = (DBLE(U(1,LPIVOT))*CLINV)**2
      DO 30 J = L1,M
   30   SM = SM + (DBLE(U(1,J))*CLINV)**2
C       CONVERT DBLE PREC SM TO SNGL. PREC. SM1
      SM1 = SM
      CL = CL*SQRT(SM1)
      IF (U(1,LPIVOT)) 50,50,40
   40 CL = -CL
   50 UP = U(1,LPIVOT) - CL
      U(1,LPIVOT) = CL
      GOTO 70
C     *** APPLY THE TRANSFORMATION I + U*(U**T)*B TO C ***
C
   60 IF (CL) 130,130,70
   70 IF (NCV.LE.0) RETURN
      B = DBLE(UP)*U(1,LPIVOT)
C       B MUST BE NONPOSITIVE HERE. IF B=0.,RETURN
C
      IF (B) 80,130,130
   80 B = ONE/B
      I2 = 1 - ICV + ICE*(LPIVOT-1)
      INCR = ICE*(L1-LPIVOT)
      DO 120 J = 1,NCV
        I2 = I2 +ICV
        I3 = I2 + INCR
        I4 = I3
        SM = C(I2)*DBLE(UP)
        DO 90 I = L1,M
          SM = SM + C(I3)*DBLE(U(1,I))
   90     I3 = I3 + ICE
          IF (SM) 100,120,100
  100     SM = SM*B
          C(I2) = C(I2) + SM*DBLE(UP)
          DO 110 I = L1,M
            C(I4) = C(I4) + SM*DBLE(U(1,I))
  110       I4 = I4 + ICE
  120 CONTINUE
  130 RETURN
      END
C--------------------------- this routine is not a joke!
C
      FUNCTION DIFF(X,Y)
C     C.L.LAWSON & R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUNE 7
C     TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974
      DIFF = X-Y
      RETURN
      END