HDU #5733 tetrahedron

tetrahedron

传送门

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

Problem Description
Given four points ABCD, if ABCD is a tetrahedron, calculate the inscribed sphere of ABCD.
 
Input
Multiple test cases .
Each test cases contains a line of 12 integers indicating 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
 
 
Solution:
求四面体的内切球的半径和球心坐标。
半径可以通过将体积算两次来求：第一次用向量算，第二次用四个面的面积和乘内切球半径算。
内心的直角坐标可用体积坐标来算。
四面体的体积坐标

$设四面体的四个顶点分别为A_1, A_2, A_3, A_4, 对于空间内任一点P, 我们用vec{P}表示vec{OP}$
$对 extbf{四面体内}任意一点P, 有$

[vec P =sum_{i=1}^{4}lambda_ivec A_i,]

[sum_{i=1}^{4}lambda_i=1]

$四面体体积坐标的几何意义:$

$其各坐标分量是以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_ivec{A_i}=frac{sumlimits_{i=1}^{4}S_ivec{A_i}}{sumlimits_{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}]

无解的情况对应着四面体四点共面, 即体积为零.

Implementation:

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;

struct point{  
    LL x,y,z; 
    int read(){
        return scanf("%lld%lld%lld", &x, &y, &z);
    }

    point operator-(const point &p){
        return {x-p.x, y-p.y, z-p.z};
    }

    point operator^(const point &p){ //cross product
        return {y*p.z-z*p.y, z*p.x-x*p.z, x*p.y-y*p.x};
    }
    double operator*(const point &p){    //dot product
        return x*p.x+y*p.y+z*p.z;
    }
    double len(){
        return sqrt(x*x+y*y+z*z);
    }

}p[4];


int main(){
    for(; ~p[0].read(); ){
        for(int i=1; i<4; i++) p[i].read();

        LL vol=abs(((p[1]-p[0])^(p[2]-p[0]))*(p[3]-p[0]));

        if(!vol){
            puts("O O O O");    //error-prone: O, not 0
            continue;
        }

        double s[4], sum=0;
        point vec[2];

        for(int i=0; i<4; i++){
           for(int j=0; j<4; j++)
                if(j!=i){
                    for(int k=0, l=0; k<4; k++)
                        if(k!=j && k!=i)
                            vec[l++]=p[k]-p[j];
                    break;
                }
            s[i]=abs((vec[0]^vec[1]).len()), sum+=s[i];
        }

        double tot=0, x, y, z;
        for(int i=0; i<4; i++) tot+=s[i]*p[i].x;
        x=tot/sum;
        tot=0;
        for(int i=0; i<4; i++) tot+=s[i]*p[i].y;
        y=tot/sum;
        tot=0;
        for(int i=0; i<4; i++) tot+=s[i]*p[i].z;
        z=tot/sum;
        printf("%.4f %.4f %.4f %.4f
", x, y, z, vol/sum);
    }
}