c# MathNet.Numerics 圆拟合使用案例
下载地址
https://www.nuget.org/packages/MathNet.Numerics/4.15.0#supportedframeworks-body-tab
原理
从标准圆方程(x-c1)2+(y-c2)2=r2中进行方程变换得到2xc1+2yc2+(r2−c12−c22)=x2+y2,其中,我们c3替换常量值r2−c12−c22,即:r2−c12−c22=c3。由此,我们得到2xc1+2yc2+c3=x2+y2,将点集带入,方程就只剩三个未知数`c1,c2 和 c3。
案例
using MathNet.Numerics;
using MathNet.Numerics.LinearAlgebra.Double;
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;namespace MathNetTest
{class Program{static void Main(string[] args){FitCircle();}private static void FitCircle(){var count = 5000;var step = 2 * Math.PI / 100;var rd = new Random();//参照圆var x0 = 204.1;var y0 = 213.1;var r0 = 98.4;//噪音绝对差var diff = (int)(r0 * 0.1);var x = new double[count];var y = new double[count];//输出点集for (int i = 0; i < count; i++){//circlex[i] = x0 + Math.Cos(i * step) * r0;y[i] = y0 + Math.Sin(i * step) * r0;//noisex[i] += Math.Cos(rd.Next() % 2 * Math.PI) * rd.Next(diff);y[i] += Math.Cos(rd.Next() % 2 * Math.PI) * rd.Next(diff);}var d = LinearAlgebraFit(x, y);}public static Tuple<double,double,double> LinearAlgebraFit(double[] X, double[] Y){if (X.Length < 3){return null;}var count = X.Length;var a = new double[count, 3];var c = new double[count, 1];for (int i = 0; i < count; i++){//matrixa[i, 0] = 2 * X[i];a[i, 1] = 2 * Y[i];a[i, 2] = 1;c[i, 0] = X[i] * X[i] + Y[i] * Y[i];}var A = DenseMatrix.OfArray(a);var C = DenseMatrix.OfArray(c);//A*B=Cvar B = A.Solve(C);double c1 = B.At(0, 0),c2 = B.At(1, 0),r = Math.Sqrt(B.At(2, 0) + c1 * c1 + c2 * c2);Tuple<double, double, double> XYR = new Tuple<double, double, double>(c1, c2, r);return XYR;}
}