HDU 5733 - tetrahedron

HDU 5733 - tetrahedron

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

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

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

 

  1 #include <iostream>
  2 #include <cstdio>
  3 #include <cmath>
  4 using namespace std;
  5 const double EPS = 1e-8;
  6 struct Point3
  7 {
  8     double x,y,z;
  9     Point3() {}
 10     Point3(double x,double y,double z):x(x),y(y),z(z) {}
 11 };
 12 typedef Point3 Vec3;
 13 Vec3 operator + (Vec3 A,Vec3 B)
 14 { 
 15     return Vec3(A.x + B.x, A.y + B.y, A.z + B.z); 
 16 }
 17 Vec3 operator - (Vec3 A,Vec3 B)
 18 { 
 19     return Vec3(A.x - B.x, A.y - B.y, A.z - B.z); 
 20 }
 21 Vec3 operator * (Vec3 A,double p)
 22 { 
 23     return Vec3(A.x * p, A.y * p, A.z * p);
 24 }
 25 Vec3 operator / (Vec3 A,double p)
 26 { 
 27     return Vec3(A.x / p, A.y / p, A.z / p); 
 28 }
 29 int dcmp(double x)//double cmp
 30 {
 31     return fabs(x) < EPS ? 0 : (x < 0? -1: 1);
 32 }
 33 double Dot3(Vec3 A,Vec3 B)//Dot mult
 34 {
 35     return A.x*B.x + A.y*B.y + A.z*B.z;
 36 }
 37 double Length3(Vec3 A)
 38 {
 39     return sqrt( Dot3(A, A) );
 40 }
 41 double Angle(Vec3 A,Vec3 B)//夹角 
 42 {
 43     return acos(Dot3(A,B) / Length3(A) / Length3(B));
 44 }
 45 Vec3 Cross3(Vec3 A,Vec3 B)//叉积 右手螺旋A->B 
 46 {
 47     return Vec3(A.y * B.z - A.z * B.y,
 48                 A.z * B.x - A.x * B.z,
 49                 A.x * B.y - A.y * B.x);
 50 }
 51 struct Plane
 52 {
 53     Vec3 n; //法线
 54     Point3 p0; 
 55     Plane() {}
 56     Plane(Vec3 nn,Point3 pp0)
 57     {
 58         n = nn / Length3(nn);
 59         p0 = pp0;
 60     }
 61     Plane(Point3 a,Point3 b,Point3 c)
 62     {
 63         Point3 nn = Cross3(b-a,c-a);
 64         n = nn/(Length3(nn));
 65         p0 = a;
 66     }
 67 };
 68 //角平分面 法向量为两平面法向量相加(内角)或相减(外角)
 69 Plane jpfPlane(Point3 a1,Point3 a2,Point3 b,Point3 c) 
 70 {
 71     Plane p1(a1,b,c),p2(a2,c,b);
 72     Vec3 temp = p1.n + p2.n;
 73     return Plane(Cross3(temp, c-b),b);
 74 }
 75 Point3 LinePlaneIntersection(Point3 p1,Point3 p2,Plane a)//线面交点 取线上任意两点 
 76 {
 77     Point3 p0 = a.p0;
 78     Vec3 n = a.n,v = p2-p1;
 79     double t = (Dot3(n,p0-p1) / Dot3(n,p2-p1)); //映射到法向量的比例 
 80     return p1 + v * t;
 81 }
 82 Point3 PlaneInsertion(Plane a,Plane b,Plane c)//三面交点 
 83 {//两面交线与两面法线均垂直，法线叉积为其方向矢量. 
 84     Vec3 nn = Cross3(a.n,b.n),use = Cross3(nn,a.n);
 85     Point3 st = LinePlaneIntersection(a.p0, a.p0+use,b);//得交线上一点 
 86     return LinePlaneIntersection(st, st+nn, c);
 87 }
 88 double DistanceToPlane(Point3 p,Plane a)
 89 {
 90     Point3 p0 = a.p0;
 91     Vec3 n = a.n;
 92     return fabs( Dot3(p-p0,n) / Length3(n) );
 93 } 
 94 bool isOnePlane(Point3 a,Point3 b,Point3 c,Point3 d)
 95 {
 96     double t = Dot3(d-a,Cross3(b-a,c-a) );
 97     return dcmp(t)==0;
 98 }
 99 int main()
100 {
101     Point3 p[4];
102     while(~scanf("%lf%lf%lf",&p[0].x,&p[0].y,&p[0].z))
103     {
104         for(int i=1;i<4;i++) scanf("%lf%lf%lf",&p[i].x,&p[i].y,&p[i].z);
105         if(isOnePlane(p[0],p[1],p[2],p[3]))
106         {
107             puts("O O O O"); continue;
108         }
109         Plane a = jpfPlane(p[3],p[2],p[1],p[0]),//三个角平分面的交点即为圆心 
110               b = jpfPlane(p[3],p[0],p[1],p[2]),
111               c = jpfPlane(p[3],p[1],p[0],p[2]);
112         Plane d(p[0],p[1],p[2]);
113         Point3 center = PlaneInsertion(a,b,c);
114         double r = DistanceToPlane(center,d);
115         printf("%.4f %.4f %.4f %.4f\n",center.x,center.y,center.z,r);
116     }
117 }

 

posted @ 2016-07-25 23:33  nicetomeetu  阅读(165)  评论(0)  编辑  收藏  举报