HDU 5733 tetrahedron （求四面体内切球模板）

转自：点击打开链接

给四个点让求内接圆心。

就求呗~

内心公式：

设四面体A1A2A3A4的顶点Ai多对的侧面积为Si（i=1,2,3,4）,顶点Ai的坐标为（xi,yi,zi）(i=1,2,3,4),四面体内心I的坐标为（xi,yi,zi）,则

x1=(s1*x1+s2*x2+s3*x3+s4*x4)/(s1+s2+s3+s4);

y1=(s1*y1+s2*y2+s3*y3+s4*y4)/(s1+s2+s3+s4);

z1=(s1*z1+s2*z2+s3*z3+s4*z4)/(s1+s2+s3+s4);

外心公式：

任取三个点对，两点间的中垂面焦点即为外接圆心。

下面上代码，求面积的时候不要用海伦公式，用向量法。不然吃精度，算不出来。

#include <iostream>#include <cstdio>#include <cstring>#include <cmath>#include <algorithm>#include <string>#include <cstdlib>#include <vector>#include <set>#include <map>using namespace std;const double EPS = 1e-6;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 - (Point3 A,Point3 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){    return fabs(x) < EPS ? 0 : (x <0 ? -1:1);}double Dot3(Vec3 A, Vec3 B){    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){    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 = Plane(a1, b, c),p2 = Plane(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 isOnPlane(Point3 a,Point3 b,Point3 c,Point3 d){    double t = Dot3(d-a,Cross3(b-a, c-a));    return dcmp(t)==0;}int main(){    //freopen("in.txt", "r", stdin);    // freopen("out.txt", "w", stdout);    int x,y,z;    Point3 p[4];    while(scanf("%d%d%d",&x,&y,&z)==3)    {        p[0] = Point3((double)x,(double)y,(double)z);        for(int i=1; i<4; i++)        {            scanf("%d%d%d",&x,&y,&z);            p[i] = Point3((double)x,(double)y,(double)z);        }        if(isOnPlane(p[0],p[1],p[2],p[3]))        {            puts("O O O O");        }        else        {            Plane a = jpfPlane(p[0],p[1],p[2],p[3]),                  b = jpfPlane(p[1],p[2],p[3],p[0]),                  c = jpfPlane(p[2],p[3],p[0],p[1]);            Point3 center = PlaneInsertion(a, b, c);            double r = DistanceToPlane(center, Plane(p[0],p[1],p[2]));            printf("%.3f %.3f %.3f %.3f\n",center.x,center.y,center.z,r);        }    }    return 0;}

//已知3点坐标，求平面ax+by+cz+d=0; void get_panel(Point p1,Point p2,Point p3,double &a,double &b,double &c,double &d){ a = ( (p2.y-p1.y)*(p3.z-p1.z)-(p2.z-p1.z)*(p3.y-p1.y) );    b = ( (p2.z-p1.z)*(p3.x-p1.x)-(p2.x-p1.x)*(p3.z-p1.z) );    c = ( (p2.x-p1.x)*(p3.y-p1.y)-(p2.y-p1.y)*(p3.x-p1.x) );    d = ( 0-(a*p1.x+b*p1.y+c*p1.z) );}// 已知三点坐标，求法向量Vec3 get_Normal(Point p1,Point p2,Point p3){ a = ( (p2.y-p1.y)*(p3.z-p1.z)-(p2.z-p1.z)*(p3.y-p1.y) );    b = ( (p2.z-p1.z)*(p3.x-p1.x)-(p2.x-p1.x)*(p3.z-p1.z) );    c = ( (p2.x-p1.x)*(p3.y-p1.y)-(p2.y-p1.y)*(p3.x-p1.x) );    return Vec3(a,b,c);}//点到平面距离 double dis_pt2panel(Point pt,double a,double b,double c,double d){ return f_abs(a*pt.x+b*pt.y+c*pt.z+d)/sqrt(a*a+b*b+c*c);}

  #include <cstdio>      #include <cstdlib>      #include <cstring>      #include <algorithm>      #include <cmath>            using namespace std;            const double EPS = 1e-9;      const int MAXN = 40;            struct Point3  //空间点      {          double x, y, z;          Point3( double x=0, double y=0, double z=0 ): x(x), y(y), z(z) { }          Point3( const Point3& a )          {              x = a.x;              y = a.y;              z = a.z;              return;          }          void showP()          {              printf("%f %f %f \n", x, y, z);          }          Point3 operator+( Point3& rhs )          {              return Point3( x+rhs.x, y+rhs.y, z+rhs.z );          }      };            struct Line3   //空间直线      {          Point3 a, b;      };            struct plane3   //空间平面      {          Point3 a, b, c;          plane3() {}          plane3( Point3 a, Point3 b, Point3 c ):              a(a), b(b), c(c) { }          void showPlane()          {              a.showP();              b.showP();              c.showP();              return;          }      };            double dcmp( double a )      {          if ( fabs( a ) < EPS ) return 0;          return a < 0 ? -1 : 1;      }            //三维叉积      Point3 Cross3( Point3 u, Point3 v )      {          Point3 ret;          ret.x = u.y * v.z - v.y * u.z;          ret.y = u.z * v.x - u.x * v.z;          ret.z = u.x * v.y - u.y * v.x;          return ret;      }            //三维点积      double Dot3( Point3 u, Point3 v )      {          return u.x * v.x + u.y * v.y + u.z * v.z;      }            //矢量差      Point3 Subt( Point3 u, Point3 v )      {          Point3 ret;          ret.x = u.x - v.x;          ret.y = u.y - v.y;          ret.z = u.z - v.z;          return ret;      }            //两点距离      double TwoPointDistance( Point3 p1, Point3 p2 )      {          return sqrt( (p1.x - p2.x)*(p1.x - p2.x) + (p1.y - p2.y)*(p1.y - p2.y) + (p1.z - p2.z)*(p1.z - p2.z) );      }            //向量的模      double VectorLenth( Point3 p )      {          return sqrt( p.x*p.x + p.y*p.y + p.z*p.z );      }            //空间直线距离      double LineToLine( Line3 u, Line3 v, Point3& tmp )      {          tmp = Cross3( Subt( u.a, u.b ), Subt( v.a, v.b ) );          return fabs( Dot3( Subt(u.a, v.a), tmp ) ) / VectorLenth(tmp);      }            //取平面法向量      Point3 pvec( plane3 s )      {          return Cross3( Subt( s.a, s.b ), Subt( s.b, s.c ) );      }            //空间平面与直线的交点      Point3 Intersection( Line3 l, plane3 s )      {          Point3 ret = pvec(s);          double t = ( ret.x*(s.a.x-l.a.x)+ret.y*(s.a.y-l.a.y)+ret.z*(s.a.z-l.a.z) )/( ret.x*(l.b.x-l.a.x)+ret.y*(l.b.y-l.a.y)+ret.z*(l.b.z-l.a.z) );          ret.x = l.a.x + ( l.b.x - l.a.x ) * t;          ret.y = l.a.y + ( l.b.y - l.a.y ) * t;          ret.z = l.a.z + ( l.b.z - l.a.z ) * t;          return ret;      }