链接:戳这里
tetrahedron
Time Limit: 2000/1000 MS (Java/Others) Memory Limit: 65536/65536 K (Java/Others)
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
思路:
s1,s2,s3,s4表示四面体的四个面的面积,V表示四面体的体积
R=3V/(s1+s2+s3+s4)
球心计算公式:
代码:
/*
S1 -- S4 是四面体的四个面积,内切球的半径为高,这样的4个体积加起来就是总体积
内切球半径的求法 r = 3V / (S1+S2+S3+S4)
算四面体内切球的球心
http://www.docin.com/p-504197705.html?qq-pf-to=pcqq.c2c
*/
#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<string>
#include<vector>
#include <ctime>
#include<queue>
#include<set>
#include<map>
#include<stack>
#include<iomanip>
#include<cmath>
#define mst(ss,b) memset((ss),(b),sizeof(ss))
#define maxn 0x3f3f3f3f
#define MAX 1000100
///#pragma comment(linker, "/STACK:102400000,102400000")
typedef long long ll;
typedef unsigned long long ull;
#define INF (1ll<<60)-1
using namespace std;
const double eps=1e-8;
#define zero(x) (((x)>0?(x):-(x))<eps)
struct point3{
double x,y,z;
point3(double x=0,double y=0,double z=0):x(x),y(y),z(z){}
};
typedef point3 vec3;
vec3 operator + (vec3 a,vec3 b){
return vec3(a.x+b.x,a.y+b.y,a.z+b.z);
}
vec3 operator - (vec3 a,vec3 b){
return vec3(a.x-b.x,a.y-b.y,a.z-b.z);
}
vec3 operator * (vec3 a,double p){
return vec3(a.x*p,a.y*p,a.z*p);
}
vec3 operator / (vec3 a,double p){
return vec3(a.x/p,a.y/p,a.z/p);
}
point3 operator * (point3 a,point3 b){
return point3(a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x);
}
struct line3{
point3 a,b;
};
struct plane3{
point3 a,b,c;
};
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 volume(point3 a, point3 b, point3 c, point3 d) {
return fabs(dmult( (b - a) * (c - a) , (d - a) ) ) / 6.0;
}
/*****平面法向量*********/
point3 pvec(point3 s1,point3 s2,point3 s3){
return xmult(subt(s1,s2),subt(s2,s3));
}
/*****两点距离***********/
double distance(point3 p1,point3 p2){
return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)+(p1.z-p2.z)*(p1.z-p2.z));
}
/*****向量长度***********/
double vlen(point3 p){
return sqrt(p.x*p.x+p.y*p.y+p.z*p.z);
}
/*****点到平面距离***************/
double ptoplane(point3 p,point3 s1,point3 s2,point3 s3){
return fabs(dmult(pvec(s1,s2,s3),subt(p,s1)))/vlen(pvec(s1 ,s2,s3));
}
/*****判四点共面*********/
int dots_onplane(point3 a,point3 b,point3 c,point3 d){
return zero(dmult(pvec(a,b,c),subt(d,a)));
}
/*****判三点共线*********/
int dots_inline(point3 p1,point3 p2,point3 p3){
return vlen(xmult(subt(p1,p2),subt(p2,p3)))<eps;
}
/*****点到平面的投影*****/
point3 shade_ptoplane(point3 p, point3 a, point3 b, point3 c) {
point3 nor = (b - a) * (c - a);
point3 nor0 = ( nor * dmult(nor, p - a) ) / vlen(nor) / vlen(nor);
return (p - nor0);
}
double Dot(point3 a,point3 b){
return a.x*b.x+a.y*b.y+a.z*b.z;
}
double Length(point3 a){
return sqrt(Dot(a,a));
}
point3 Cross(point3 a,point3 b){
return point3(a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x);
}
double area2(point3 a,point3 b,point3 c){
return Length(Cross(b-a,c-a))*0.5;
}
point3 A,B,C,D,O;
int main(){
while(scanf("%lf%lf%lf",&A.x,&A.y,&A.z)!=EOF){
scanf("%lf%lf%lf",&B.x,&B.y,&B.z);
scanf("%lf%lf%lf",&C.x,&C.y,&C.z);
scanf("%lf%lf%lf",&D.x,&D.y,&D.z);
if(dots_onplane(A,B,C,D) || dots_inline(A,B,C)
|| dots_inline(A,C,D) || dots_inline(A,B,D) || dots_inline(B,C,D)){
printf("O O O O\n");
continue;
}
double s1=area2(B,C,D);
double s2=area2(A,C,D);
double s3=area2(A,B,D);
double s4=area2(A,B,C);
double R=3.0*volume(A,B,C,D)/(s1+s2+s3+s4);
double x1=(s1*A.x+s2*B.x+s3*C.x+s4*D.x)/(s1+s2+s3+s4);
double y1=(s1*A.y+s2*B.y+s3*C.y+s4*D.y)/(s1+s2+s3+s4);
double z1=(s1*A.z+s2*B.z+s3*C.z+s4*D.z)/(s1+s2+s3+s4);
printf("%.4f %.4f %.4f %.4f\n",x1,y1,z1,R);
}
return 0;
}