HDU 5733 tetrahedron（计算几何）

原题链接

Problem Description

Given four points ABCD, if ABCD is a tetrahedron, calculate the inscribed sphere of ABCD.

Input

Multiple test cases (test cases ≤100).

Each test cases contains a line of 12 integers [−1e6,1e6] indicate the coordinates of four vertices of ABCD.

Input ends by EOF.

Output

Print the coordinate of the center of the sphere and the radius, rounded to 4 decimal places.

If there is no such sphere, output “O O O O”.

Sample Input

0 0 0 2 0 0 0 0 2 0 2 0
0 0 0 2 0 0 3 0 0 4 0 0

Sample Output

0.4226 0.4226 0.4226 0.4226
O O O O

题目大意

给出四个坐标，首先判断是否是四面体，如不是，直接输出O O O O，否则的话输出它的内切球坐标及球半径。

解题思路

有公式r=v*3/s（v为四面体体积，s为四面体表面积），所以半径是比较好求的。求坐标的话可以利用公式

ansx=( Sabc*a[4].x+Sabd*a[3].x + Sacd*a[2].x+Sbcd*a[1].x)/S; 
ansy=( Sabc*a[4].y+Sabd*a[3].y + Sacd*a[2].y+Sbcd*a[1].y)/S; 
ansz=( Sabc*a[4].z+Sabd*a[3].z + Sacd*a[2].z+Sbcd*a[1].z)/S;

似乎用VJ拉题的时候会把’O’搞成’0’？（于是WA了一万发）

AC代码

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cctype>
#include<algorithm>
#include<cmath>
#include<vector>
#include<string>
#include<queue>
#include<list>
#include<stack>
#include<set>
#include<map>
#define ll long long
#define ull unsigned long long
#define rep(i,a,b) for (int i=(a),_ed=(b);i<_ed;i++)
#define fil(a,b) memset((a),(b),sizeof(a))
#define cl(a) fil(a,0)
#define pb push_back
#define mp make_pair
#define PI 3.1415927
#define inf 0x3f3f3f3f
#define fi first
#define se second
#define eps 1e-8
#define sign(x) ((x)>eps?1:((x)<-eps?(-1):(0)))
using namespace std;

struct point3
{
    double x,y,z;
    point3(){}
    point3(double _x,double _y,double _z)
    {
        x=_x,y=_y,z=_z;
    }
    point3 operator-(const point3 &ne)
    {
        return point3(x-ne.x,y-ne.y,z-ne.z);
    }
    point3 operator+(const point3 &ne)
    {
        return point3(x+ne.x,y+ne.y,z+ne.z);
    }
    point3 operator*(const double t)
    {
        return point3(x*t,y*t,z*t);
    }
};
struct line3
{
    point3 a,b;
    line3(){}
    line3(point3 _a,point3 _b)
    {
        a=_a;
        b=_b;
    }
};
struct plane3
{
    point3 a,b,c;
    plane3(){}
    plane3(point3 _a,point3 _b,point3 _c)
    {
        a=_a;
        b=_b;
        c=_c;
    }
};
point3 xmult(point3 a,point3 b)
{
    return point3(a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y -a.y*b.x);
}
double dmult(point3 a,point3 b)
{
    return a.x*b.x+a.y*b.y+a.z*b.z;
}
double length(point3 v)
{
    return sqrt(v.x*v.x+v.y*v.y+v.z*v.z);
}
double dist(point3 a,point3 b)
{
    return length(a-b);
}
double dist2(point3 a,point3 b)
{
    return dmult(a-b,a-b);
}
//平面法向量
point3 pvec(plane3 s)   
{     
    return xmult(s.b-s.a,s.c-s.a); 
} 
//判断点是否在平面上
bool point_on_plane(point3 p,plane3 s)
{     
    return sign(dmult(p-s.a,pvec(s)))==0; 
}
//四面体体积
double volume(point3 a,point3 b,point3 c,point3 d)
{     //abc 面方向与 d一致时为正     
    return fabs(dmult(xmult(b-a,c-a),d-a))/6.0;
} 
double area3(point3 a,point3 b,point3 c)
{
    double q=dist(a,b);
    double w=dist(b,c);
    double e=dist(a,c);
    double s=(q+w+e)/2;
    return sqrt(s*(s-q)*(s-w)*(s-e));
}
int main(void)
{
    point3 a[5];
    double v;
    double s1,s2,s3,s4;
    double r;
    double rx,ry,rz;
    while(scanf("%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",&a[1].x,&a[1].y,&a[1].z,&a[2].x,&a[2].y,&a[2].z,
        &a[3].x,&a[3].y,&a[3].z,&a[4].x,&a[4].y,&a[4].z)!=EOF)
    {
        plane3 p1(a[1],a[2],a[3]);
        if(point_on_plane(a[4],p1))
        {
            printf("O O O O\n");
        }
        else
        {
            v=abs(volume(a[1],a[2],a[3],a[4]));
            s1=area3(a[1],a[2],a[3]);
            s2=area3(a[1],a[2],a[4]);
            s3=area3(a[1],a[3],a[4]);
            s4=area3(a[2],a[3],a[4]);
            r=v*3/(s1+s2+s3+s4);
            if(r<1e-9) printf("O O O O\n");
            else{
                rx=(s1*a[4].x+s2*a[3].x+s3*a[2].x+s4*a[1].x)/(s1+s2+s3+s4);
                ry=(s1*a[4].y+s2*a[3].y+s3*a[2].y+s4*a[1].y)/(s1+s2+s3+s4);
                rz=(s1*a[4].z+s2*a[3].z+s3*a[2].z+s4*a[1].z)/(s1+s2+s3+s4);
                printf("%.4f %.4f %.4f %.4f\n",rx,ry,rz,r);
            }
        }
    }
    return 0;
}