求三维空间中的三角形外接圆圆心坐标的算法

基本思想：过外接圆的圆心与每条边的中点的直线垂直于每条边，对每条边建立这样的方程，联立方程组，然后解之即得。这里要注意Cramer法则中必须排除系数行列式为零的情况，即要排除三个顶点在一条直线上的情形。

下面是实现代码

//两个向量之差
static void get_vector_diff( float3& aimV, const float3 a, const float3 b )
{
    aimV[0] = b[0] - a[0];

    aimV[1] = b[1] - a[1];

    aimV[2] = b[2] - a[2];
}


//两个向量的内积
static float get_inner_product(const float3 v1,const float3 v2)
{
return v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2];
}


//三阶矩阵的行列式
static float get_vector3_det( const float3 v1, const float3 v2, const float3 v3)
{
    float a[3][3];
    for ( int i = 0; i != 3; ++i )
    {
        a[0][i] = v1[i];
        a[1][i] = v2[i];
        a[2][i] = v3[i];
    }

    return a[0][0] * a[1][1] * a[2][2] + 
a[0][1] * a[1][2] * a[2][0] + 
a[0][2] * a[1][0] * a[2][1] - 
a[0][2] * a[1][1] * a[2][0] - 
a[0][1] * a[1][0] * a[2][2] - 
a[0][0] * a[1][2] * a[2][1];
}

//欧氏范数
static float get_Euclidean_distance(const float3 v1,const float3 v2)
{
return sqrt(pow(v1[0]-v2[0],2)+pow(v1[1]-v2[1],2)+pow(v1[2]-v2[2],2));
}

static void get_vector_center( float3& aimV, const float3 a, const float3 b )
{
aimV[0] = (a[0] + b[0])/2;

    aimV[1] = (a[1] + b[1])/2;

    aimV[2] = (a[2] + b[2])/2;
}


///////////////////////////////////////////  
//求三维空间的三角形外接圆圆心坐标及半径  
///////////////////////////////////////////  
static bool triangle_circumcentre(float3 center,float *radiu,float3 pt[3])
{
float3 line_center1,line_center2,line_center3;

float3 normal1,normal2,normal3;

float3 result,rep1,rep2,rep3;

float d_det;

get_vector_center(line_center1,pt[0],pt[1]);

get_vector_center(line_center2,pt[0],pt[2]);

get_vector_center(line_center3,pt[1],pt[2]);

get_vector_diff(normal1,pt[0],pt[1]);

get_vector_diff(normal2,pt[0],pt[2]);

get_vector_diff(normal3,pt[1],pt[2]);


result[0]=get_inner_product(line_center1,normal1);

result[1]=get_inner_product(line_center2,normal2);

result[2]=get_inner_product(line_center3,normal3);

rep1[0]=normal1[0];

rep1[1]=normal2[0];

rep1[2]=normal3[0];


rep2[0]=normal1[1];

rep2[1]=normal2[1];

rep2[2]=normal3[1];


rep3[0]=normal1[2];

rep3[1]=normal2[2];

rep3[2]=normal3[2];

d_det=get_vector3_det(rep1,rep2,rep3);

if(d_det==0)
{
return false;
}

//Cramer法则
center[0]=get_vector3_det(result,rep2,rep3)/d_det;

center[1]=get_vector3_det(rep1,result,rep3)/d_det;

center[2]=get_vector3_det(rep1,rep2,result)/d_det;

*radiu=get_Euclidean_distance(center,pt[0]);

return true;
}