录十六

已知点P₀(x₀, y₀)，P₁(x₁, y₁)，P₂(x₂, y₂)是三维空间中不共线的三点，且按照顺时针或者逆时针排序。求过三点的圆弧。

一、求圆弧的圆心及半径

（1）平面方程

在三维空间中给定了一点M₀，与两个不共线的向量a、b，那么通过点M₀且与a、b平行的平面π就唯一地被确定，向量a、b叫做平面π的方位向量。

如果设点M₀的坐标为(x₀, y₀, z₀)，并设向量a=(X₁, Y₁, Z₁) 、向量b=(X₂, Y₂, Z₂)，那么平面π的点位式方程为：

已知P₀，P_1，P₂是空间中不共线的三点，则向量a=P₁-P₀和b=P₂-P₀必然不共线，并且是圆弧所在平面的方位向量；取P₀为圆弧所在平面的一个固定点，代入平面的点位式方程有：

设

利用行列式的拆分性质，以及代数余子式展开公式，并整理有：

由于圆心坐标O(o_x, o_y, o_z)位于圆弧所在的平面内，可将圆心坐标代入平面方程，可得：

（2）圆方程

圆心坐标为O(o_x, o_y, o_z)，圆的半径为r，则圆的方程可写为：

将圆弧上的三点P₀，P₁，P₂分别代入圆的方程，可得方程组：

以上方程组消去r后，并化简可得：

其中等式(4)和方程组(7)组合，并写成矩阵相乘的形式，可得：

（3）求解方程组

struct Point3d{
    double x;
    double y;
    double z;

    Point3d()
        :x(0.0),y(0.0),z(0.0)
    {
    }

    Point3d(double xx,double yy,double zz)
        :x(xx),y(yy),z(zz)
    {    
    }
};

bool SloveCircleCenter(const Point3d& p0,const Point3d& p1,const Point3d& p2,
                       Point3d& ptCenter)
{
    double ax = p1.x - p0.x;
    double ay = p1.y - p0.y;
    double az = p1.z - p0.z;

    double bx = p2.x - p0.x;
    double by = p2.y - p0.y;
    double bz = p2.z - p0.z;

    double A = ay*bz - by*az;
    double B = ax*bz - bx*az;
    double C = ax*by - bx*ay;
    double D = A*p0.x - B*p0.y + C*p0.z;

    double p0LenSquare = p0.x*p0.x+ p0.y*p0.y+ p0.z*p0.z;
    double p1LenSquare = p1.x*p1.x+ p1.y*p1.y+ p1.z*p1.z;
    double p2LenSquare = p2.x*p2.x+ p2.y*p2.y+ p2.z*p2.z;

    double a1,b1,c1,d1;
    double a2,b2,c2,d2;
    double a3,b3,c3,d3;

    a1 = A; b1=-B;  c1= C;  d1 = -D;
    a2= ax; b2= ay; c2= az; d2 = 0.5*(p1LenSquare-p0LenSquare);
    a3= bx; b3= by; c3= bz; d3 = 0.5*(p2LenSquare-p0LenSquare);
    
    double dInv = a1*(b2*c3-b3*c2)-a2*(b1*c3-b3*c1)+a3*(b1*c2-b2*c1);
    if(fabs(dInv) > 0.00000000001)
    {
         dInv = 1.0/dInv;
    }
    else
    {
        return false;
    }

    ptCenter.x = dInv*((b2*c3-b3*c2)*d1 + (b3*c1-b1*c3)*d2 + (b1*c2-c1*b2)*d3); 
    ptCenter.y = dInv*((c2*a3-c3*a2)*d1 + (c3*a1-c1*a3)*d2 + (c1*a2-a1*c2)*d3); 
    ptCenter.z = dInv*((a2*b3-a3*b2)*d1 + (a3*b1-a1*b3)*d2 + (a1*b2-b1*a2)*d3); 

    return true;
}