四面体内接圆_______tetrahedron( hdu 5733 2016多校第一场 )

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.4226O O O O

 

题意：

给出四个点，问四个点组成的四面体的内接圆圆心和半径，如果没有内接圆就输出O O O O。

分析：

很显然如果没有内接圆一定是无法构成四面体，也就是四点共面，

至于内接圆公式：

o.x = ( Sabc * d.x + Sabd * c.x + Sacd * b.x + Sbcd * a.x) / S

o.y = ( Sabc * d.y + Sabd * c.y + Sacd * b.y + Sbcd * a.y) / S

o.z = ( Sabc * d.z + Sabd * c.z + Sacd * b.z + Sbcd * a.z) / S

S = Sabc + Sbcd + Sabd + Sacd.

r 等于圆心到任何一面的距离。

代码：

#include<stdio.h>#include<string.h>#include<algorithm>#include<cmath>#include<iostream>using namespace std;#define eps 1e-8#define zero(x)(((x)>0?(x):-(x))<eps)struct point3{    double x,y,z;};struct line3{    line3(){}    line3(point3 ta,point3 tb){a = ta,b = tb;}    point3 a,b;    double len(){ return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)+(a.z-b.z)*(a.z-b.z));}};struct plane3{    plane3(){}    plane3(point3 ta,point3 tb,point3 tc){ a = ta,b = tb,c = tc;}    point3 a,b,c;    double s()    {        line3 ab(a,b),ac(a,c),bc(b,c);        double h = (ab.len() + bc.len() + ac.len())/2;        return sqrt(h*(h-ab.len())*(h-ac.len())*(h-bc.len()));    }};point3 xmult(point3 u,point3 v){    point3 ret;    ret.x = u.y*v.z - v.y*u.z;    ret.y = u.z*v.x - u.x*v.z;    ret.z = u.x*v.y - u.y*v.x;    return ret;}double dmult(point3 u,point3 v){    return u.x*v.x+u.y*v.y+u.z*v.z;}point3 subt(point3 u,point3 v){    point3 ret;    ret.x = u.x - v.x;    ret.y = u.y - v.y;    ret.z = u.z - v.z;    return ret;}double vlen(point3 v){    return sqrt(v.x*v.x + v.y*v.y + v.z*v.z);}point3 pvec(plane3 s){    return xmult(subt(s.a,s.b),subt(s.b,s.c));}point3 pvec(point3 s1,point3 s2,point3 s3){    return xmult(subt(s1,s2),subt(s2,s3));}int dots_onplane(point3 a,point3 b,point3 c,point3 d){    return zero(dmult(pvec(a,b,c),subt(d,a)));}double ptoplane(point3 p,plane3 s){    return fabs(dmult(pvec(s),subt(p,s.a)))/vlen(pvec(s));}point3 p[5]; //a,b,c,dline3 ab,ac,ad,bc,bd,cd;int main(){    while(cin >> p[1].x >> p[1].y >> p[1].z)    {        for(int i = 2 ; i <= 4 ; i ++)            cin >> p[i].x >> p[i].y >> p[i].z;        if(dots_onplane(p[1],p[2],p[3],p[4]))        {            printf("O O O O\n");            continue;        }        plane3 abc(p[1],p[2],p[3]),abd(p[1],p[2],p[4]),acd(p[1],p[3],p[4]),bcd(p[2],p[3],p[4]);        double s = abc.s() + abd.s() + acd.s() + bcd.s();        point3 o;        o.x = (abc.s()*p[4].x + abd.s()*p[3].x + acd.s()*p[2].x + bcd.s()*p[1].x)/s;        o.y = (abc.s()*p[4].y + abd.s()*p[3].y + acd.s()*p[2].y + bcd.s()*p[1].y)/s;        o.z = (abc.s()*p[4].z + abd.s()*p[3].z + acd.s()*p[2].z + bcd.s()*p[1].z)/s;        double r = ptoplane(o,abc);        printf("%.4lf %.4lf %.4lf %.4lf\n",o.x,o.y,o.z,r);    }    return 0;}