给定三个点时找到圆心

Question

我研究了this link 并进行了相应的编码，但在链接中解释的示例中得到了错误的答案， 在求解方程的过程中，我从方程 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;
}

为什么我的代码给出了错误的答案？

Answer 1

我不得不说，如果您关注您列出的链接，这将有助于保持变量名称相同。看到 x1、y1、x2、y2、x3、y3 而不是 p1、p2、p3、p4、alpha 和 beta，我们可以更好地理解算法。实际上，我在您的算法中看不到与链接匹配的内容。我不想像 cmets 那样苛刻（如果您担心将 float 转换为 double，这对于 typedef 来说是一个非常好的案例），但是当您不必转换时，调试算法是最简单的变量名。

我建议在链接中简单地使用他们为 h 和 k 提供的内容，即通过计算 3x3 矩阵的行列式来完成。你可以找到很多references。

我会做两个函数，如下：

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);
}

有很多方法可以优化这个，我有可能在行列式计算中打错了，但它应该能让你继续前进。

Answer 2

链接中给出的解决方案是一个“盲”解决方案，即，您知道方程式，请快速求解。

但是，如果您更深入地了解幕后情况，您将能够：

1. 编写更易读、更可靠、更灵活的代码。
2. 轻松调试。

从等式 2 中减去等式 1 会发生什么？实际上，您尝试找到描述与点 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 以图形方式检查。