关于拟合圆的方法有很多，主要演示以下三种：
首先是准备工作，模拟一组随便数组：
模拟点集，由于现实中的数据采集存在着精度、数据记录等众多不确定因素的影像。模拟点集中也将加入一定程度的噪音。以下代码中 x 与 y
中存储着我们的点集数据：
const int count = 100;
double[] x = new double[count + 1];
double[] y = new double[count + 1];
double step = 2 * Math.PI / count;
Random rd = new Random();
//参照圆
double centerX = 204.1;
double centerY = 213.1;
double radius = 98.4;
//噪音绝对差
var diff = (int)(radius * 0.1);
//输出点集
for (int i = 0; i < count; i++)
{
//circle
x[i] = centerX + Math.Cos(i * step) * radius;
y[i] = centerY + Math.Sin(i * step) * radius;
//noise
x[i] += Math.Cos(rd.Next() % 2 * Math.PI) * rd.Next(diff);
y[i] += Math.Cos(rd.Next() % 2 * Math.PI) * rd.Next(diff);
}
x[count] = x[0];//闭合圆
y[count] = y[0];
以下是一个基于最小二乘法和线性代数的圆拟合类的示例。它通过输入一组点的坐标，返回这些点的最佳拟合圆（在平面上）。
//多个点拟合圆并给出圆心坐标。
public class FitCircle
{
/// <summary>
/// 拟合圆
/// </summary>
private double centerX;
private double centerY;
private double radius;
public double CenterX { get { return centerX; } }
public double CenterY { get { return centerY; } }
public double Radius { get { return radius; } }
/// <summary>
/// 最小二乘法拟合圆
/// </summary>
public void LeastSquaresFit(double[] X, double[] Y)
{
if (X.Length != Y.Length || X.Length < 3)
{
return;
}
double sum_x = 0.0f, sum_y = 0.0f;
double sum_x2 = 0.0f, sum_y2 = 0.0f;
double sum_x3 = 0.0f, sum_y3 = 0.0f;
double sum_xy = 0.0f, sum_x1y2 = 0.0f, sum_x2y1 = 0.0f;
int N = X.Length;
double x, y, x2, y2;
for (int i = 0; i < N; i++)
{
x = X[i];
y = Y[i];
x2 = x * x;
y2 = y * y;
sum_x += x;
sum_y += y;
sum_x2 += x2;
sum_y2 += y2;
sum_x3 += x2 * x;
sum_y3 += y2 * y;
sum_xy += x * y;
sum_x1y2 += x * y2;
sum_x2y1 += x2 * y;
}
double C, D, E, G, H;
double a, b, c;
C = N * sum_x2 - sum_x * sum_x;
D = N * sum_xy - sum_x * sum_y;
E = N * sum_x3 + N * sum_x1y2 - (sum_x2 + sum_y2) * sum_x;
G = N * sum_y2 - sum_y * sum_y;
H = N * sum_x2y1 + N * sum_y3 - (sum_x2 + sum_y2) * sum_y;
a = (H * D - E * G) / (C * G - D * D);
b = (H * C - E * D) / (D * D - G * C);
c = -(a * sum_x + b * sum_y + sum_x2 + sum_y2) / N;
centerX = a / (-2);
centerY = b / (-2);
radius = Math.Sqrt(a * a + b * b - 4 * c) / 2;
}
/// <summary>
/// 线性代数拟合圆
/// </summary>
public void FLinearAlgebraFit(double[] X, double[] Y)
{
if (X.Length != Y.Length || X.Length < 3)
{
return;
}
var count = X.Length;
var a = new double[count, 3];
var c = new double[count, 1];
for (int i = 0; i < count; i++)
{
//matrix
a[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];
}
//using MathNet.Numerics.LinearAlgebra.Double;
var A = DenseMatrix.OfArray(a);
var C = DenseMatrix.OfArray(c);
//A*B=C
var 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);
centerX = c1;
centerY = c2;
radius = r;
}
/// <summary>
/// 拟合圆程序
/// </summary>
/// <param name="pPointList">要拟合点集</param>
/// <returns>返回圆对象</returns>
public void circle_fitting_2D(double[] X, double[] Y)
{
if (X.Length != Y.Length || X.Length < 3)
{
Console.WriteLine("最少需要三个点进行拟合");
return;
}
double X1 = 0;
double Y1 = 0;
double X2 = 0;
double Y2 = 0;
double X3 = 0;
double Y3 = 0;
double X1Y1 = 0;
double X1Y2 = 0;
double X2Y1 = 0;
for (int i = 0; i <Math.Min(X.Length,Y.Length); i++)
{
X1 = X1 + X[i];
Y1 = Y1 + Y[i];
X2 = X2 + X[i] * X[i];
Y2 = Y2 + Y[i] * Y[i];
X3 = X3 + X[i] * X[i] * X[i];
Y3 = Y3 + Y[i] * Y[i] * Y[i];
X1Y1 = X1Y1 + X[i] * Y[i];
X1Y2 = X1Y2 + X[i] * Y[i] * Y[i];
X2Y1 = X2Y1 + X[i] * X[i] * Y[i];
}
double C, D, E, G, H, N;
double a, b, c;
N = X.Length;
C = N * X2 - X1 * X1;
D = N * X1Y1 - X1 * Y1;
E = N * X3 + N * X1Y2 - (X2 + Y2) * X1;
G = N * Y2 - Y1 * Y1;
H = N * X2Y1 + N * Y3 - (X2 + Y2) * Y1;
a = (H * D - E * G) / (C * G - D * D);
b = (H * C - E * D) / (D * D - G * C);
c = -(a * X1 + b * Y1 + X2 + Y2) / N;
centerX = a / (-2);
centerY = b / (-2);
radius = Math.Sqrt(a * a + b * b - 4 * c) / 2;
}
}
使用示例：
FitCircle fitCircle=new FitCircle();
fitCircle.LeastSquaresFit(x, y);
string str1 = string.Format("Center: ({0}, {1})",
fitCircle.CenterX.ToString("F2"), fitCircle.CenterY.ToString("F2"));
string str2 = string.Format("Radius: {0}", fitCircle.Radius.ToString("F2"));
Console.WriteLine(str1);
Console.WriteLine(str2);
fitCircle.FLinearAlgebraFit(x, y);
string str3 = string.Format("Center: ({0}, {1})",
fitCircle.CenterX.ToString("F2"), fitCircle.CenterY.ToString("F2"));
string str4 = string.Format("Radius: {0}", fitCircle.Radius.ToString("F2"));
Console.WriteLine(str3);
Console.WriteLine(str4);
fitCircle.circle_fitting_2D(x, y);
string str5 = string.Format("Center: ({0}, {1})",
fitCircle.CenterX.ToString("F2"), fitCircle.CenterY.ToString("F2"));
string str6 = string.Format("Radius: {0}", fitCircle.Radius.ToString("F2"));
Console.WriteLine(str5);
Console.WriteLine(str6);
输出：
Center: (204.28, 213.14)
Radius: 98.00
Center: (204.28, 213.14)
Radius: 98.00
Center: (204.28, 213.14)
Radius: 98.00
绘制圆点坐标：
//picGraphics是我用picturePox类中的CreateGraphics方法创建的Graphics对象
Graphics picGraphics = this.picbox.CreateGraphics();
picGraphics.Clear(SystemColors.Control);
RectangleF rectangle2 = new RectangleF((float)((fitCircle.CenterX - fitCircle.Radius)), (float)((fitCircle.CenterY - fitCircle.Radius)),
(float)(2 * fitCircle.Radius), (float)(2 * fitCircle.Radius));
picGraphics.TranslateTransform(150f, 200f);
picGraphics.DrawEllipse(new Pen(Color.GreenYellow, 2), rectangle2);
int fillSize = 4;
for (int i = 0; i < Math.Min(x.Length, y.Length); i++)
{
picGraphics.FillRectangle(new SolidBrush(Color.Red), (int)(x[i] - fillSize / 2), (int)(y[i] - fillSize / 2), fillSize, fillSize);// 即用一个像素填充方法.
}
最后绘制结果如下图所示：
在最小二乘法中，只有一个及其简单的 for 循环，很少涉及内存写。但在线性代数中，需要进行矩阵的生成DenseMatrix.OfArray，以及矩阵运算，这二者都需要内存写。再者，矩阵计算有着繁重的计算量，这些都在影响着线性代数拟合圆的效率。最终的胜利还是属于最小二乘法。
扫一扫
6174
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
