找到一个具有最大点数的圆 ((x,y,r))；给定二维平面中的一组点（x，y）

Posted 2023-02-22

【中文标题】找到一个具有最大点数的圆 ((x,y,r))；给定二维平面中的一组点（x，y）

【英文标题】：Find a Circle ((x,y,r)) that has maximum number of points 'on' it; given a set of points(x,y) in a 2D plane

【发布时间】：2019-05-12 04:16:07

【问题描述】：

给定平面上的 N 个点（形式为 (x, y)），找到其上点数最多的圆？ 附： ：该点应位于圆的圆周上。 解决这个问题的最有效算法是什么？它是如何工作的？您将使用哪些数据结构来解决此问题。这是在一次 FANG 编码采访中提出的。

【问题讨论】：

在当前公式中：选择任意x,y 并将r 设置为足够大以覆盖所有点。

@MBo，我理解的问题不同：找到圆线上点最多的圆。

@Aldert 也许你是对的，问题是关于圆周上的点。希望作者把事情说的更清楚。

如果问题的意思是“在曲线本身上”，那么这里有一个 O(n^4) 解决方案的提示：有多少个圆经过 3 个非共线点？顺便说一句，任何解决方案都面临一些精度问题。

@MBo 圆就是曲线。圆内的点集是一个圆盘

【参考方案1】：

作为起点，简单的 O(N³) 解决方案是找到对应于每个唯一三元组的圆，同时计算出现的次数您找到的每个圈子。

如果一个圆圈上有 N 个点，那么你会找到它 N-choose-3 次，所以你最常找到的那个圆圈是最多分。

在任何实际实施中都有复杂性，但它们是不同的复杂性，具体取决于您的观点是如何表示的以及您想要准确还是近似的答案。

【讨论】：

即使我也有类似的方法，找到所有可能的圆，取三个点并计算每个圆的出现次数 (x,y,R)。 O(N^3) 此时，每个圆的出现实际上与该圆上的点数无关。所以要找到圆上的实际点数；我需要再次遍历所有点（O(N)）并验证哪些点在这个特定的圆圈上。虽然这里也可能存在一些并发症。我想知道这是否可以以较低的复杂性来完成。

【参考方案2】：

案例 1：霍夫变换

在计算机视觉问题解决中，我们经常在边缘信息中搜索圆圈。这个问题的特点是有很多数据点，可能来自许多不同的圆圈，存在大量噪音。

解决此问题的常用方法是 霍夫变换 https://en.wikipedia.org/wiki/Circle_Hough_Transform。基本思想是对可以通过每个点 (x, y) 的圆的证据求和。

我们创建了一个名为 Hough [a, b, r] 的整数数组，它参数化了所有可能通过您的点 (x,y) 的圆。这相当于在以(x,y)为中心的r=1平面上画一个半径为1的圆；以 (x,y) 等为中心的 r=2 平面中半径为 2 的圆。

每次通过 [a, b, r] 中的一个点绘制一个圆时，我们将相应的值加 1。有些点积累了很多的证据。这些点对应于感兴趣的圆圈。

来自 cis.rit.edu 的图片说明了其中一个 r 平面中发生的情况。

对每个点 (x,y) 执行此操作将生成指向 [a,b,r] 中与您寻找的圆圈相对应的每个点的证据。因此，只需扫描此数组即可找到具有最大证据的点。那是你的圈子。

霍夫变换示例

知道圆的半径可以将这个问题从 O(n^3) 问题减少到 O(n^2) 问题，因为只需要构建和扫描一个平面。我在对数空间中绘制半径以给出（不太准确）O（n^2 log n）算法也有很好的结果。

案例 2：圆形拟合

如果已知点位于单个圆的边界附近，和/或如果点不是很多和/或我们非常确定噪声非常小，则霍夫变换是一个糟糕的解决方案因为它是计算密集型的，内存占用很大，并且由于累加器数组的光栅特性可能不是很准确。

在这种情况下，我们可能希望通过类似于使用线性回归的线拟合技术来拟合圆。圆拟合技术的讨论可以在https://pdfs.semanticscholar.org/faac/44067f04abf10af7dd583fca0c35c5937f95.pdf找到。

这个算法的一个（相当简单的）实现如下所示。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

typedef struct  
    float x;
    float y;
     point;

/*
 * Function using a modified least squares approach to circle fitting.
 *
 * Reference :
 *
 * Umbach, D. and Jones, K. N. "A few methods for fitting circles to data", 
 * IEEE Trans Instrumentation and Measurement
 * vol XX(Y) 2000
 *
 * https://pdfs.semanticscholar.org/faac/44067f04abf10af7dd583fca0c35c5937f95.pdf
 *
 * NOTES 
 *
 * The code below has not been checked for numerical stability or conditioning.
 */

int circle_fit_MLS (point P[], int n, double *x_pos, double *y_pos, double *radius)

 int i;
 double sum_x=0.0, sum_y=0.0, sum_xx=0.0, sum_xy=0.0, sum_yy=0.0, sum_xyy=0.0, sum_yxx=0.0, sum_xxx=0.0, sum_yyy=0.0;
 double A, B, C, D, E;
 double x2, y2, xy, F, xdif, ydif;


 for (i=0; i<n; i++)
 
    sum_x += P[i].x;
    sum_y += P[i].y;
    x2 = P[i].x * P[i].x;
    y2 = P[i].y * P[i].y;
    xy = P[i].x * P[i].y;
    sum_xx += x2;
    sum_yy += y2;
    sum_xy += xy;
    sum_xyy += xy*P[i].y;
    sum_yxx += P[i].y*x2;
    sum_xxx += x2*P[i].x;
    sum_yyy += y2*P[i].y;
 

 A = n * sum_xx - sum_x * sum_x;
 B = n * sum_xy - sum_x * sum_y;
 C = n * sum_yy - sum_y * sum_y;
 D = 0.5 * ( n * (sum_xyy + sum_xxx) - sum_x * sum_yy - sum_x * sum_xx);
 E = 0.5 * ( n * (sum_yxx + sum_yyy) - sum_y * sum_xx - sum_y * sum_yy);

 F = A*C - B*B;
 *x_pos = (D*C - B*E) / F;
 *y_pos = (A*E - B*D) / F;

 *radius = 0;
 for (i=0; i<n; i++)
 
     xdif = P[i].x - *x_pos;
     ydif = P[i].y - *y_pos;
     *radius += sqrt(xdif * xdif + ydif * ydif);
 
 *radius /= n;

 return 0;

下面的主程序可用于测试代码。请发回任何结果/观察/建议以改进 cmets。

int main()

 point *P;
 int    n, i;
 double xpos, ypos, radius;

 printf ("Please enter the number of points \n> ");
 scanf ("%d", &n);
 P = malloc (n * sizeof(point));

 for (i=0; i<n; i++)
 
     printf ("x%d = ", i);
     scanf ("%f", &P[i].x);
     printf ("y%d = ", i);
     scanf ("%f", &P[i].y);
 

 circle_fit_MLS (P, n, &xpos, &ypos, &radius);

 printf (" a = %f\n", xpos);
 printf (" b = %f\n", ypos);
 printf (" r = %f\n", radius);

【讨论】：

答案已更新，包括圆拟合的算法和代码。