The best circle fitting algorithm

I need a very precise algorithm for fitting a circle to the set of data points (actually I need to determine the center). The data comes after the binarization and segmentation of the image. I tried simple center of mass and weighted center of mass algorithms and also pre-made OpenCv::fitEllipse function. I have obtained the best results from OpenCV function, but still, the accuracy is not sufficient. The results are impaired significantly, when the center is anchored in subpixel regions. The accuracy I obtained is not sufficient even when dealing with the modeled data, which is bad, as eventually, the procedure will have to deal with the data captured by comeras. Have you got any suggestions what kind of algorithm should I look for or have you got any ready solution? I would rather refrain from linking any external libs. Thanks for your help.

edited: Calibration target can be localizaed in any region of the field of view. Hereunder the best results I acieved with OpenCV procedure:

169,367 748,345  
167,557 820,788  
165,690 893,158  
164,047 965,197  
162,715 1036,729  
161,575 1108,089  
160,477 1179,552  
233,297 1015,313  
232,076 1086,965  
220,359 1229,578  
268,494 1160,275  
339,544 1162,980  
362,017 1235,669  
433,390 1238,491  
482,754 1168,299  
505,233 1241,039  
554,856 1170,664  
577,302 1243,439  
627,331 1172,795  
649,507 1245,665  
713,572 588,896  
711,995 661,853  
710,440 735,034  
708,722 808,856  
707,018 882,674  
705,377 956,169  
703,609 1029,211  
701,716 1101,950  
699,760 1174,689  
721,895 1247,620  
785,829 614,754  
784,344 687,750  
782,819 761,315  
781,292 835,225  
779,389 908,975  
777,619 982,335  
775,688 1055,275  
773,672 1128,091  
771,603 1200,724  

Editted: The numerically generated model and the real coordinates of the centres:

51,1    79,8
51,1    179,8
51,1    279,8
51,1    379,8
51,1    479,8
51,1    579,8
51,1    679,8
51,1    779,8
51,1    879,8
51,1    979,8
51,1    1079,8
51,1    1179,8
51,1    1279,8
51,1    1379,8
51,1    1479,8
151,1   79,8
151,1   179,8
151,1   279,8
151,1   379,8
151,1   479,8
151,1   579,8
151,1   679,8
151,1   779,8
151,1   879,8
151,1   979,8
151,1   1079,8
151,1   1179,8
151,1   1279,8
151,1   1379,8
151,1   1479,8
251,1   79,8
251,1   179,8
251,1   279,8
251,1   379,8
251,1   479,8
251,1   579,8
251,1   679,8
251,1   779,8
251,1   879,8
251,1   979,8
251,1   1079,8
251,1   1179,8
251,1   1279,8
251,1   1379,8
251,1   1479,8
351,1   79,8
351,1   179,8
351,1   279,8
351,1   379,8
351,1   479,8
351,1   579,8
351,1   679,8
351,1   779,8
351,1   879,8
351,1   979,8
351,1   1079,8
351,1   1179,8
351,1   1279,8
351,1   1379,8
351,1   1479,8
451,1   79,8
451,1   179,8
451,1   279,8
451,1   379,8
451,1   479,8
451,1   579,8
451,1   679,8
451,1   779,8
451,1   879,8
451,1   979,8
451,1   1079,8
451,1   1179,8
451,1   1279,8
451,1   1379,8
451,1   1479,8
551,1   79,8
551,1   179,8
551,1   279,8
551,1   379,8
551,1   479,8
551,1   579,8
551,1   679,8
551,1   779,8
551,1   879,8
551,1   979,8
551,1   1079,8
551,1   1179,8
551,1   1279,8
551,1   1379,8
551,1   1479,8
651,1   79,8
651,1   179,8
651,1   279,8
651,1   379,8
651,1   479,8
651,1   579,8
651,1   679,8
651,1   779,8
651,1   879,8
651,1   979,8
651,1   1079,8
651,1   1179,8
651,1   1279,8
651,1   1379,8
651,1   1479,8
751,1   79,8
751,1   179,8
751,1   279,8
751,1   379,8
751,1   479,8
751,1   579,8
751,1   679,8
751,1   779,8
751,1   879,8
751,1   979,8
751,1   1079,8
751,1   1179,8
751,1   1279,8
751,1   1379,8
751,1   1479,8
851,1   79,8
851,1   179,8
851,1   279,8
851,1   379,8
851,1   479,8
851,1   579,8
851,1   679,8
851,1   779,8
851,1   879,8
851,1   979,8
851,1   1079,8
851,1   1179,8
851,1   1279,8
851,1   1379,8
851,1   1479,8
951,1   79,8
951,1   179,8
951,1   279,8
951,1   379,8
951,1   479,8
951,1   579,8
951,1   679,8
951,1   779,8
951,1   879,8
951,1   979,8
951,1   1079,8
951,1   1179,8
951,1   1279,8
951,1   1379,8
951,1   1479,8
1051,1  79,8
1051,1  179,8
1051,1  279,8
1051,1  379,8
1051,1  479,8
1051,1  579,8
1051,1  679,8
1051,1  779,8
1051,1  879,8
1051,1  979,8
1051,1  1079,8
1051,1  1179,8
1051,1  1279,8
1051,1  1379,8
1051,1  1479,8
1151,1  79,8
1151,1  179,8
1151,1  279,8
1151,1  379,8
1开发者_高级运维151,1  479,8
1151,1  579,8
1151,1  679,8
1151,1  779,8
1151,1  879,8
1151,1  979,8
1151,1  1079,8
1151,1  1179,8
1151,1  1279,8
1151,1  1379,8
1151,1  1479,8

---

An Algorithm using Image Transformations and Clustering

---

I made up a small Algorithm using Image Transformations and some Statistic to detect your circles. Let's see if it is up to your error expectation.
Any good image and statistics library will do, I implemented it using Mathematica.

Run as follows:

1. Import your image and run a Bottom Hat Transform

We start trying to isolate the circles. The Bottom Hat Transform with a Box Matrix kernel helps. Almost any image library comes with the algorithm already implemented.

a = Import@"http://i.stack.imgur.com/hiSjj.png";   
b = BottomHatTransform[Binarize@a, BoxMatrix[30]]  

The result is

2. Run a Hit Miss Transform to isolate the circles

The Hit Miss Transform excels in finding well defined geometrical objects. It is also easy to program and is almost always present in image libraries.

c = Binarize@HitMissTransform[b, DiskMatrix[20]]

The result is:

And our circles are already isolated and reduced to their central core.

3. Get just the white pixels from image

This is an implementation-dependent step, so I'll not comment on this one.

ttflat = Flatten[Table[{i, j, ImageData[c][[i, j]]}, {i, 1232}, {j, 1624}], 1];  
ttfilter = Select[ttflat, #[[3]] == 1 &];  

Let's see how many pixels are left

Dimensions@ttfilter  
{3684, 3}   

So 3684 pixels left, almost 82 per circle. Enough to do some statistics.

3. Use Cluster Analysis to pick each circle

Cluster Analysis may be an overkill here, but as I have it already implemented, is easier to use it than program something new :). You may do your own or use a stats library.

ttc = FindClusters[ttfilter, 45, Method -> {"Agglomerate", "Linkage" -> "Complete"}];

With our clusters already found, let's find the mean for x and y in each cluster. Those are the centers of the circles:

means = N[Mean /@ ttc, 5]  

The result is a list of 45 coordinates like:

{{161.67, 1180.1}, {162.75, 1108.9}, 
 {164.11, 1037.6}, {165.47, 966.19} .....  

We are almost done.

Let's check our result. We superimposed both images, drawing crosses and circles around the detected centers.

Click to enlarge, so you may get an idea of the errors involved.

HTH!

Edit

I compared the results from your table with my results.

Assumming the circles are in straight lines, I used Least Squares Fit to trace a line and calculated the residuals.

Form the graph below, you may see that "M"y line fit better than "Y"ours. But that is assuming the circles aligned ...

Edit 2

These are the calculated coordinates for the first 45 circles in your second image. I have a systematic offset of 1 pixel. Probably due to some image manipulation I did, but is easy to correct :) ... just subtracted one pixel on X and Y ...

{{51.135, 79.692}, {51.135, 179.69}, {51.135, 279.69},{51.135, 379.69}, {51.135, 479.69},
 {51.135, 579.69}, {51.135, 679.69}, {51.135, 779.69},{51.135, 879.69}, {51.135, 979.69}, 
 {51.135, 1079.7}, {51.135, 1179.7}, {51.135, 1279.7},{51.135, 1379.7}, {51.135, 1479.7}, 
 {151.13, 79.692}, {151.13, 179.69}, {151.13, 279.69},{151.13, 379.69}, {151.13, 479.69},
 {151.13, 579.69}, {151.13, 679.69}, {151.13, 779.69},{151.13, 879.69}, {151.13, 979.69}, 
 {151.13, 1079.7}, {151.13, 1179.7}, {151.13, 1279.7},{151.13, 1379.7}, {151.13, 1479.7}, 
 {251.13, 79.692}, {251.13, 179.69}, {251.13, 279.69},{251.13, 379.69}, {251.13, 479.69}, 
 {251.13, 579.69}, {251.13, 679.69}, {251.13, 779.69},{251.13, 879.69}, {251.13, 979.69}, 
 {251.13, 1079.7}, {251.13, 1179.7}, {251.13, 1279.7},{251.13, 1379.7}, {251.13, 1479.7}}

And this is the image:

---

"Best" depends on the kind of noise in the input data. The problem is trivial if there is no noise in the source data points: Just pick 3 points and calculate the circle.

If you expect normal distributed, independent translations of each data point, then a least mean squares algorithm should be optimal. The data points should fit the equation:

(x - xm)^2 + (y - ym)^2 = r^2

where xm, ym, r are unknown, so:

x^2 - 2*x*xm + xm^2 + y^2 - 2*y*ym + ym^2 = r^2

substitute c for r^2-xm^2-ym^2 and you have an overdetermined system of linear equations:

2*x*xm + 2*y*ym = c - x^2 - y^2

Any good linear algebra library (e.g. IPP) can solve that for you.

If you expect outliers in the data, I would suggest using a RANSAC-strategy to find the set of non-outlier points, then use the algorithm above to find the exact center for that set.

---

To find circles' centers with subpixel accuracy if the radius of the circle is known (and constant), I use this approach:

1. Take an image with a single circle (refferred to as marker below). Its radius should be the same as the radius of the circles you want to find.
2. Detect edges (magnitude of the Sobel gradient, then some threshold to remove low-intensity edges), both in the test image and in the flipped marker image. I do not apply edge thinning at this stage and do not identify exact points on the edge.
3. Cross-correlate edges of the test image with the edges of the flipped marker. You get some peaks where the centers are located.
4. Find peaks' centers with subpixel accuracy. Center of mass or fitting a 2D Gaussian bell may work well.
5. Add shifts which correspond to known position of the center of the circle in the marker.

Otherwise, if points on the circle are known with sufficient accuracy, least mean squares fitting should solve the problem of finding the center (see @nikie's answer).

---

Take a look at Hough transform. It can be used to detect circles.

---

OpenCV 2.4.6.0 has findCirclesGrid function to find the centers of circles in a grid.