- 题目链接
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()
{
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;
}