HDU 5733 tetrahedron (2016 Multi-University Training Contest 1 计算几何)

传送门
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
题目大意：
给定4个点，如果是四面体的话求内接球，输出球心及其半径，如果不存在内接球，输出 Ｏ Ｏ O O。
解题思路：

都是公式，首先我们判断是不是四面体，也就是判断这四个点是不是共面，我们需要求一个3个点构成平面的方程，然后判断另一个点是不是在这个平面上。然后我们再来求这个四面体的内接球的半径r，根据公式：r=V∗3/S(其中S=s1+s2+s3+s4,si表示的是4个面的面积，也就是S是四面体的表面积)那么我们来求一下四面体的体积，根据公式，四面体的体积:V=[AB,AD,AC]/6,其中[AB,AD,AC]表示的是AB,AC,AD的混合积，也就是AB与AC的叉积在与AD的点积

然后再来求每个面的面积，我们可以先求一下每条边的边长，然后再根据海伦公式：S = sqrt(s1* (s1-a) * (s1-b)*(s1-c))，其中s1表示三角形周长的一半 求得每个面的面积。然后我们再求一下球心O：内切球的球心坐标就是根据公式：
ansx=( Sabc*a[4].x+Sabd*a[3].x + Sacd*a[2].x+Sbcd*a[1].x)/S;
ansy=( Sabc*a[4].y+Sabd*a[3].y + Sacd*a[2].y+Sbcd*a[1].y)/S;
ansz=( Sabc*a[4].z+Sabd*a[3].z + Sacd*a[2].z+Sbcd*a[1].z)/S;
S为表面积。
My Code:

#include <iostream>#include <cstdio>#include <cstdlib>#include <cmath>using namespace std;const double eps = 1e-6;struct point{    double x, y, z;} a[10];struct V{    point st;    point end;};point Cross(point a, point b){    point s;    s.x = a.y*b.z - a.z*b.y;    s.y = a.z*b.x - a.x*b.z;    s.z = a.x*b.y - b.x*a.y;    return s;}double dot(point a, point b){    double s = a.x*b.x+a.y*b.y+a.z*b.z;    return s;}int main(){    while(cin>>a[1].x>>a[1].y>>a[1].z)    {        for(int i=2; i<=4; i++)            cin>>a[i].x>>a[i].y>>a[i].z;        ///求平面方程 A*x + B*y + C*z + D = 0        double A = ((a[2].y-a[1].y)*(a[3].z-a[1].z)-(a[2].z-a[1].z)*(a[3].y-a[1].y));        double B = ((a[2].z-a[1].z)*(a[3].x-a[1].x)-(a[2].x-a[1].x)*(a[3].z-a[1].z));        double C = ((a[2].x-a[1].x)*(a[3].y-a[1].y)-(a[2].y-a[1].y)*(a[3].x-a[1].x));        double D = -(A * a[1].x + B * a[1].y + C * a[1].z);        ///cout<<A<<" "<<B<<" "<<C<<" "<<D<<endl;        double ret = A*a[4].x+B*a[4].y+a[4].z*C+D;        if(fabs(ret) < eps)        {            puts("O O O O");            continue;        }        V AB, AC, BC, AD, BD, CD;        AB.st = a[1], AB.end = a[2];        AC.st = a[1], AC.end = a[3];        BC.st = a[2], BC.end = a[3];        AD.st = a[1], AD.end = a[4];        BD.st = a[2], BD.end = a[4];        CD.st = a[3], CD.end = a[4];        double aa = sqrt( (AB.st.x-AB.end.x)*(AB.st.x-AB.end.x)+(AB.st.y-AB.end.y)*(AB.st.y-AB.end.y)+                          (AB.st.z-AB.end.z)*(AB.st.z-AB.end.z) );        double b = sqrt( (AC.st.x-AC.end.x)*(AC.st.x-AC.end.x)+(AC.st.y-AC.end.y)*(AC.st.y-AC.end.y)+                         (AC.st.z-AC.end.z)*(AC.st.z-AC.end.z) );        double c = sqrt( (BC.st.x-BC.end.x)*(BC.st.x-BC.end.x)+(BC.st.y-BC.end.y)*(BC.st.y-BC.end.y)+                         (BC.st.z-BC.end.z)*(BC.st.z-BC.end.z) );        double d = sqrt( (AD.st.x-AD.end.x)*(AD.st.x-AD.end.x)+(AD.st.y-AD.end.y)*(AD.st.y-AD.end.y)+                         (AD.st.z-AD.end.z)*(AD.st.z-AD.end.z) );        double e = sqrt( (BD.st.x-BD.end.x)*(BD.st.x-BD.end.x)+(BD.st.y-BD.end.y)*(BD.st.y-BD.end.y)+                         (BD.st.z-BD.end.z)*(BD.st.z-BD.end.z) );        double f = sqrt( (CD.st.x-CD.end.x)*(CD.st.x-CD.end.x)+(CD.st.y-CD.end.y)*(CD.st.y-CD.end.y)+                         (CD.st.z-CD.end.z)*(CD.st.z-CD.end.z) );        double s1 = 0.5*(aa+b+c);        double Sabc = sqrt(s1*(s1-aa)*(s1-b)*(s1-c));        double s2 = 0.5*(aa+e+d);        double Sabd = sqrt(s2*(s2-aa)*(s2-e)*(s2-d));        double s3 = 0.5*(d+b+f);        double Sacd = sqrt(s3*(s3-d)*(s3-b)*(s3-f));        double s4 = 0.5*(e+f+c);        double Sbcd = sqrt(s4*(s4-c)*(s4-e)*(s4-f));        double S = Sabc + Sabd + Sacd + Sbcd;        point AD1,BD1,CD1;        AD1.x = AD.end.x-AD.st.x;        AD1.y = AD.end.y-AD.st.y;        AD1.z = AD.end.z-AD.st.z;        BD1.x = BD.end.x-BD.st.x;        BD1.y = BD.end.y-BD.st.y;        BD1.z = BD.end.z-BD.st.z;        CD1.x = CD.end.x-CD.st.x;        CD1.y = CD.end.y-CD.st.y;        CD1.z = CD.end.z-CD.st.z;        double r = fabs(dot(Cross(AD1,BD1),CD1))/2.0/S;        double ansx, ansy, ansz;        ansx=( Sabc*a[4].x+Sabd*a[3].x + Sacd*a[2].x+Sbcd*a[1].x)/S;        ansy=( Sabc*a[4].y+Sabd*a[3].y + Sacd*a[2].y+Sbcd*a[1].y)/S;        ansz=( Sabc*a[4].z+Sabd*a[3].z + Sacd*a[2].z+Sbcd*a[1].z)/S;        printf("%.4lf %.4lf %.4lf %.4lf\n",ansx, ansy, ansz, r);    }    return 0;}