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

1. 题目链接
   HDU 5733
2. 题意：
   给出一个四面体的四个点，求内切球的半径和圆心。
3. 题解：
   设四面体的四个顶点分别为$A_1$, $A_2$, $A_3$, $A_4$。
   
   • 四面体内切球半径：
     四面体的总体积：
     $$V=V_{PA_2 A_3A_4}+V_{PA_1A_3A_4}+V_{PA_1A_2A_4}+V_{PA_1A_2A_3}$$
     $$=\frac{R}{3}(S_{A_2 A_3A_4}+S_{A_1A_3A_4}+S_{A_1A_2A_4}+S_{A_1A_2A_3})$$
     $$=\frac{R}{3}(S_1+S_2+S_3+S_4)$$
     内切球半径：
     $$R=\frac{3*V}{(S_{A_2 A_3A_4}+S_{A_1A_3A_4}+S_{A_1A_2A_4}+S_{A_1A_2A_3})}$$
     $$=\frac{3*V}{(S_1+S_2+S_3+S_4)}$$
     由任意四面体的体积公式可有:
     $$V=\frac { |(A_1-A_4) \cdot [({A_2}-{A_4}) \times ({A_3}-{A_4})]|}{6}$$
     又由三个向量混合积的几何意义：
     $$R= \frac { |(A_1-A_4) \cdot [({A_2}-{A_4}) \times ({A_3}-{A_4})]|/2}{(S_1+S_2+S_3+S_4)}$$
     又由两个向量叉积的几何意义：
     $$S_1=(A_2-A_4)\cdot(A_3-A_4)/2$$
     $$S_2=(A_4-A_1)\cdot(A_3-A_1)/2$$
     $$S_3=(A_4-A_1)\cdot(A_2-A_1)/2$$
     $$S_4=(A_3-A_1)\cdot(A_2-A_4)/2$$
     所以：
     $R= \frac { |(A_1-A_4) \cdot [({A_2}-{A_4}) \times ({A_3}-{A_4})]|}{((A_2-A_4)\cdot(A_3-A_4)+(A_4-A_1)\cdot(A_3-A_1)+(A_4-A_1)\cdot(A_2-A_1)+(A_3-A_1)\cdot(A_2-A_4))}$
   • 四面体内心坐标：
     对四面体内任意一点$P$, 我们用$\vec{P}$表示$\overrightarrow{OP}$， 有
     $$\vec P =\sum_{i=1}^{4}\lambda_i\vec A_i=\lambda_1\vec A_1+\lambda_2\vec A_2+\lambda_3\vec A_3+\lambda_4\vec A_4$$
     其中$0 ≤\lambda_1,\lambda_2,\lambda_3,\lambda_4 ≤1$ ,且$\sum_{i=1}^{4}\lambda_i=\lambda_1+\lambda_2+\lambda_3+\lambda_4=1$ ,则称$(\lambda_1,\lambda_2,\lambda_3,\lambda_4)$为 P 点关于四面体A1A2A3A4 的体积坐标 。
     体积坐标具有明显的几何意义, 其各坐标分量是以 P 为顶点, 以各底面为底的四面体体积与原四面体体积之比。即:
     $$\lambda_1=\frac{V_{PA_2 A_3A_4}}{V_{A_1A_2A_3A_4}}$$
     $$\lambda_2=\frac{V_{PA_1A_3A_4}}{V_{A_1A_2A_3A_4}}$$
     $$\lambda_3=\frac{V_{PA_1A_2A_4}}{V_{A_1A_2A_3A_4}}$$
     $$\lambda_4=\frac{V_{PA_1A_2A_3}}{V_{A_1A_2A_3A_4}}$$
     记四面体$A_1A_2A_3A_4$的四个底面的面积分别为$S_1, S_2, S_3, S_4$, 若P是四面体$A_1A_2A_3A_4$的内心$I$,因为内心到四个面的距离相等，则有
     $$\lambda_i = \frac{S_i}{S_1+S_2+S_3+S_4}, \quad i=1, 2, 3, 4$$
     故
     $$\vec{OI}=\sum_{i=1}^{4}\lambda_i\vec{A_i}=\frac{\sum\limits_{i=1}^{4}S_i\vec{A_i}}{\sum\limits_{i=1}^{4}S_i}$$
     从而I的直角坐标(x, y, z)为:
     $$x=\frac{\sum\limits_{i=1}^{4}S_ix_i} {\sum\limits_{i=1}^{4}S_i}$$
     $$y=\frac{\sum\limits_{i=1}^{4}S_iy_i} {\sum\limits_{i=1}^{4}S_i}$$
     $$z=\frac{\sum\limits_{i=1}^{4}S_iz_i} {\sum\limits_{i=1}^{4}S_i}$$
     与上面一样面积可以表示成向量叉积的形式：
     $$S_1=(A_2-A_4)\cdot(A_3-A_4)/2$$
     $$S_2=(A_4-A_1)\cdot(A_3-A_1)/2$$
     $$S_3=(A_4-A_1)\cdot(A_2-A_1)/2$$
     $$S_4=(A_3-A_1)\cdot(A_2-A_4)/2$$
     代入计算即可。
4. 代码

#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;
}