算法设计与优化编程 第十一讲 tetrahedron

tetrahedron

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

思路：

首先用叉积判断是否四点共面，若四点共面，则不能形成四面体。随后，利用四面体体积公式，海伦公式（求三角形面积），四面体内心公式，求出四面体的内切圆半径和内心。

设四个点是A、B、C、D，那么内切球的圆心可以表示为

||||代表欧几里德的范式，本质是求距离。
注意： 在||(b-a)(c-a)||d中，b-a是求两个点b和a的差，返回值ba点={b.x-a.x,
b.y-a.y,b.z-a.z}。同理，c-a是求两个点c和a的差，返回值ca点={c.x-a.x,
c.y-a.y,c.z-a.z}。(b-a)(c-a)是两个向量ba点和ca点的叉积，返回值baca点=
{(b.y-a.y)(c.z-a.z)-(c.y-a.y)(b.z-a.z),
(b.z-a.z)(c.x-a.x)-(b.x-a.x)(c.z-a.z),
(b.x-a.x)(c.y-a.y)-(b.y-a.y)(c.x-a.x)}
||(b-a)(c-a)||是求向量baca点的长度，返回值=，
= sqrt（((b.y-a.y)(c.z-a.z)-(c.y-a.y)(b.z-a.z))^2
+((b.z-a.z)(c.x-a.x)-(b.x-a.x)(c.z-a.z))^2
+( (b.x-a.x)(c.y-a.y)-(b.y-a.y)(c.x-a.x))^2），结果是一个常数e
||(b-a)(c-a)||d是求向量d的倍数乘积，返回值={e*d.x, e*d.y, e*d.z}
还差的函数：常数乘向量，向量加法。

半径r可以表示为

程序代码：

#include <iostream>#include <cstdio>#include <cstdlib>#include <cmath>#include <map>#include <vector>#include <queue>#include <cstring>#include <string>#include <algorithm>using namespace std;typedef  long  ll;typedef unsigned long  ull;#define MM(a,b) memset(a,b,sizeof(a));#define inf 0x7f7f7f7f#define FOR(i,n) for(int i=1;i<=n;i++)#define CT continue;#define PF printf#define SC scanfconst int mod=1000000007;const int N=1e6+100;int n,m,c[N],pre[N],sum[N],a[N],ans[N];struct Point{  ll x,y,z;  void read()  {      scanf("%lld%lld%lld",&x,&y,&z);  }}p[6];/*两点相减，在||（b-a）*(c-a)||*d中，b-a是求两个点b和a的差，返回值ba点={b.x-a.x,b.y-a.y,b.z-a.z}Point operator-(Point a,Point b){    return  (Point){b.x-a.x,b.y-a.y,b.z-a.z};}*/Point operator-(Point a,Point b){     Point m;    m.x=b.x-a.x;    m.y=b.y-a.y;    m.z=b.z-a.z;    return  m;}/*两点叉积，（b-a）*(c-a)，ba，ca的叉积，{a.y*b.z-b.y*a.z,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x}Point cross(Point a,Point b){    return (Point){a.y*b.z-b.y*a.z,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x};}*/Point cross(Point a,Point b){     Point n;    n.x=a.y*b.z-b.y*a.z;    n.y=a.z*b.x-a.x*b.z;    n.z=a.x*b.y-a.y*b.x;    return  n;} //点到原点的距离，||（b-a）*(c-a)||是求向量baca点的长度，结果是一个常数，||·||欧几里得的范式，本质是求距离double dis(Point a){   return sqrt(a.x*a.x+a.y*a.y+a.z*a.z);}//两个向量的点积double dot(Point a,Point b){    return a.x*b.x+a.y*b.y+a.z*b.z;} //判断是否共面，共面dot(m,c-a)为0double pointtoface(Point c,Point a,Point b,Point d){    Point m=cross(b-a,d-a);    return dot(m,c-a)/dis(m);}//三维几何中点到面的距离利用 向量a*向量b=|a|*|b|cos(ang)int main(){    long double s[5];    //输入四个点的坐标    while(~scanf("%lld%lld%lld",&p[1].x,&p[1].y,&p[1].z))    {        p[2].read();        p[3].read();        p[4].read();        //判断是否共面，（b-a）*(c-a)与d求点积，点积为0则共面        if(dot(cross(p[2]-p[1],p[3]-p[1]),p[4])==0)        {printf("O O O O\n");        CT;}        double ts=0;        //求面积        //(c-b)*(d-b)        s[1]=dis(cross(p[3]-p[2],p[4]-p[2]))/2;        //(c-a)*(d-a)        s[2]=dis(cross(p[3]-p[1],p[4]-p[1]))/2;        //(b-a)*(d-a)        s[3]=dis(cross(p[2]-p[1],p[4]-p[1]))/2;        //(c-a)*(b-a)        s[4]=dis(cross(p[3]-p[1],p[2]-p[1]))/2;        printf("%f %f %f %f \n",s[1],s[2],s[3],s[4]);        //求面积和        for(int i=1;i<=4;i++)        ts+=s[i];        //求高        double h=pointtoface(p[3],p[1],p[2],p[4]);        double r=fabs(s[3]/ts*h);        //求坐标        double x=(s[1]*p[1].x+s[2]*p[2].x+s[3]*p[3].x+s[4]*p[4].x)/ts;        double y=(s[1]*p[1].y+s[2]*p[2].y+s[3]*p[3].y+s[4]*p[4].y)/ts;        double z=(s[1]*p[1].z+s[2]*p[2].z+s[3]*p[3].z+s[4]*p[4].z)/ts;        printf("%.4f %.4f %.4f %.4f %.4f\n",ts,x,y,z,r);    }    return 0;}