Linear Least Square Fitting - Circle

[pal...@gmail.com]

I managed to fit a circle to a Points, although the code is using List<IntPoint>, but the same algorithm should work for System.Drawing.Point, System.Windows.Point, AForge.Point, AForge.DoublePoint, or any other point types.

The idea is simple, here are the maths:

Given a set of points:
{ (x1, y1), (x2, y2), ... , (xn, yn) }
find the circle that best-fit these points.

First we know that the equation of a circle is:
(x - x0)^2 + (y - y0)^2 = r^2
Expand it we get:
(x^2 - 2 * x0 * x + x0^2) + (y^2 - 2 * y0 * y + y0^2) = r^2
Reordering them to form "ax + by = c":
(2 * x0 * x) + (2 * y0 * y) + (r^2 - x0^2 - y0^2) = x^2 + y^2
So we get:
c1 = 2 * x0;
c2 = 2 * y0;
c3 = r^2 - x0^2 - y0^2

So the matrix is:
A c Y
============================================
[ 2 * x1, 2 * y1, 1 ] [ x1^2 + y1^2 ]
[ 2 * x2, 2 * y2, 1 ] [ c1 ] [ x2^2 + y2^2 ]
[ 2 * x3, 2 * y3, 1 ] * [ c2 ] = [ x3^2 + y3^2 ]
[ ... ] [ c3 ] [ ... ]
[ 2 * xn, 2 * yn, 1 ] [ xn^2 + yn^2 ]
now solve for c we get the coefficients of the best fit circle equation.

first multiply A and Y by AT
B = AT * M;
Z = AT * Y;

then solve for c:
c = Solve(B, Z);

public Circle FitCircle(List<IntPoint> points)
{
double[,] A = new double[points.Count, 3];
double[,] Y = new double[points.Count, 1];

// setup the matrices
for (int i = 0;i < points.Count;++i)
{
// we solve for [ 2*c1 2*c2 c3 ] here,
// avoid doing 2*xn / 2*yn in the loop.
A[i, 0] = points[i].X;
A[i, 1] = points[i].Y;
A[i, 2] = 1;
Y[i, 0] = points[i].X * points[i].X + points[i].Y * points[i].Y;
}

// get AT * A and AT * Y
double[,] AT = A.Transpose();
double[,] B = AT.Multiply(A);
double[,] Z = AT.Multiply(Y);

// solve for c
double[,] c = Matrix.Solve(B, Z, true);

// okay now we get the circle :-D
double x = c[0, 0] * 0.5;
double y = c[1, 0] * 0.5;
double r = System.Math.Sqrt(c[2, 0] + x * x + y * y);

return new Circle(x, y, r);
}

Elegant! Isn't it?

ps. Line fitting and Ellipse fitting can also be done, but I just don't have the need to write them... XD
Maybe I'll implement them when I got more time.

[César]

Indeed, elegant :-) 

Thanks for sharing the code. Can I add it as a user code contribution, extending one of the shape fitting algorithms?

You might find a (robust) line fitting algorithm on the Accord.MachineLearning.Geometry namespace, at https://code.google.com/p/accord/source/browse/trunk/Sources/Accord.MachineLearning/Geometry/RansacLine.cs

Hope you will find it interesting!

Best regards,

Cesar

[pal...@gmail.com]

Sure, I wouldn't post it here if I would just like to keep that code secret. XD

I made a RansacCircle class to fit circle using the above algorithm, but it usually give wrong result (because of the randomness i guess...).
Like 7/10 the time it doesn't fit my circle (my datasets are pretty noisy, tho)...

This is pretty much just a clone of the RansacLine class.

=== cut here ===
using System;
using Accord.MachineLearning;
using Accord.Math;
using Accord.Math.Decompositions;
using AForge;
using AForge.Math.Geometry;
using System.Collections.Generic;
using System.Linq;
using FuseAOI.Inspector;
using FuseAOI.Inspector.Math.Geometry;

public class RansacCircle
{
private RANSAC<Circle> ransac;
private int[] inliers;

private Point[] points;
private double[] d2;

public RANSAC<Circle> Ransac { get { return this.ransac; } }
public int[] Inliers { get { return inliers; } }

public RansacCircle(double threshold, double probability)
{
ransac = new RANSAC<Circle>(3, threshold, probability);
ransac.Fitting = definecircle;
ransac.Distances = distance;
ransac.Degenerate = degenerate;
}

public Circle Estimate(IEnumerable<IntPoint> points)
{
return Estimate(points.Select(p => new Point(p.X, p.Y)).ToArray());
}

public Circle Estimate(IEnumerable<Point> points)
{
return Estimate(points.ToArray());
}

public Circle Estimate(IntPoint[] points)
{
return Estimate(points.Apply(p => new Point(p.X, p.Y)));
}

public Circle Estimate(Point[] points)
{
if (points.Length < 3)
throw new ArgumentException("At least three points are required to fit a circle");

this.d2 = new double[points.Length];
this.points = points;

ransac.Compute(points.Length, out inliers);
Circle circle = fitting(points.Submatrix(inliers));

return circle;
}

private Circle definecircle(int[] x)
{
System.Diagnostics.Debug.Assert(x.Length == 3);
return new Circle(points[x[0]], points[x[1]], points[x[2]]);
}

private int[] distance(Circle c, double t)
{
for (int i = 0; i < points.Length; i++)
d2[i] = c.DistanceToPoint(points[i]);

return Matrix.Find(d2, z => z < t);
}

private bool degenerate(int[] indices)
{
System.Diagnostics.Debug.Assert(indices.Length == 3);

Point p1 = points[indices[0]];
Point p2 = points[indices[1]];
Point p3 = points[indices[2]];

return p1 == p2 || p2 == p3 || p1 == p3;
}

static Circle fitting(Point[] points)
{
if (points.Length == 3)
return new Circle(points[0], points[1], points[2]);

return ShapeFitter.FitCircle(points);
}
}

public class Circle
{
public double Area { get { return Radius * Radius * System.Math.PI; } }
public double Radius { get; set; }
public AForge.Point Origin { get; set; }

public Circle()
{
Origin = new AForge.Point(0, 0);
Radius = 0;
}

public Circle(float x, float y, double radius)
{
Origin = new AForge.Point(x, y);
Radius = radius;
}

public Circle(AForge.Point origin, double radius)
{
Origin = origin;
Radius = radius;
}

public Circle(AForge.Point p1, AForge.Point p2, AForge.Point p3)
{
// ya = ma * (x - x1) + y1
// yb = mb * (x - x2) + y2
//
// ma = (y2 - y1) / (x2 - x1)
// mb = (y3 - y2) / (x3 - x2)
double ma = (p2.Y - p1.Y) / (p2.X - p1.X);
double mb = (p3.Y - p2.Y) / (p3.X - p2.X);

// (ma * mb * (y1 - y3) + mb * (x1 + x2) - ma * (x2 + x3)
// x = ----------------------------------------------------------
// 2 * (mb - ma)
double x = (ma * mb * (p1.Y - p3.Y) + mb * (p1.X + p2.Y) - ma * (p2.X + p3.X)) / (2 * (mb - ma));
double y = ma * (x - p1.X) + p1.Y;

Origin = new AForge.Point((float)x, (float)y);
Radius = Origin.DistanceTo(p1);
}

public double DistanceToPoint(AForge.Point point)
{
return System.Math.Abs(Origin.DistanceTo(point) - Radius);
}
}
=== cut here ===

[pal...@gmail.com]

My data set looks a lot like this:
https://docs.google.com/drawings/d/1K4vYGJxPbb4wgUaQabPE_WAPjwzj_YQCdBnD4_lF0HA/edit?usp=sharing

the blue line is the circle I want to fit, but there are lots of noise (black lines) around it.

my RansacCircle hardly fits the circle, tho...
like 70% of the time it gives either super huge circle whos origin is even outside the image, or super small circle which only covers the noises in the middle.

So I'm currently using other circle fitting algorithm (inspired by EDCircles, but a lot simpler in segmentation).

[César]

Thanks!

Indeed, I am afraid there could be something odd with the Ransac<Shape> classes. But by seeing your example, it could be because it would likely work better to cases where there could be multiple points around a not necessarily continuous/smooth ring. By the way, feel free to continue contributing whenever you wish!

Best regards,

Cesar

[anh...@gmail.com]

There is a comparison of different fitting algorithms:
http://people.cas.uab.edu/~mosya/cl/CPPcircle.html