hdu5733 （计算几何，全是公式）

tetrahedron

Time Limit: 2000/1000 MS (Java/Others)    Memory Limit: 65536/65536 K (Java/Others)

Total Submission(s): 408    Accepted Submission(s): 157

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
  
  

 

Author

HIT

 

Source

2016 Multi-University Training Contest 1

 

题目大意：给四面体四个顶点的左边，求该四面体的内切球的球心坐标及其半径，如果不存在内切球，则输出O O O O。

解题思路：各种公式的套用，首先要判断是否存在内切球，可以转换成给的四个顶点是否可容易构成四面体。我们先要求出三个点构成的平面方程，然后判断另一点是否在这个平面上。

接着要求每个面的面积，再此之前还要求一下每条边的边长，再根据海伦公式求面积。

海伦公式为： S = sqrt(s1* (s1-a) * (s1-b)*(s1-c))，其中s1表示的是三角形的周长的一半。

最后求内切球的球心，公式为：

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; 

其中S为表面积。

#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<vector>
#include<map>
#include<set>
#include<queue>
#include<stack>
#include<iostream>
using namespace std;
typedef long long LL;

#define eps 1e-8
#define zero(x) (((x)>0?(x):-(x))<eps)
struct point3
{
    double x,y,z;
};
struct line3
{
    point3 a,b;
};
struct plane3
{
    point3 a,b,c;
};

//计算 cross product product product product U x V
point3 xmult(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;
}
//计算 dot product product product product U . V
double dmult(point3 u,point3 v)
{
    return u.x*v.x+u.y*v.y+u.z*v.z;
}
//矢量差 U - V
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;
}
//取平面法向量
point3 pvec(plane3 s)
{
    return xmult(subt(s.a,s.b),subt(s.b,s.c));
}
point3 pvec(point3 s1,point3 s2,point3 s3)
{
    return xmult(subt(s1,s2),subt(s2,s3));
}
//两点距离,单参数取向量大小
double dis(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 vlen(point3 p)
{
    return sqrt(p.x*p.x+p.y*p.y+p.z*p.z);
}

int dots_onplane(point3 a,point3 b,point3 c,point3 d)
{
    return zero(dmult(pvec(a,b,c),subt(d,a)));
}

//点到平面距离
double ptoplane(point3 p,point3 s1,point3 s2,point3 s3)
{
    return fabs(dmult(pvec(s1,s2,s3),subt(p,s1)))/vlen(pvec(s1,s2,s3));
}

double Area(point3 a,point3 b,point3 c)
{
    double ab=dis(a,b),bc=dis(b,c),ac=dis(a,c);
    double p=(ab+bc+ac)/2;
    return sqrt(p* (p-ab) * (p-bc)*(p-ac));
}

int main()
{
    point3 p[5];
    while(~scanf("%lf%lf%lf",&p[1].x,&p[1].y,&p[1].z))
    {
        for(int i=2;i<=4;i++) scanf("%lf%lf%lf",&p[i].x,&p[i].y,&p[i].z);
        if(dots_onplane(p[1],p[2],p[3],p[4])) { printf("O O O O\n"); continue; }
        double Sabc=Area(p[1],p[2],p[3]);
        double Sabd=Area(p[1],p[2],p[4]);
        double Sacd=Area(p[1],p[3],p[4]);
        double Sbcd=Area(p[2],p[3],p[4]);
        double S=Sabc+Sabd+Sacd+Sbcd;
        point3 heart;
        heart.x=( Sabc*p[4].x+Sabd*p[3].x + Sacd*p[2].x+Sbcd*p[1].x)/S;
        heart.y=( Sabc*p[4].y+Sabd*p[3].y + Sacd*p[2].y+Sbcd*p[1].y)/S;
        heart.z=( Sabc*p[4].z+Sabd*p[3].z + Sacd*p[2].z+Sbcd*p[1].z)/S;
        printf("%.4lf %.4lf %.4lf %.4lf\n",heart.x,heart.y,heart.z,ptoplane(heart,p[1],p[2],p[3]));
    }
    return 0;
}