四面体内接圆

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  

题意：

给出四个点，问四个点组成的四面体的内接圆圆心和半径，如果没有内接圆就输出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,d line3 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; }