算法设计与优化编程 第十一讲 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 scanf
const 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)为0

double 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;
}