hdu5733内切球

题意：给定任意四面体的四个顶点坐标，求四面体内接球的中心坐标和半径，内接球可以不存在。
设四个点是A、B、C、D ,那么内切球的圆心可以表示为

#include<cstdio>
#include<cmath>

struct point  //存点
{
    double x, y, z;
}p[5];

struct plane   //存平面方程系数
{
    double a, b, c, d;
} pm[5];

double area[5];  //记录每个面的面积

bool check(point p0, point p1, point p2)  //检查两向量是否平行
{
    p1.x -= p0.x, p1.y -= p0.y, p1.z -= p0.z;
    p2.x -= p0.x, p2.y -= p0.y, p2.z -= p0.z;
    if((p1.x * p2.y == p1.y * p2.x) && (p1.x * p2.z==p1.z * p2.x) && (p1.y * p2.z==p1.z * p2.y))  return false;
    return true;
}

double pointoplaneDist(double x, double y, double z, double a, double b, double c, double d)  //点到平面的距离
{
    return fabs(a * x + b * y + c * z + d) / sqrt(a * a + b * b + c * c);
}

double dist(point p0, point p1)  //两点距离
{
    return sqrt((p0.x - p1.x) * (p0.x - p1.x) + (p0.y - p1.y) * (p0.y - p1.y) + (p0.z - p1.z) * (p0.z - p1.z));
}

point cross(point p0,point p1)  //叉乘
{
    point p2;
    p2.x = p0.y * p1.z - p0.z * p1.y;
    p2.y = -p0.x * p1.z+p0.z * p1.x;
    p2.z = p0.x * p1.y - p0.y * p1.x;
    return p2;
}

void cal(int cnt,point p0,point p1,point p2)  //计算平面方程系数
{
    double a = dist(p0,p1);
    double b = dist(p1,p2);
    double c = dist(p0,p2);
    double d = (a + b + c) / 2.0;
    area[cnt] = sqrt(d * (d - a) * (d - b) * (d - c));
    p1.x -= p0.x, p1.y -= p0.y, p1.z -= p0.z;
    p2.x -= p0.x, p2.y -= p0.y, p2.z -= p0.z;
    point p3 = cross(p1,p2);
    pm[cnt].a = p3.x, pm[cnt].b = p3.y, pm[cnt].c = p3.z;
    pm[cnt].d = -(pm[cnt].a * p0.x + pm[cnt].b * p0.y + pm[cnt].c * p0.z);
}

int main()
{
    while(scanf("%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",&p[0].x,&p[0].y,&p[0].z,&p[1].x,&p[1].y,&p[1].z,&p[2].x,&p[2].y,&p[2].z,&p[3].x,&p[3].y,&p[3].z)!=EOF)
    {
        if(!(check(p[0],p[1],p[2]) && check(p[0],p[2],p[3]) && check(p[0],p[1],p[3]) && check(p[1],p[2],p[3])))  //判断四点是否共面
        {
            printf("O O O O\n");
            continue;
        }
        cal(3, p[0], p[1], p[2]);
        cal(2, p[0], p[1], p[3]);
        cal(1, p[0], p[2], p[3]);
        cal(0, p[1], p[2], p[3]);
        double r = area[3] * pointoplaneDist(p[3].x, p[3].y, p[3].z, pm[3].a, pm[3].b, pm[3].c, pm[3].d)/(area[0] + area[1] + area[2] + area[3]);
        if(r < 1e-9)
        {
            printf("O O O O\n");
            continue;
        }
        point ans;
        ans.x = (area[0] * p[0].x + area[1] * p[1].x + area[2] * p[2].x + area[3] * p[3].x) / (area[0] + area[1] + area[2] + area[3]);
        ans.y = (area[0] * p[0].y + area[1] * p[1].y + area[2] * p[2].y + area[3] * p[3].y) / (area[0] + area[1] + area[2] + area[3]);
        ans.z = (area[0] * p[0].z + area[1] * p[1].z + area[2] * p[2].z + area[3] * p[3].z) / (area[0] + area[1] + area[2] + area[3]);
        printf("%.4lf %.4lf %.4lf %.4lf\n",ans.x, ans.y, ans.z, r);
    }
    return 0;
}