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