Solving equation to find center point of circle from 3 points

Question

I'm looking for a high precision solution to find the center point of a circle from 3 data points on a canvas (x,y). I found this example in the attached screenshot above, now I'm using the Math.NET package to solve the equation and I'm comparing the results against this online tool: https://planetcalc.com/8116/.

However, when I calculate the radius its completely off and often a negative number???

using MathNet.Numerics.LinearAlgebra.Double.Solvers;
using MathNet.Numerics.LinearAlgebra.Double;
using System;

namespace ConsoleAppTestBed
{
  class Program
  {
    static void Main(string[] args)
    {
        var dataPoints = new double[,]
        {

            { 5, 80 },
            { 20, 100 },
            { 40, 140 }
        };


        var fitter = new CircleFitter();
        var result = fitter.Fit(dataPoints);

        var x = -result[0];
        var y = -result[1];
        var c = result[2];

        Console.WriteLine("Center Point:");
        Console.WriteLine(x);
        Console.WriteLine(y);
        Console.WriteLine(c);

        //// (x^2 + y^2 - c^2)
        var radius = Math.Pow(x, 2) + Math.Pow(y, 2) - Math.Pow(c, 2);
        //// sqrt((x^2 + y^2 - c^2))
        radius =  Math.Sqrt(radius);
        Console.WriteLine("Radius:");
        Console.WriteLine(radius);
        Console.ReadLine();
    }

    public class CircleFitter
    {
        public double[] Fit(double[,] v)
        {
            var xy1 = new double[] { v[0,0], v[0,1] };
            var xy2= new double[] { v[1, 0], v[1, 1] };
            var xy3 = new double[] { v[2, 0], v[2, 1] };

            // Create Left Side Matrix of Equation
            var a = CreateLeftSide_(xy1);
            var b = CreateLeftSide_(xy2);
            var c = CreateLeftSide_(xy3);
            var matrixA = DenseMatrix.OfArray(new[,] 
            {
                { a[0], a[1], a[2] },
                { b[0], b[1], b[2] },
                { c[0], c[1], c[2] }
            });

            // Create Right Side Vector of Equation
            var d = CreateRightSide_(xy1);
            var e = CreateRightSide_(xy2);
            var f = CreateRightSide_(xy3);
            double[] vector = { d, e, f };
            var vectorB =  Vector<double>.Build.Dense(vector);

            // Solve Equation
            var r = matrixA.Solve(vectorB);
            var result = r.ToArray();

            return result;
        }

        //2x, 2y, 1
        public double[] CreateLeftSide_(double[] d) 
        {
            return new double[] { (2 * d[0]), (2 * d[1]) , 1};
        
        }

        // -(x^2 + y^2)
        public double CreateRightSide_(double[] d) 
        { 
            return -(Math.Pow(d[0], 2) + Math.Pow(d[1], 2));

        }
    }
  }

}

Any ideas?

Thanks in advance.

Answer 1

The solution to your problem is here: The NumberDecimalDigits property

Code:

using System;
using System.Globalization;
namespace ConsoleApp1
{
class Program
{
    static void Main()
    {
        double x1 = 1, y1 = 1;
        double x2 = 2, y2 = 4;
        double x3 = 5, y3 = -3;
        findCircle(x1, y1, x2, y2, x3, y3);
        Console.ReadKey();
    }
    static void findCircle(double x1, double y1,
                           double x2, double y2,
                           double x3, double y3)
    {
        NumberFormatInfo setPrecision = new NumberFormatInfo();
        setPrecision.NumberDecimalDigits = 3; // 3 digits after the double point

        double x12 = x1 - x2;
        double x13 = x1 - x3;

        double y12 = y1 - y2;
        double y13 = y1 - y3;

        double y31 = y3 - y1;
        double y21 = y2 - y1;

        double x31 = x3 - x1;
        double x21 = x2 - x1;

        double sx13 = (double)(Math.Pow(x1, 2) -
                        Math.Pow(x3, 2));

        double sy13 = (double)(Math.Pow(y1, 2) -
                        Math.Pow(y3, 2));

        double sx21 = (double)(Math.Pow(x2, 2) -
                        Math.Pow(x1, 2));

        double sy21 = (double)(Math.Pow(y2, 2) -
                        Math.Pow(y1, 2));

        double f = ((sx13) * (x12)
                + (sy13) * (x12)
                + (sx21) * (x13)
                + (sy21) * (x13))
                / (2 * ((y31) * (x12) - (y21) * (x13)));
        double g = ((sx13) * (y12)
                + (sy13) * (y12)
                + (sx21) * (y13)
                + (sy21) * (y13))
                / (2 * ((x31) * (y12) - (x21) * (y13)));

        double c = -(double)Math.Pow(x1, 2) - (double)Math.Pow(y1, 2) -
                                    2 * g * x1 - 2 * f * y1;
        double h = -g;
        double k = -f;
        double sqr_of_r = h * h + k * k - c;

        // r is the radius
        double r = Math.Round(Math.Sqrt(sqr_of_r), 5);

        Console.WriteLine("Center of a circle: x = " + h.ToString("N", setPrecision) + 
        ", y = " + k.ToString("N", setPrecision));
        Console.WriteLine("Radius: " + r.ToString("N", setPrecision));
        }
        }
        }