Find circle center from a segment of surface points

Question

I need to find the center of a circle based on a number of surface measurements of a cylindrical object.

I currently use a simplified algorithm (written in C#) based on three points to find the center (taken from How to determine the radius and center of a circle when only three noncollinear points are known?):

private Point CircleCenter(List<Point> points, double p1Skip, double p2Skip, double p3Skip)
{
    var p1 = points.Skip((int)(points.Count * p1Skip)).First();
    var p2 = points.Skip((int)(points.Count * p2Skip)).First();
    var p3 = points.Skip((int)(points.Count * p3Skip)).First();

    double mr = (p2.Y - p1.Y) / (p2.X - p1.X);
    double mt = (p3.Y - p2.Y) / (p3.X - p2.X);

    double centerX = (mr * mt * (p3.Y - p1.Y) + mr * (p2.X + p3.X) - mt * (p1.X + p2.X)) / (2 * (mr - mt));

    double centerY = (-1 / mr) * (centerX - ((p1.X + p2.X) / 2)) + ((p1.Y + p2.Y) / 2);

    return new Point(centerX, centerY, p1.Z);
}

The problem is that this method is very sensitive to noise. If one of the points is off, it will, of course, affect the center point. I have 640 surface points available, for each cross section, and figure it should be possible to use much more than 3 points.

I'm guessing that the existing algorithm should be possible to extend with more points, but I can't figure out how.

Answer 1

With the help of my math guru colleague I've arrived at a solution using least square matrix calculation.

The algorithm used is described here.

public static void FitCircle(IEnumerable<Point> points, out double x0, out double y0, out double r)
{
    setMklProvider();
    DenseMatrix A = DenseMatrix.Create(3, 3, (i, j) => 0);
    DenseMatrix b = DenseMatrix.Create(3, 1, (i, j) => 0);
    A[0, 0] = points.Sum(point => point.X * point.X);
    A[0, 1] = points.Sum(point => point.X * point.Y);
    A[0, 2] = points.Sum(point => point.X);
    A[1, 0] = A[0, 1];
    A[1, 1] = points.Sum(point => point.Y * point.Y);
    A[1, 2] = points.Sum(point => point.Y);
    A[2, 0] = A[0, 2];
    A[2, 1] = A[1, 2];
    A[2, 2] = points.Count();
    b[0, 0] = points.Sum(point => point.X * (point.X * point.X + point.Y * point.Y));
    b[1, 0] = points.Sum(point => point.Y * (point.X * point.X + point.Y * point.Y));
    b[2, 0] = points.Sum(point => point.X * point.X + point.Y * point.Y);
    var x = A.QR().Solve(b);
    x0 = x[0, 0] / 2;
    y0 = x[1, 0] / 2;
    r = Math.Sqrt(x[2, 0] + x0 * x0 + y0 * y0);
}

private static void setMklProvider()
{
    if (!_mklProviderSet) MathNet.Numerics.Control.LinearAlgebraProvider = new MathNet.Numerics.Algorithms.LinearAlgebra.Mkl.MklLinearAlgebraProvider();
}

This solution produces really nice repeatable results, at least for my data. The DenseMatrix is a part of the MathNet library.

EDIT

To further reduce noise, as user samgak suggested, I added an iterative reduction approach that improved the accuracy:

double x0, y0, r;
FitCircle(surfacePoints, out x0, out y0, out r);
var center = new Point(x0, y0, surfacePoints.First().Z);

int reductionIterations = 10;
var reducedSet = surfacePoints;

for (int i = 1; i < reductionIterations; i++)
{
    var orderedByDistanceToCenter = reducedSet.OrderBy(p => (p-center).GetRho()).ToList();

    reducedSet = orderedByDistanceToCenter
        .Skip((int)(orderedByDistanceToCenter.Count * (i / 10f)))
        .Take((int)(orderedByDistanceToCenter.Count - orderedByDistanceToCenter.Count * (i / 10f)*2))
        .ToList();

    // Reduced to zero, abort
    if (reducedSet.Count < 3)
        break;

    FitCircle(reducedSet, out x0, out y0, out r);
    center = new Point(x0, y0, reducedSet.First().Z);
}

public static double GetRho(this Point p) => Math.Sqrt(p.X * p.X + p.Y * p.Y);