HDU 5733 求四面体 内心 外心 内接圆圆心 外接圆圆心

给四个点让求内接圆心。

就求呗~

内心公式：

设四面体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 <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <string>
#include <cmath>
#include <map>
#include <cstdlib>
using namespace std;

typedef long long LL;
const int MOD = int(1e9)+7;
const int INF = 0x3f3f3f3f;
const double EPS = 1e-9;
const double PI = acos(-1.0); //M_PI;
const int maxn = 100001;


class Point_3
{
public:
    double x , y , z;
    Point_3() {}
    Point_3(double xx , double yy , double zz) : x(xx)  , y(yy) , z(zz) {}
    int input()
    {
        return scanf("%lf%lf%lf",&x,&y,&z);
    }
    friend Point_3 operator  - (const Point_3 &a , const Point_3 &b)
    {
        return Point_3(a.x - b.x , a.y - b.y , a.z - b.z);
    }
};

Point_3 det(const Point_3 &a , const Point_3 &b)
{
    return Point_3(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 dot(const Point_3 &a , const Point_3 &b)
{
    return a.x * b.x + a.y * b.y + a.z * b.z;
}

Point_3 pvec(Point_3 &s1 , Point_3 &s2 , Point_3 s3)
{
    return det((s1 - s2) , (s2 - s3));
}

bool zreo(double x)
{
    return fabs(x) < EPS;
}

int dots_onplane(Point_3  a , Point_3 b  , Point_3 c , Point_3 d )
{
    return zreo(dot(pvec(a , b , c ) , d - a));
}
double dis(Point_3 a,Point_3 b)
{
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)+(a.z-b.z)*(a.z-b.z));
}
void get_panel(Point_3 p1,Point_3 p2,Point_3 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) );
}

double dis_pt2panel(Point_3 pt,double a,double b,double c,double d)
{
    return fabs(a*pt.x+b*pt.y+c*pt.z+d)/sqrt(a*a+b*b+c*c);
}

double fun(Point_3 a,Point_3 b,Point_3 c)
{
    Point_3 u(b.x-a.x,b.y-a.y,b.z-a.z);
    Point_3 v(c.x-a.x,c.y-a.y,c.z-a.z);

    return sqrt(dot(det(u,v),det(u,v)))/2.0;
}

int main()
{
#ifdef DoubleQ
    freopen("in.txt","r",stdin);
#endif
    Point_3 a , b , c , d;
    while(scanf ("%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",&a.x,&a.y,&a.z,&b.x,&b.y,&b.z,&c.x,&c.y,&c.z,&d.x,&d.y,&d.z)!=EOF)
    {
        if(dots_onplane(a , b , c , d))
            printf("O O O O\n");
        else
        {
            double ab=dis(a,b);
            double ac=dis(a,c);
            double ad=dis(a,d);
            double bc=dis(b,c);
            double cd=dis(c,d);
            double db=dis(b,d);

            double sa=fun(b,c,d);
            double sb=fun(a,c,d);
            double sc=fun(a,b,d);
            double sd=fun(a,b,c);

            double x1=(sa*a.x+sb*b.x+sc*c.x+sd*d.x)/(sa+sb+sc+sd);
            double y1=(sa*a.y+sb*b.y+sc*c.y+sd*d.y)/(sa+sb+sc+sd);
            double z1=(sa*a.z+sb*b.z+sc*c.z+sd*d.z)/(sa+sb+sc+sd);

            double aa,bb,cc,dd;
            get_panel(b,c,d,aa,bb,cc,dd);
            double dis=dis_pt2panel(a,aa,bb,cc,dd);
            double r=(sa*dis)/(sa+sb+sc+sd);
            printf ("%.4f %.4f %.4f %.4f\n",x1,y1,z1,r);
        }
    }
    return 0;
}

以下是标程，顺便存一发模板

/*************************************************************************
    > File Name: ZNL.cpp
    > Author: znl1087
    > Mail: loveCJNforever@gmail.com
    > Created Time: 三  7/29 11:14:13 2015
 ************************************************************************/

#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;  
    }  
      
    /************以上模板*************/