HDU 573 tetrahedron 三维计算几何（叉积， 四面体内切球）

tetrahedron

Problem Description

Given four points ABCD, if ABCD is a tetrahedron, calculate the inscribed sphere of ABCD.

 

Input

Multiple test cases (test cases  ≤100 ).

Each test cases contains a line of 12 integers  [−1e6,1e6]  indicate 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
  
  

 

参考：http://blog.csdn.net/mikezhw/article/details/51960523

设四面体四个顶点坐标分别为  (xi,yi,zi) ，与之相对的面的面积为 si ，四面体体积为 V ，则有 
四面体内心（内切球球心）坐标 

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪x = s1x1+s2x2+s3x3+s4x4s1+s2+s3+s4y = s1y1+s2y2+s3y3+s4y4s1+s2+s3+s4z = s1z1+s2z2+s3z3+s4z4s1+s2+s3+s4

四面体内切球半径 

r = 3Vs1+s2+s3+s4

面积可以用三维叉积求，体积则用混合积 
若体积为0，则内切球不存在

#include <bits/stdc++.h>
#define endl "\n"
using namespace std;
typedef long long ll;

const int MAXN = 500000 + 7;
const double eps = 1e-7;
const double PI = acos(-1.0);

int sgn(double x) {
    if(fabs(x) < eps) return 0;
    return x > 0 ? 1 : -1;
}

struct Point {
    double x, y, z;
    Point() {}
    Point(double x, double y, double z) : x(x), y(y), z(z) {}
};

typedef Point Vec;

Vec operator + (Vec v1, Vec v2) {return Vec(v1.x + v2.x, v1.y + v2.y, v1.z + v2.z);}
Vec operator - (Vec v1, Vec v2) {return Vec(v1.x - v2.x, v1.y - v2.y, v1.z - v2.z);}
Vec operator * (Vec v, double p) {return Vec(v.x * p, v.y * p, v.z * p);}
Vec operator / (Vec v, double p) {return Vec(v.x / p, v.y / p, v.z / p);}
bool operator == (Vec v1, Vec v2) {
    return sgn(v1.x - v2.x) == 0 && sgn(v1.y - v2.y) == 0;
}

double vDiot(Vec v1, Vec v2) {return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z;}
double vLen(Vec v) {return sqrt(vDiot(v, v));}
Vec vCross(Vec v1, Vec v2) {return Vec(v1.y * v2.z - v1.z * v2.y, v1.z * v2.x - v1.x * v2.z, v1.x * v2.y - v1.y * v2.x);}
double getArea(Vec a, Vec b, Vec c) {
    return vLen(vCross(b - a, c - a)) / 2;
}
double getVolum(Vec v1, Vec v2, Vec v3, Vec v4) {
    return fabs(vDiot(v4 - v1, vCross(v2 - v1, v3 - v1)) / 6.0);
}

int main() {
    ios::sync_with_stdio(false);
    Point p[4];
    double X, Y, Z, V, R, sa, sb, sc, sd;
    while(cin >> p[0].x >> p[0].y >> p[0].z >> p[1].x >> p[1].y >> p[1].z >> p[2].x
              >> p[2].y >> p[2].z >> p[3].x >> p[3].y >> p[3].z) {
        sa = getArea(p[1], p[2], p[3]);
        sb = getArea(p[0], p[2], p[3]);
        sc = getArea(p[0], p[1], p[3]);
        sd = getArea(p[0], p[1], p[2]);
        X = (sa * p[0].x + sb * p[1].x + sc * p[2].x + sd * p[3].x) / (sa + sb + sc + sd);
        Y = (sa * p[0].y + sb * p[1].y + sc * p[2].y + sd * p[3].y) / (sa + sb + sc + sd);
        Z = (sa * p[0].z + sb * p[1].z + sc * p[2].z + sd * p[3].z) / (sa + sb + sc + sd);
        V = getVolum(p[0], p[1], p[2], p[3]);
        R = 3 * V / (sa + sb + sc + sd);
        if(V < eps) {
            cout << "O O O O" << endl;
        }
        else {
            cout << fixed << setprecision(4) << X << " " << Y << " " << Z << " " << R << endl;
        }
    }
    return 0;
}