hdu 5733 （计算几何）

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
  
  

 

参考队长的题解写的。

题意：给定三维空间中四个点，如果能构成一个四面体，则输出这个四面体的内切圆的球心和半径，否则输出四个0;

直接给出公式：

设这四个点分别为ABCD，这四个点对面的面积为SA,SB,SC,SD。其球心为

X=(Xa*SA+Xb*SB+Xc*SC+Xd*SD)/(SA+SB+SC+SD)。

Y=(Ya*SA+Yb*SB+Yc*SC+Yd*SD)/(SA+SB+SC+SD)。

Z=(Za*SA+Zb*SB+Zc*SC+Zd*SD)/(SA+SB+SC+SD)。

半径为

R=3V/(SA+SB+SC+SD)

V是这个四面体，可以用任意一个点引出的三个向量求混合积/6求得。

#include<iostream>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
const double eps = 1e-8;
int cmp(double x)
{
    if (fabs(x) < eps) return 0;
    if (x < 0) return -1;
    return 1;
}
struct Point
{
    double x, y, z;
    Point()
    {
        x = y = z = 0;
    }
    Point(double x, double y, double z) :x(x), y(y), z(z) {}
    int read()
    {
        return scanf("%lf %lf %lf", &x, &y, &z);
    }
    Point operator - (const Point &b) const
    {
        return Point(x - b.x, y - b.y, z - b.z);
    }
    double norm()
    {
        return sqrt(x*x + y*y + z*z);
    }
};

double dot(Point a, Point b)
{
    return a.x*b.x + a.y*b.y + a.z*b.z;
}
Point det(Point a, Point b)
{
    return Point(a.y*b.z - a.z*b.y, a.z*b.x - a.x*b.z, a.x*b.y - a.y*b.x);
}

int main()
{
    Point A[4];
    while (~A[0].read())
    {
        for (int i = 1; i <= 3; i++) A[i].read();
        Point v[3];
        for (int i = 0; i < 3; i++)
            v[i] = A[i + 1] - A[0];
        double V = dot(v[2], det(v[0], v[1]));
        if (!cmp(V))
        {
            puts("O O O O");
            continue;
        }
        V /= 6;
        V = fabs(V);
        double s[4];
        for (int i = 0; i < 4; i++)
        {
            int tmp[4];
            int cnttmp = 0;
            for (int j = 0; j < 4; j++)
                if (j != i) tmp[cnttmp++] = j;
            s[i] = fabs(det(A[tmp[1]] - A[tmp[0]], A[tmp[2]] - A[tmp[0]]).norm()) / 2;
        }
        Point ans;
        for (int i = 0; i < 4; i++)
        {
            ans.x += s[i] * A[i].x;
            ans.y += s[i] * A[i].y;
            ans.z += s[i] * A[i].z;
        }
        ans.x /= s[0] + s[1] + s[2] + s[3];
        ans.y /= s[0] + s[1] + s[2] + s[3];
        ans.z /= s[0] + s[1] + s[2] + s[3];
        double R =3* V / (s[0] + s[1] + s[2] + s[3]);
        printf("%.4lf %.4lf %.4lf %.4lf\n", ans.x, ans.y, ans.z, R);
    }
    return 0;
}