hdu5733内切球

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