簡介

還是太菜從解析官方程式碼開始。
官方程式碼共有三個核心函式。

參考連結

http://blog.sina.com.cn/s/blog_6029f0330101irlh.html
http://paulbourke.net/papers/triangulate/triangulate.c 官方程式碼 ？
https://www.geeksforgeeks.org/equation-of-circle-when-three-points-on-the-circle-are-given/ 求圓心座標的程式碼

首先關於二維的知道三點如何求解圓心

設定方程為 \(A(x^2 + y^2) + Dx + Ey + F = 0\)
如果知道不在一條直線上的三點\((x_1,y_1),(x_2,y_2),(x_3,y_3)\)

,那麼我們可以使用克拉默法則來求解聯立方程組

\[\left\{\begin{array}{l} A\left(x^{2}+y^{2}\right)+D x+E y+F=0 \\ A\left(x_{1}^{2}+y_{1}^{2}\right)+D x_{1}+E y_{1}+F=0 \\ A\left(x_{2}^{2}+y_{2}^{2}\right)+D x_{2}+E y_{2}+F=0 \\ A\left(x_{3}^{2}+y_{3}^{2}\right)+D x_{3}+E y_{3}+F=0 \end{array}\right. \]

則個構成關於A,D,E,F的四元線性方程組，有非零解的充要條件是係數行列式的值為0

\[\left|\begin{array}{llll} x^{2}+y^{2} & x & y & 1 \\ x_{1}^{2}+y_{1}^{2} & x_{1} & y_{1} & 1 \\ x_{2}^{2}+y_{2}^{2} & x_{2} & y_{2} & 1 \\ x_{3}^{2}+y_{3}^{2} & x_{3} & y_{3} & 1 \end{array}\right|=0 \]

（至於為什麼，有人知道可以發在評論區裡面，我忘了）。

師弟的解答：如果不為0的話那麼ADEF都為0一定是其中的解，那麼就沒有意義了。

當任意點(x,y)與已知三點的某一點重合和，則有兩行完全相同行列式為0，表明已知三點滿足這個方程。如果這樣感覺應該有無窮多個解，四個未知量，只有三個方程組。
若將係數行列式按照第一行展開可得

\[\left|\begin{array}{lll} x_{1} & y_{1} & 1 \\ x_{2} & y_{2} & 1 \\ x_{3} & y_{3} & 1 \end{array}\right|\left(x^{2}+y^{2}\right)-\left|\begin{array}{lll} x_{1}^{2}+y_{1}^{2} & y_{1} & 1 \\ x_{2}^{2}+y_{2}^{2} & y_{2} & 1 \\ x_{3}^{2}+y_{3}^{2} & y_{3} & 1 \end{array}\right| x+\left|\begin{array}{lll} x_{1}^{2}+y_{1}^{2} & x_{1} & 1 \\ x_{2}^{2}+y_{2}^{2} & x_{2} & 1 \\ x_{3}^{2}+y_{3}^{2} & x_{3} & 1 \end{array}\right| y-\left|\begin{array}{lll} x_{1}^{2}+y_{1}^{2} & x_{1} & y_{1} \\ x_{2}^{2}+y_{2}^{2} & x_{2} & y_{2} \\ x_{3}^{2}+y_{3}^{2} & x_{3} & y_{3} \end{array}\right|=0 \]

四個係數行列式依次對應著A，-D，E，-F的值。最特別的是三個方程求解出了四個未知量？其實不是的ADEF是呈現一個內在的比例的

\[\left|\begin{array}{lll} x_{1} & y_{1} & 1 \\ x_{2} & y_{2} & 1 \\ x_{3} & y_{3} & 1 \end{array}\right| \neq 0 \]

對應著三點不共線。

QU

1. 太神奇了為啥三個函式可以解除四個未知量。
   最特別的是三個方程求解出了四個未知量？其實不是的ADEF是呈現一個內在的比例的關係。

那如何求解圓心座標和半徑呢？

發現官網的程式碼有點問題。

code

/*
 * @Author: your name
 * @Date: 2020-12-09 15:47:38
 * @LastEditTime: 2020-12-09 16:03:39
 * @LastEditors: Please set LastEditors
 * @Description: In User Settings Edit
 * @FilePath: /delaunay/test.cc
 */

// C++ implementation of the approach 
#include <bits/stdc++.h> 
using namespace std; 
  
// Function to find the circle on 
// which the given three points lie 
void findCircle(int x1, int y1, int x2, int y2, int x3, int y3) 
{ 
    int x12 = x1 - x2; 
    int x13 = x1 - x3; 
  
    int y12 = y1 - y2; 
    int y13 = y1 - y3; 
  
    int y31 = y3 - y1; 
    int y21 = y2 - y1; 
  
    int x31 = x3 - x1; 
    int x21 = x2 - x1; 
  
    // x1^2 - x3^2 
    int sx13 = pow(x1, 2) - pow(x3, 2); 
  
    // y1^2 - y3^2 
    int sy13 = pow(y1, 2) - pow(y3, 2); 
  
    int sx21 = pow(x2, 2) - pow(x1, 2); 
    int sy21 = pow(y2, 2) - pow(y1, 2); 
  
    int f = ((sx13) * (x12) 
             + (sy13) * (x12) 
             + (sx21) * (x13) 
             + (sy21) * (x13)) 
            / (2 * ((y31) * (x12) - (y21) * (x13))); 
    int g = ((sx13) * (y12) 
             + (sy13) * (y12) 
             + (sx21) * (y13) 
             + (sy21) * (y13)) 
            / (2 * ((x31) * (y12) - (x21) * (y13))); 
  
    int c = -pow(x1, 2) - pow(y1, 2) - 2 * g * x1 - 2 * f * y1; 
  
    // eqn of circle be x^2 + y^2 + 2*g*x + 2*f*y + c = 0 
    // where centre is (h = -g, k = -f) and radius r 
    // as r^2 = h^2 + k^2 - c 
    int h = -g; 
    int k = -f; 
    int sqr_of_r = h * h + k * k - c; 
  
    // r is the radius 
    float r = sqrt(sqr_of_r); 
  
    cout << "Centre = (" << h << ", " << k << ")" << endl; 
    cout << "Radius = " << r; 
} 
const double EPSILON=0.00001;
int CircumCircle(double xp,double yp,
   double x1,double y1,double x2,double y2,double x3,double y3,
   double *xc,double *yc,double *rsqr)
{
   double m1,m2,mx1,mx2,my1,my2;
   double dx,dy,drsqr;
   double fabsy1y2 = fabs(y1-y2);
   double fabsy2y3 = fabs(y2-y3);

   /* Check for coincident points */
   if (fabsy1y2 < EPSILON && fabsy2y3 < EPSILON)
       return(false);

   if (fabsy1y2 < EPSILON) {
      m2 = - (x3-x2) / (y3-y2);
      mx2 = (x2 + x3) / 2.0;
      my2 = (y2 + y3) / 2.0;
      *xc = (x2 + x1) / 2.0;
      *yc = m2 * (*xc - mx2) + my2;
   } else if (fabsy2y3 < EPSILON) {
      m1 = - (x2-x1) / (y2-y1);
      mx1 = (x1 + x2) / 2.0;
      my1 = (y1 + y2) / 2.0;
      *xc = (x3 + x2) / 2.0;
      *yc = m1 * (*xc - mx1) + my1;
   } else {
      m1 = - (x2-x1) / (y2-y1);
      m2 = - (x3-x2) / (y3-y2);
      mx1 = (x1 + x2) / 2.0;
      mx2 = (x2 + x3) / 2.0;
      my1 = (y1 + y2) / 2.0;
      my2 = (y2 + y3) / 2.0;
      *xc = (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2);
      if (fabsy1y2 > fabsy2y3) {
         *yc = m1 * (*xc - mx1) + my1;
      } else {
         *yc = m2 * (*xc - mx2) + my2;
      }
   }

   dx = x2 - *xc;
   dy = y2 - *yc;
   *rsqr = dx*dx + dy*dy;

   dx = xp - *xc;
   dy = yp - *yc;
   drsqr = dx*dx + dy*dy;

   // Original
   //return((drsqr <= *rsqr) ? TRUE : FALSE);
   // Proposed by Chuck Morris
   return((drsqr - *rsqr) <= EPSILON ? true : false);
}


// Driver code 
int main() 
{ 
    int x1 = 1, y1 = 1; 
    int x2 = 1, y2 = 4; 
    int x3 = 5, y3 = 3; 
    findCircle(x1, y1, x2, y2, x3, y3); // 明顯上面的是對的，下面的出現了錯誤
    std::cout << "\n----------\n";
    double x4, y4, r;
    CircumCircle(0,0, x1, y1, x2, y2, x3, y3, &x4, &y4, &r);
    cout << "Centre = (" << x4 << ", " << y4 << ")" << endl; 
    cout << "Radius = " << r; 
    return 0; 
}