起因
网友问：“如何通过一系列的点计算出圆？”
俺的回答：“可以问AI”
先上结果
原始数据
计算的结果
传入32万点，计算耗时3毫秒
心路历程
其实这是一个常见的题目，直接搜代码应该也可以搜到。但是俺支持使用AI进行了创作。
俺的第一个提示词是：
你是一个有着良好开发习惯的C#程序员。请帮我写一个类，用于基于一系列的点，拟合出一个圆，入参是List<Point>,返回圆心和坐标。并给出算法的说明。
AI给俺了一个类 ，没啥错误直接编译。AI还给出了推理过程。
public class CircleFitter
{
/// <summary>
/// 拟合圆的结果
/// </summary>
public class FitResult
{
public PointF Center { get; set; }
public double Radius { get; set; }
public double Error { get; set; } // 平均误差
}
/// <summary>
/// 使用最小二乘法拟合圆
/// </summary>
/// <param name="points">点集合</param>
/// <returns>拟合结果，包含圆心、半径和误差</returns>
public static FitResult FitCircle(List<Point> points)
{
if (points == null || points.Count < 3)
{
throw new ArgumentException("至少需要3个点来拟合圆");
}
double sumX = 0, sumY = 0;
double sumX2 = 0, sumY2 = 0, sumXY = 0;
double sumX3 = 0, sumY3 = 0, sumXY2 = 0, sumX2Y = 0;
int n = points.Count;
foreach (var p in points)
{
double x = p.X;
double y = p.Y;
sumX += x;
sumY += y;
sumX2 += x * x;
sumY2 += y * y;
sumXY += x * y;
sumX3 += x * x * x;
sumY3 += y * y * y;
sumXY2 += x * y * y;
sumX2Y += x * x * y;
}
double A = n * sumX2 - sumX * sumX;
double B = n * sumXY - sumX * sumY;
double C = n * sumY2 - sumY * sumY;
double D = 0.5 * (n * sumXY2 - sumX * sumY2 + n * sumX3 - sumX * sumX2);
double E = 0.5 * (n * sumX2Y - sumY * sumX2 + n * sumY3 - sumY * sumY2);
// 解线性方程组
// A*a + B*b = D
// B*a + C*b = E
double det = A * C - B * B;
if (Math.Abs(det) < 1e-10)
{
throw new InvalidOperationException("点共线或太接近共线，无法拟合圆");
}
double a = (D * C - E * B) / det;
double b = (A * E - B * D) / det;
// 计算半径
double radius = 0;
double error = 0;
foreach (var p in points)
{
double dx = p.X - a;
double dy = p.Y - b;
radius += Math.Sqrt(dx * dx + dy * dy);
}
radius /= n;
计算平均误差
//foreach (var p in points)
//{
// double dx = p.X - a;
// double dy = p.Y - b;
// double actualRadius = Math.Sqrt(dx * dx + dy * dy);
// error += Math.Abs(actualRadius - radius);
//}
//error /= n;
return new FitResult
{
Center = new PointF((float)a, (float)b),
Radius = radius,
Error = error
};
}
}测试
生成测试数据
基于圆的轨迹随机生成一系列的点。32万个点 ，不耗什么时间。
private unsafe void makeTestData(Size s)
{
bmp = new Bitmap(s.Width, s.Height, PixelFormat.Format24bppRgb);
Graphics g = Graphics.FromImage(bmp);
g.Clear(Color.White);
g.Dispose();
int d = Math.Min(s.Width, s.Height);
int r = d / 3;
List<Point> list= CircleGenerator.GenerateCircle(new Point(s.Width/2, s.Height/2), r);
Random rand = new Random(DateTime.Now.Millisecond);
bool[] flag = new bool[s.Width * s.Height];
pts = new List<Point>();
int c,x,y,idx,rd_max;
rd_max = r * 9 / 10;
foreach (Point pt in list)
{
for (int rd = 3; rd < rd_max; rd+=10)
{
c = rand.Next(rd_max - rd);
for (int i = 0; i < c; i++)
{
x = pt.X + rand.Next(-rd,rd);
y = pt.Y + rand.Next(-rd,rd);
idx = s.Width * y + x;
if (idx < 0)
continue;
if (idx < flag.Length)
{
if (flag[idx])
continue;
flag[idx] = true;
pts.Add(new Point(x, y));
}
}
}
}
BitmapData bmpdata = bmp.LockBits(new Rectangle(0, 0, bmp.Width, bmp.Height), ImageLockMode.ReadWrite, PixelFormat.Format24bppRgb);
int len = bmpdata.Stride;
foreach (Point pt in pts)
{
IntPtr pline = bmpdata.Scan0 + len * pt.Y;
byte* p = (byte*)pline;
p[pt.X * 3] = 0;
p[pt.X * 3+1] = 0;
p[pt.X * 3+2] = 0;
}
bmp.UnlockBits(bmpdata);
}
计算
CircleFitter.FitResult r = CircleFitter.FitCircle(pts); 计算出圆 。
var stopwatch = Stopwatch.StartNew();
CircleFitter.FitResult r = CircleFitter.FitCircle(pts);
stopwatch.Stop();
Graphics g = Graphics.FromImage(bmp);
g.SmoothingMode = System.Drawing.Drawing2D.SmoothingMode.AntiAlias;
Pen pen = new Pen(Color.Red, 7);
g.DrawEllipse(pen, new RectangleF((float)(r.Center.X - r.Radius), (float)(r.Center.Y - r.Radius), (float)(r.Radius * 2), (float)(r.Radius * 2)));
g.Dispose();
pictureBox1.Image = bmp;
label1.Text = "基于" + pts.Count.ToString() + "个点计算圆，耗时" + stopwatch.ElapsedMilliseconds + "毫秒";
全部代码
CircleFitter
/*
提示词：
你是一个有着良好开发习惯的C#程序员。请帮我写一个类，用于基于一系列的点，拟合出一个圆，入参是List<Point>,返回圆心和坐标。并给出算法的说明。
*/
public class CircleFitter
{
/// <summary>
/// 拟合圆的结果
/// </summary>
public class FitResult
{
public PointF Center { get; set; }
public double Radius { get; set; }
public double Error { get; set; } // 平均误差
}
/// <summary>
/// 使用最小二乘法拟合圆
/// </summary>
/// <param name="points">点集合</param>
/// <returns>拟合结果，包含圆心、半径和误差</returns>
public static FitResult FitCircle(List<Point> points)
{
if (points == null || points.Count < 3)
{
throw new ArgumentException("至少需要3个点来拟合圆");
}
double sumX = 0, sumY = 0;
double sumX2 = 0, sumY2 = 0, sumXY = 0;
double sumX3 = 0, sumY3 = 0, sumXY2 = 0, sumX2Y = 0;
int n = points.Count;
foreach (var p in points)
{
double x = p.X;
double y = p.Y;
sumX += x;
sumY += y;
sumX2 += x * x;
sumY2 += y * y;
sumXY += x * y;
sumX3 += x * x * x;
sumY3 += y * y * y;
sumXY2 += x * y * y;
sumX2Y += x * x * y;
}
double A = n * sumX2 - sumX * sumX;
double B = n * sumXY - sumX * sumY;
double C = n * sumY2 - sumY * sumY;
double D = 0.5 * (n * sumXY2 - sumX * sumY2 + n * sumX3 - sumX * sumX2);
double E = 0.5 * (n * sumX2Y - sumY * sumX2 + n * sumY3 - sumY * sumY2);
// 解线性方程组
// A*a + B*b = D
// B*a + C*b = E
double det = A * C - B * B;
if (Math.Abs(det) < 1e-10)
{
throw new InvalidOperationException("点共线或太接近共线，无法拟合圆");
}
double a = (D * C - E * B) / det;
double b = (A * E - B * D) / det;
// 计算半径
double radius = 0;
double error = 0;
foreach (var p in points)
{
double dx = p.X - a;
double dy = p.Y - b;
radius += Math.Sqrt(dx * dx + dy * dy);
}
radius /= n;
计算平均误差
//foreach (var p in points)
//{
// double dx = p.X - a;
// double dy = p.Y - b;
// double actualRadius = Math.Sqrt(dx * dx + dy * dy);
// error += Math.Abs(actualRadius - radius);
//}
//error /= n;
return new FitResult
{
Center = new PointF((float)a, (float)b),
Radius = radius,
Error = error
};
}
}CircleGenerator
/*
提示词：
你是一个有着良好开发习惯的C#程序员。请帮我写一个类，基于Bresenham算法输出圆轨迹，入参是圆心和坐标，返回 List<Point> 坐标点。并给出算法的说明。
*/
public class CircleGenerator
{
/// <summary>
/// 使用Bresenham算法生成圆轨迹
/// </summary>
/// <param name="center">圆心坐标</param>
/// <param name="radius">半径</param>
/// <returns>圆轨迹上的所有点</returns>
public static List<Point> GenerateCircle(Point center, int radius)
{
var points = new List<Point>();
if (radius <= 0)
return points;
int x = 0;
int y = radius;
int d = 3 - 2 * radius;
// 生成1/8圆弧的点，然后通过对称性得到完整圆
while (x <= y)
{
AddMirroredPoints(points, center, x, y);
if (d < 0)
{
d = d + 4 * x + 6;
}
else
{
d = d + 4 * (x - y) + 10;
y--;
}
x++;
}
return points;
}
private static void AddMirroredPoints(List<Point> points, Point center, int x, int y)
{
// 利用圆的八对称性添加所有对称点
points.Add(new Point(center.X + x, center.Y + y));
points.Add(new Point(center.X - x, center.Y + y));
points.Add(new Point(center.X + x, center.Y - y));
points.Add(new Point(center.X - x, center.Y - y));
points.Add(new Point(center.X + y, center.Y + x));
points.Add(new Point(center.X - y, center.Y + x));
points.Add(new Point(center.X + y, center.Y - x));
points.Add(new Point(center.X - y, center.Y - x));
}
}总结
其实就写好提示词
后记
和网友又交流了一下，网友是要是计算出包含在圆内。
所以俺又问了一下AI。
提示词是：
你是一个有着良好开发习惯的C#程序员。请帮我写一个类适用于.Net framework4.8，用于基于一系列的点，计算出一个圆，要求这些点都位于圆内，并且圆的半径最小，入参是List<Point>, 坐标为整数。 返回圆心和坐标。并给出算法的说明
结果： 一万个点 77毫秒 。
makeTestData2(new Size(1200, 700));
pictureBox1.Image = bmp;
var stopwatch = Stopwatch.StartNew();
var result = MinimumEnclosingCircle.Calculate(pts);
stopwatch.Stop();
Graphics g = Graphics.FromImage(bmp);
g.SmoothingMode = System.Drawing.Drawing2D.SmoothingMode.AntiAlias;
Pen pen = new Pen(Color.Red, 3);
float x = result.Item1.X;
float y = result.Item1.Y;
float r= result.Item2;
g.DrawEllipse(pen, new RectangleF(x-r,y-r,r*2,r*2 ));
g.Dispose();
pictureBox1.Image = bmp;
label1.Text = "基于" + pts.Count.ToString() + "个点计算最小外接圆，耗时" + stopwatch.ElapsedMilliseconds + "毫秒";
public static class MinimumEnclosingCircle
{
public static Tuple<PointF, float> Calculate(List<Point> points)
{
if (points == null || points.Count == 0)
throw new ArgumentException("点集不能为空");
var workingPoints = new List<Point>(points);
var boundary = new List<Point>();
var result = WelzlAlgorithm(workingPoints, boundary);
return Tuple.Create(result.Center, result.Radius);
}
private struct CircleResult
{
public PointF Center;
public float Radius;
}
private static CircleResult WelzlAlgorithm(List<Point> points, List<Point> boundary)
{
if (points.Count == 0)
{
if (boundary.Count == 0)
return new CircleResult { Center = new PointF(0, 0), Radius = 0 };
if (boundary.Count == 1)
return new CircleResult { Center = boundary[0], Radius = 0 };
if (boundary.Count == 2)
return GetCircleFromTwoPoints(boundary[0], boundary[1]);
return GetCircleFromThreePoints(boundary[0], boundary[1], boundary[2]);
}
var randomIndex = new Random().Next(points.Count);
var p = points[randomIndex];
points.RemoveAt(randomIndex);
var circle = WelzlAlgorithm(points, boundary);
if (!IsPointInCircle(p, circle.Center, circle.Radius))
{
boundary.Add(p);
circle = WelzlAlgorithm(points, boundary);
boundary.RemoveAt(boundary.Count - 1);
}
points.Insert(randomIndex, p);
return circle;
}
private static bool IsPointInCircle(Point point, PointF center, float radius)
{
var dx = point.X - center.X;
var dy = point.Y - center.Y;
return dx * dx + dy * dy <= radius * radius + float.Epsilon;
}
private static CircleResult GetCircleFromTwoPoints(Point a, Point b)
{
var center = new PointF((a.X + b.X) / 2f, (a.Y + b.Y) / 2f);
var radius = (float)Math.Sqrt(Math.Pow(a.X - b.X, 2) + Math.Pow(a.Y - b.Y, 2)) / 2;
return new CircleResult { Center = center, Radius = radius };
}
private static CircleResult GetCircleFromThreePoints(Point a, Point b, Point c)
{
var d = 2f * (a.X * (b.Y - c.Y) + b.X * (c.Y - a.Y) + c.X * (a.Y - b.Y));
if (Math.Abs(d) < float.Epsilon)
return GetCircleFromTwoPoints(a, b);
var aSq = a.X * a.X + a.Y * a.Y;
var bSq = b.X * b.X + b.Y * b.Y;
var cSq = c.X * c.X + c.Y * c.Y;
var centerX = (aSq * (b.Y - c.Y) + bSq * (c.Y - a.Y) + cSq * (a.Y - b.Y)) / d;
var centerY = (aSq * (c.X - b.X) + bSq * (a.X - c.X) + cSq * (b.X - a.X)) / d;
var center = new PointF(centerX, centerY);
var radius = (float)Math.Sqrt(Math.Pow(center.X - a.X, 2) + Math.Pow(center.Y - a.Y, 2));
return new CircleResult { Center = center, Radius = radius };
}
}
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
