给定三个点，如何找到圆的中心

Question

我翻译的结果如下：

我阅读了这个链接，并编写了相应的代码，但在链接中解释的示例中得到了错误的答案。在解方程时，我从方程1中减去方程2，从方程2中减去方程3，然后继续进行。请查看链接以获取更多澄清。

我的代码是:

include<stdio.h>
int is_formCircle(float a1,float b1,float a2,float b2,float a3,float b3) {
  float p1=a2-a1;
  float p2=a3-a2;
  float p3=b2-b1;
  float p4=b3-b2;
  float alpha=(a1+a2)*(a1-a2) + (b1+b2)*(b1-b2);
  float beta =(a2+a3)*(a2-a3) + (b2+b3)*(b2-b3);
  float y1=p1*beta - p2*alpha;
  float y2=p2*p3 - p1*p4;
  if(y2==0 || y1==0) return 1;
  float y=y1/y2;
  float x1 = 2*p4*y + beta;
  float x2 = 2*p2;
  float x = x1/x2;
  printf("x=%f  y=%f\n",x,y);
  return 0;
}
int main() {
 float a1,a2,a3,a4,b1,b2,b3,b4;
 a1=4.0;
 b1=1.0;
 a2=-3.0;
 b2=7.0;
 a3=5.0;
 b3=-2.0;
 is_formCircle(a1,b1,a2,b2,a3,b3);
 return 0;
}

我的另一个代码：

#include<stdio.h>
int is_formCircle(float a1,float b1,float a2,float b2,float a3,float b3) {                  
  float mid1,mid2,mid3,mid4,m1,m2,D,Dx,Dy,x,y;
  mid1 = a1+(a2-a1)/2;
  mid2 = b1+(b2-b1)/2;
  mid3 = a2+(a3-a2)/2;
  mid4 = b2+(b3-b2)/2;
  m1=(b2-b1)/(a2-a1);
  m2=(b3-b2)/(a3-a2);
  m1=-1*m1;
  m2=-1*m2;
  D=m2-m1;
  Dx=mid2-(m1*mid1) + (mid3*m2) - mid4;
  Dy=(m1*(mid3*m2-mid4))-(m2*(mid1*m1-mid2));
  x=Dx/D;
  y=Dy/D;
  printf("%f %f",x,y);    
  return 0;
}
int main() {
 float a1,a2,a3,b1,b2,b3;      
 a1=4.0;
 b1=1.0;
 a2=-3.0;
 b2=7.0;
 a3=5.0;
 b3=-2.0;
 is_formCircle(a1,b1,a2,b2,a3,b3);      
 return 0;
}

为什么我的代码会出现错误答案？

- avinashse

---

请编写易读的代码。即使是比赛，您如何调试它？ - zw324

typedef double f;太糟糕了！毫无道理。 - Daniel Fischer

@avinashse typedef float f; 不是更好的选择。类型有名称是有原因的。请使用它们。 - Daniel Fischer

@Aravind 对于可读性不佳，我很抱歉。如果我必须更改数据类型，我使用了那个typedef，这样我只需要在一行中更改，如果它造成了不便，那么很抱歉:( - avinashse

@PatriciaShanahan 是的，typedef double float_type;（或者如果你喜欢的话，可以用 real）是可以的，你可以给类型取一个有意义的名字。typedef whatever a; 则不同。 - Daniel Fischer

Answer 1

在链接中给出的解决方案是一种“盲目”的解决方案，也就是说，你知道方程，就可以解决它。但是，如果你更深刻地理解了背后的原理，你就能够：

1.编写更易读、可靠、灵活的代码。

2.轻松调试。

当你从方程1中减去方程2时会发生什么？实际上，你试图找到描述那些与点1和点2等距离的点的直线方程。然后，你对点2和3进行相同的操作。最后，在这些线的交点处，你找到了圆的中心。

如何描述与点1和2等距离的点的直线？你取在两者之间的点，并沿着垂直于点1和2之间方向的方向前进。

如果这不是非常清楚，请拿起一张纸画一个例子：放置点1、2和3，找到两条线并找到交点。

现在，既然你理解了所有的内容，请用两个函数重新编写你的代码，一个用于查找两点之间的等距离线，另一个用于计算两条线之间的交点……

修改后，你的代码看起来更好了，尽管理解起来并不简单。我认为错误在于你求解两条直线交点时，没有忘记它是参数方程式。

Dx = (mid4-mid2) - m2*(mid3-mid1);

lambda=Dx/D;

x = mid1 + lambda*m1;
y = mid2 + lambda*1.0;

使用Matlab进行图形化检查。

- Dr_Sam

---

谢谢回复。我编辑了我的代码，也许做了你说的同样的事情，找到了两点的中点和垂直于连接这两点的线段的斜率，然后形成了该垂直平分线的方程，接下来解决了它以找到圆的中心。 - avinashse

谢谢您的回复，我有一个问题。我的任务是要检查第四个点是在圆内、外还是在圆上。例如，假设有四个点A、B、C、D。我正在寻找圆心A、B、C，并计算半径(R)，然后计算|AD|。如果(|AD| > R)，则点D在圆外，那么是否有其他更快的方法来完成它呢？ - avinashse

@avinashse 注意，你应该计算圆心K的|KD|，而不是|AD|。否则，看起来很好。你也可以测试|KD|^2 > R^2，这样就不需要平方根，但这不是一个很大的优化。 - Dr_Sam

哦，是的，应该是KD树...抱歉我的错误。我想问的是，这是找出一个点是否在圆内的唯一方法吗？如果假设有30个点，并且我必须检查每三个点中有多少个点在圆上或圆内，那么这将是解决问题的最佳方式吗？即，可以形成30C3 = 4060个圆，并且对于每个圆，我必须检查其他27个未包含在圆形成中的点是否在圆内，如果在圆内，则增加计数。是否有其他更快的方法来做到这一点？ - avinashse

@avinashse 嗯，这是一个比较困难的问题。我没有看到任何更好的方法或捷径可以采取。这与特定问题有关还是只是一个测试用例？ - Dr_Sam

Answer 2

我必须说，如果你按照你所列的链接进行操作，保持变量名相同将会更有帮助。我们能够更好地理解算法，看到x1、y1、x2、y2、x3和y3，而不是p1、p2、p3、p4、alpha和beta。实际上，我在你的算法中没有看到与链接匹配的内容。我不想像评论那样苛刻（如果你担心从float转换为double，这是一个完全合适的typedef案例），但是当你不必转换变量名时，调试算法最容易。

我建议仅使用链接中提供的h和k，这是通过计算3x3矩阵的行列式来完成的。你可以在这里找到很多参考资料。

我会编写两个函数，如下所示：

float calculateH(float x1, float y1, float x2, float y2, float x3, float y3) {
    float numerator = (x2*x2+y2*y2)*y3 - (x3*x3+y3*y3)*y2 - 
                      ((x1*x1+y1*y1)*y3 - (x3*x3+y3*y3)*y1) +
                      (x1*x1+y1*y1)*y2 - (x2*x2+y2*y2)*y1;
    float denominator = (x2*y3-x3*y2) -
                        (x1*y3-x3*y1) +
                        (x1*y2-x2*y1);
    denominator *= 2;
    return numerator / denominator;
}
float calculateK(float x1, float y1, float x2, float y2, float x3, float y3) {
    float numerator = x2*(x3*x3+y3*y3) - x3*(x2*x2+y2*y2) -
                      (x1*(x3*x3+y3*y3) - x3*(x1*x1+y1*y1)) +
                      x1*(x2*x2+y2*y2) - x2*(x1*x1+y1*y1);
    float denominator = (x2*y3-x3*y2) -
                        (x1*y3-x3*y1) +
                        (x1*y2-x2*y1);
    denominator *= 2;
    return numerator / denominator;
}

那么你的 is_formCircle 就会变成这样：

float is_formCircle(float x1, float y1, float x2, float y2, float x3, float y3) {
    float h = calculateH(x1, y1, x2, y2, x3, y3);
    float k = calculateK(x1, y1, x2, y2, x3, y3);
    printf("x=%f  y=%f\n",h,k);
}

有很多方法可以优化这个过程，我可能会在行列式计算中出现笔误，但这应该能帮助您开始。

- Scott Mermelstein

---

谢谢。我之前用纸和笔做了同样的事情，并将其公式化到我的代码中，但是它并没有起作用，而这个方法却完美地解决了问题。 - avinashse

不用谢。我注意到你的代码实际上是一个布尔检查，用于判断是否可以从三个点形成一个圆。答案（如Dr_Sam所解释的那样）总是“是”-除非这三个点共线。如果它们共线，在此计算中分母的行列式将计算为0。 - Scott Mermelstein

是的，我考虑过这个问题，除非三个点共线，否则它们总是形成一个外接圆。我的做法是要检查第四个点是在圆内，圆外还是圆上。例如，假设有四个点A、B、C、D。我先找到圆心A、B、C，然后计算半径(R)，再计算|AD|的长度。如果(|AD| > R)，那么点D就在圆外。那么有没有更快的方法呢..？ - avinashse

@avinashse 这是一个更有趣的问题，你应该将其作为另一个问题提出 - 不是因为我在追求声望，而是因为这个问题已经得到解答，很少有人会关注它。（这可能在math.stackexchange.com上更好。）顺便说一下，我有一个相对快速的答案，但我不确定它是否足够。只需遍历您的 4060 种三元组组合，并计算它们之间的距离。我没有证明，但我认为你想要的是使 |AB| + |BC| 最大化的圆 ABC。那应该包括最多的点。 - Scott Mermelstein

使用（分子/分母）时会有精度问题吗？ - avinashse

@avinashse 我可能会惹恼任何关注的人，但我的答案是：根据您应用程序的需求，很可能没有问题。有些人说float一直很糟糕，但它对于大约6-8个有效数字是好的。在大多数情况下，这已经足够了。这就是你开始使用的typedef double f会派上用场的地方 - 如果你发现精度不够，把所有的浮点数改成双精度。 - Scott Mermelstein