HDU 5733

HDU 5733 - tetrahedron

题意：
    给定四点，求是否能够成四面体，若能则求出其内接圆心和半径

是否能构成四面体： 三点成面的法线和另一点与三点中任一点相连的向量是否垂直？

四面体内接球
    球心： 任意三个角平分面的交点
    半径： 交点到任意面的距离

#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
const double EPS = 1e-8;
struct Point3
{
    double x,y,z;
    Point3() {}
    Point3(double x,double y,double z):x(x),y(y),z(z) {}
};
typedef Point3 Vec3;
Vec3 operator + (Vec3 A,Vec3 B)
{ 
    return Vec3(A.x + B.x, A.y + B.y, A.z + B.z); 
}
Vec3 operator - (Vec3 A,Vec3 B)
{ 
    return Vec3(A.x - B.x, A.y - B.y, A.z - B.z); 
}
Vec3 operator * (Vec3 A,double p)
{ 
    return Vec3(A.x * p, A.y * p, A.z * p);
}
Vec3 operator / (Vec3 A,double p)
{ 
    return Vec3(A.x / p, A.y / p, A.z / p); 
}
int dcmp(double x)//double cmp
{
    return fabs(x) < EPS ? 0 : (x < 0? -1: 1);
}
double Dot3(Vec3 A,Vec3 B)//Dot mult
{
    return A.x*B.x + A.y*B.y + A.z*B.z;
}
double Length3(Vec3 A)
{
    return sqrt( Dot3(A, A) );
}
double Angle(Vec3 A,Vec3 B)//夹角 
{
    return acos(Dot3(A,B) / Length3(A) / Length3(B));
}
Vec3 Cross3(Vec3 A,Vec3 B)//叉积 右手螺旋A->B 
{
    return Vec3(A.y * B.z - A.z * B.y,
                A.z * B.x - A.x * B.z,
                A.x * B.y - A.y * B.x);
}
struct Plane
{
    Vec3 n; //法线
    Point3 p0; 
    Plane() {}
    Plane(Vec3 nn,Point3 pp0)
    {
        n = nn / Length3(nn);
        p0 = pp0;
    }
    Plane(Point3 a,Point3 b,Point3 c)
    {
        Point3 nn = Cross3(b-a,c-a);
        n = nn/(Length3(nn));
        p0 = a;
    }
};
//角平分面 法向量为两平面法向量相加(内角)或相减(外角)
Plane jpfPlane(Point3 a1,Point3 a2,Point3 b,Point3 c) 
{
    Plane p1(a1,b,c),p2(a2,c,b);
    Vec3 temp = p1.n + p2.n;
    return Plane(Cross3(temp, c-b),b);
}
Point3 LinePlaneIntersection(Point3 p1,Point3 p2,Plane a)//线面交点 取线上任意两点 
{
    Point3 p0 = a.p0;
    Vec3 n = a.n,v = p2-p1;
    double t = (Dot3(n,p0-p1) / Dot3(n,p2-p1)); //映射到法向量的比例 
    return p1 + v * t;
}
Point3 PlaneInsertion(Plane a,Plane b,Plane c)//三面交点 
{//两面交线与两面法线均垂直，法线叉积为其方向矢量. 
    Vec3 nn = Cross3(a.n,b.n),use = Cross3(nn,a.n);
    Point3 st = LinePlaneIntersection(a.p0, a.p0+use,b);//得交线上一点 
    return LinePlaneIntersection(st, st+nn, c);
}
double DistanceToPlane(Point3 p,Plane a)
{
    Point3 p0 = a.p0;
    Vec3 n = a.n;
    return fabs( Dot3(p-p0,n) / Length3(n) );
} 
bool isOnePlane(Point3 a,Point3 b,Point3 c,Point3 d)
{
    double t = Dot3(d-a,Cross3(b-a,c-a) );
    return dcmp(t)==0;
}
int main()
{
    Point3 p[4];
    while(~scanf("%lf%lf%lf",&p[0].x,&p[0].y,&p[0].z))
    {
        for(int i=1;i<4;i++) scanf("%lf%lf%lf",&p[i].x,&p[i].y,&p[i].z);
        if(isOnePlane(p[0],p[1],p[2],p[3]))
        {
            puts("O O O O"); continue;
        }
        Plane a = jpfPlane(p[3],p[2],p[1],p[0]),//三个角平分面的交点即为圆心 
              b = jpfPlane(p[3],p[0],p[1],p[2]),
              c = jpfPlane(p[3],p[1],p[0],p[2]);
        Plane d(p[0],p[1],p[2]);
        Point3 center = PlaneInsertion(a,b,c);
        double r = DistanceToPlane(center,d);
        printf("%.4f %.4f %.4f %.4f
",center.x,center.y,center.z,r);
    }
}