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:
求四面体的内切球的半径和球心坐标。
半径可以通过将体积算两次来求：第一次用向量算，第二次用四个面的面积和乘内切球半径算。
内心的直角坐标可用体积坐标来算。
四面体的体积坐标

$\mathrm{设}\mathrm{四}\mathrm{面}\mathrm{体}\mathrm{的}\mathrm{四}\mathrm{个}\mathrm{顶}\mathrm{点}\mathrm{分}\mathrm{别}\mathrm{为}{A}_1,A_2,A_3,A_4,\mathrm{对}\mathrm{于}\mathrm{空}\mathrm{间}\mathrm{内}\mathrm{任}\mathrm{一}\mathrm{点}P,\mathrm{我}\mathrm{们}\mathrm{用}\overrightarrow{P}\mathrm{表}\mathrm{示}\overrightarrow{OP}$
$\mathrm{对}\mathbf{\text{四面体内}}\mathrm{任}\mathrm{意}\mathrm{一}\mathrm{点}P,\mathrm{有}$

$$\overrightarrow{P}=\sum _{i=1}^4{\lambda }_i{\overrightarrow{A}}_i,$$

$$\sum _{i=1}^4{\lambda }_i=1$$

$\mathrm{四}\mathrm{面}\mathrm{体}\mathrm{体}\mathrm{积}\mathrm{坐}\mathrm{标}\mathrm{的}\mathrm{几}\mathrm{何}\mathrm{意}\mathrm{义}:$

$\mathrm{其}\mathrm{各}\mathrm{坐}\mathrm{标}\mathrm{分}\mathrm{量}\mathrm{是}\mathrm{以}P\mathrm{为}\mathrm{顶}\mathrm{点},\mathrm{以}\mathrm{各}\mathrm{底}\mathrm{面}\mathrm{为}\mathrm{底}\mathrm{的}\mathrm{四}\mathrm{面}\mathrm{体}\mathrm{的}\mathrm{体}\mathrm{积}\mathrm{与}\mathrm{原}\mathrm{四}\mathrm{面}\mathrm{体}\mathrm{的}\mathrm{体}\mathrm{积}\mathrm{之}\mathrm{比}.\mathrm{即}:$

$${\lambda }_1=\frac{V_{P{A}_2{A}_3{A}_4}}{V_{A_1{A}_2{A}_3{A}_4}}$$

$${\lambda }_2=\frac{V_{P{A}_1{A}_3{A}_4}}{V_{A_1{A}_2{A}_3{A}_4}}$$

$${\lambda }_3=\frac{V_{P{A}_1{A}_2{A}_4}}{V_{A_1{A}_2{A}_3{A}_4}}$$

$${\lambda }_4=\frac{V_{P{A}_1{A}_2{A}_3}}{V_{A_1{A}_2{A}_3{A}_4}}$$

$\mathrm{记}\mathrm{四}\mathrm{面}\mathrm{体}{A}_1{A}_2{A}_3{A}_4\mathrm{的}\mathrm{四}\mathrm{个}\mathrm{底}\mathrm{面}\mathrm{的}\mathrm{面}\mathrm{积}\mathrm{分}\mathrm{别}\mathrm{为}{S}_1,S_2,S_3,S_4,\mathrm{若}P\mathrm{是}\mathrm{四}\mathrm{面}\mathrm{体}{A}_1{A}_2{A}_3{A}_4\mathrm{的}\mathrm{内}\mathrm{心}I,\mathrm{则}\mathrm{有}$

$${\lambda }_i=\frac{S_i}{S_1+S_2+S_3+S_4},i=1,2,3,4$$

$\mathrm{故}$

$$\overrightarrow{OI}=\sum _{i=1}^4{\lambda }_i\overrightarrow{A_i}=\frac{\sum _{i=1}^4{S}_i\overrightarrow{A_i}}{\sum _{i=1}^4{S}_i}$$

$\mathrm{从}\mathrm{而}I\mathrm{的}\mathrm{直}\mathrm{角}\mathrm{坐}\mathrm{标}(x,y,z)\mathrm{为}:$

$$x=\frac{\sum _{i=1}^4{S}_i{x}_i}{\sum _{i=1}^4{S}_i}$$

$$y=\frac{\sum _{i=1}^4{S}_i{y}_i}{\sum _{i=1}^4{S}_i}$$

$$z=\frac{\sum _{i=1}^4{S}_i{z}_i}{\sum _{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\n", x, y, z, vol/sum);
    }
}