2016 Multi-University Training Contest 1-1011---HDU 5733 tetrahedron（计算几何）

题目链接

HDU 5733

题意：

给出一个四面体的四个点，求内切球的半径和圆心。

题解：

设四面体的四个顶点分别为A1, A2, A3, A4。

四面体内切球半径：

四面体的总体积：

V=VPA2A3A4+VPA1A3A4+VPA1A2A4+VPA1A2A3

=R3(SA2A3A4+SA1A3A4+SA1A2A4+SA1A2A3)

=R3(S1+S2+S3+S4)

内切球半径：

R=3∗V(SA2A3A4+SA1A3A4+SA1A2A4+SA1A2A3)

=3∗V(S1+S2+S3+S4)

由任意四面体的体积公式可有:

V=|(A1−A4)⋅[(A2−A4)×(A3−A4)]|6

又由三个向量混合积的几何意义：

R=|(A1−A4)⋅[(A2−A4)×(A3−A4)]|/2(S1+S2+S3+S4)

又由两个向量叉积的几何意义：

S1=(A2−A4)⋅(A3−A4)/2

S2=(A4−A1)⋅(A3−A1)/2

S3=(A4−A1)⋅(A2−A1)/2

S4=(A3−A1)⋅(A2−A4)/2

所以：

R=|(A1−A4)⋅[(A2−A4)×(A3−A4)]|((A2−A4)⋅(A3−A4)+(A4−A1)⋅(A3−A1)+(A4−A1)⋅(A2−A1)+(A3−A1)⋅(A2−A4))

四面体内心坐标：

对四面体内任意一点P, 我们用P⃗ 表示OP−→−， 有

P⃗ =∑i=14λiA⃗ i=λ1A⃗ 1+λ2A⃗ 2+λ3A⃗ 3+λ4A⃗ 4

其中0≤λ1,λ2,λ3,λ4≤1 ,且∑4i=1λi=λ1+λ2+λ3+λ4=1 ,则称(λ1,λ2,λ3,λ4)为 P 点关于四面体A1A2A3A4 的体积坐标 。

体积坐标具有明显的几何意义, 其各坐标分量是以 P 为顶点, 以各底面为底的四面体体积与原四面体体积之比。即:

λ1=VPA2A3A4VA1A2A3A4

λ2=VPA1A3A4VA1A2A3A4

λ3=VPA1A2A4VA1A2A3A4

λ4=VPA1A2A3VA1A2A3A4

记四面体A1A2A3A4的四个底面的面积分别为S1,S2,S3,S4, 若P是四面体A1A2A3A4的内心I,因为内心到四个面的距离相等，则有

λi=SiS1+S2+S3+S4,i=1,2,3,4

故

OI→=∑i=14λiAi→=∑i=14SiAi→∑i=14Si

从而I的直角坐标(x, y, z)为:

x=∑i=14Sixi∑i=14Si

y=∑i=14Siyi∑i=14Si

z=∑i=14Sizi∑i=14Si

与上面一样面积可以表示成向量叉积的形式：

S1=(A2−A4)⋅(A3−A4)/2

S2=(A4−A1)⋅(A3−A1)/2

S3=(A4−A1)⋅(A2−A1)/2

S4=(A3−A1)⋅(A2−A4)/2

代入计算即可。

代码

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
const double eps=1e-8;
int dcmp(double a)
{
if(fabs(a)<eps) return 0;
return a>0?1:-1;
}
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;
}
};
typedef point3 vector3;
vector3 operator +(vector3 a,vector3 b)
{
return vector3(a.x+b.x,a.y+b.y,a.z+b.z);
}
vector3 operator -(vector3 a,vector3 b)
{
return vector3(a.x-b.x,a.y-b.y,a.z-b.z);
}

vector3 operator *(vector3 a,double k)
{
return vector3(a.x*k,a.y*k,a.z*k);
}

vector3 operator /(vector3 a,double k)
{
return vector3(a.x/k,a.y/k,a.z/k);
}

vector3 Cross3(vector3 u,vector3 v)
{
point3 ret;
ret.x = u.y * v.z - u.z * v.y;
ret.y = u.z * v.x - u.x * v.z;
ret.z = u.x * v.y - u.y * v.x;
return ret;
}

double Dot3( vector3 u, vector3 v )
{
return u.x * v.x + u.y * v.y + u.z * v.z;
}
double len(vector3 a)
{
return sqrt(Dot3(a,a));
}

int main()
{
//freopen("in.txt","r",stdin);
//freopen("out.txt","w",stdout);
point3 a,b,c,d;
while(cin>>a.x>>a.y>>a.z>>b.x>>b.y>>b.z>>c.x>>c.y>>c.z>>d.x>>d.y>>d.z)
{
if(dcmp(Dot3(b-a,Cross3(c-a,d-a)))==0) printf("O O O O\n");
else
{
vector3 v1=Cross3(b-a,c-a);
vector3 v2=Cross3(b-a,d-a);
vector3 v3=Cross3(c-a,d-a);
vector3 v4=Cross3(c-b,d-b);
double r=fabs(Dot3(a-d,Cross3(b-d,c-d)))/(len(v1)+len(v2)+len(v3)+len(v4));
double x=(d.x*len(v1)+c.x*len(v2)+b.x*len(v3)+a.x*len(v4))/(len(v1)+len(v2)+len(v3)+len(v4));
double y=(d.y*len(v1)+c.y*len(v2)+b.y*len(v3)+a.y*len(v4))/(len(v1)+len(v2)+len(v3)+len(v4));
double z=(d.z*len(v1)+c.z*len(v2)+b.z*len(v3)+a.z*len(v4))/(len(v1)+len(v2)+len(v3)+len(v4));
printf("%.4lf %.4lf %.4lf %.4lf\n",x,y,z,r);
}
}
return 0;
}