AstroPy - fitcircle.py


Description

Fit a circle (cx, cy, radius) to a set of points using an iterative, converging, method that does not use derivatives. Uses a 3d grid of candidates and is deterministic, so it gives the same results every time. Can fit a circle to an arbitrary precision (although this does not mean that the error will reach zero). This program is used in the alignment programs that follow.

Usage

fitcircle.py x1 y1 x2 y2 x3 y3 [x4 y4 ...]

The program contains a function which can be called from other programs:

fitCirc(pts, maxErr=0.001, maxIters=100, dz=10, fact=0.5, printout=False)
   The function exits when either maxErr or maxIters have been reached.
   dz is the initial change in coordinates.
   fact is the change in dz whenever the error does not change.

Examples

When run from the command line, this program fits a circle to a set of input points. The initial circle is at the center of the points, with the average radius. The iterations will stop after a fixed number, or when the error falls below a set value. At each step the program prints out the index of the best grid candidate (0-26, 13 at center), cx, cy, radius, error, and current grid size. The following arguments are 4 points on a semi-circle at cx=123.45, cy=54.321, rad=67.89:

> fitcircle.py 191.340 54.321 175.457 97.960 135.239 121.180 89.505 113.115
Initial: cx=147.885250      cy=96.644000       rad=44.131258       err=16.528032      

iter idx cx              cy              rad             err             dz
---- --- --------------- --------------- --------------- --------------- ---------------
1    2   127.885250      76.644000       54.131258       8.043911        10.000000      
2    11  127.885250      66.644000       64.131258       6.113788        10.000000      
3    10  127.885250      56.644000       64.131258       2.210322        10.000000      
4    2   117.885250      46.644000       74.131258       2.133990        10.000000      
5    13  117.885250      46.644000       74.131258       2.133990        10.000000      
6    24  122.885250      51.644000       69.131258       1.019784        5.000000       
7    13  122.885250      51.644000       69.131258       1.019784        5.000000       
8    16  122.885250      54.144000       69.131258       0.967667        2.500000       
9    13  122.885250      54.144000       69.131258       0.967667        2.500000       
10   10  122.885250      52.894000       69.131258       0.386292        1.250000       
11   13  122.885250      52.894000       69.131258       0.386292        1.250000       
12   15  122.885250      53.519000       68.506258       0.238218        0.625000       
13   24  123.510250      54.144000       67.881258       0.137164        0.625000       
14   13  123.510250      54.144000       67.881258       0.137164        0.625000       
15   16  123.510250      54.456500       67.881258       0.102990        0.312500       
16   13  123.510250      54.456500       67.881258       0.102990        0.312500       
17   10  123.510250      54.300250       67.881258       0.041142        0.156250       
18   13  123.510250      54.300250       67.881258       0.041142        0.156250       
19   4   123.432125      54.300250       67.881258       0.029063        0.078125       
20   13  123.432125      54.300250       67.881258       0.029063        0.078125       
21   14  123.432125      54.300250       67.920321       0.012989        0.039062       
22   13  123.432125      54.300250       67.920321       0.012989        0.039062       
23   21  123.451656      54.300250       67.900790       0.008808        0.019531       
24   13  123.451656      54.300250       67.900790       0.008808        0.019531       
25   7   123.441891      54.310016       67.900790       0.003573        0.009766       
26   24  123.451656      54.319781       67.891024       0.001499        0.009766       
27   13  123.451656      54.319781       67.891024       0.001499        0.009766       
28   15  123.451656      54.324664       67.886141       0.001273        0.004883       
29   13  123.451656      54.324664       67.886141       0.001273        0.004883       
30   11  123.451656      54.322223       67.888583       0.000476        0.002441       

Final: cx=123.451656      cy=54.322223       rad=67.888583       err=0.000476        npts=4  

The following arguments are 4 points on the same semi-circle, but perturbed by random values of +-0.5:

> fitcircle.py 191.794 54.352 175.833 98.121 135.317 120.707 89.938 113.330
Initial: cx=148.220500      cy=96.627500       rad=44.077921       err=16.592483      

iter idx cx              cy              rad             err             dz
---- --- --------------- --------------- --------------- --------------- ---------------
1    2   128.220500      76.627500       54.077921       8.214430        10.000000      
2    11  128.220500      66.627500       64.077921       6.202184        10.000000      
3    10  128.220500      56.627500       64.077921       2.238150        10.000000      
4    2   118.220500      46.627500       74.077921       2.040214        10.000000      
5    13  118.220500      46.627500       74.077921       2.040214        10.000000      
6    24  123.220500      51.627500       69.077921       0.955905        5.000000       
7    13  123.220500      51.627500       69.077921       0.955905        5.000000       
8    13  123.220500      51.627500       69.077921       0.955905        2.500000       
9    16  123.220500      52.877500       69.077921       0.305358        1.250000       
10   13  123.220500      52.877500       69.077921       0.305358        1.250000       
...
90   13  123.527117      53.637896       68.381259       0.215035        0.000000       
91   13  123.527117      53.637896       68.381259       0.215035        0.000000       
92   12  123.527117      53.637896       68.381259       0.215035        0.000000       
93   13  123.527117      53.637896       68.381259       0.215035        0.000000       
94   12  123.527117      53.637896       68.381259       0.215035        0.000000       
95   13  123.527117      53.637896       68.381259       0.215035        0.000000       
96   13  123.527117      53.637896       68.381259       0.215035        0.000000       
97   13  123.527117      53.637896       68.381259       0.215035        0.000000       
98   13  123.527117      53.637896       68.381259       0.215035        0.000000       
99   13  123.527117      53.637896       68.381259       0.215035        0.000000       
100  12  123.527117      53.637896       68.381259       0.215035        0.000000       

Final: cx=123.527117      cy=53.637896       rad=68.381259       err=0.215035        npts=4  

In this case the error doesn't fall to zero, although it stops changing (to 6 places) at iteration 30.


©Sky Coyote 2007