HDU 5733 (三维几何)

题目链接：点击这里

题意：给出四个点， 判断是否构成四面体， 并且求出内切球球心和半径。

判断可以三点定平面然后判断第四个点在不在平面上。

四面体的体积比较好求，直接求一个面的面积然后点到面的距离乱搞。然后半径也好求了。求球心可以先把两个平面沿着法向量的方向推进半径的长度求出此时两个新的平面的交线，然后再用第三个面也按照这个方法推进半径的长度，求出和交线的交点，这个交点就是球心。

这么一通下来接下来只需要堆板子了。

#include <cstdio>
#include <cmath>
#include <algorithm>
#include <iostream>
#include <vector>
using namespace std;
#define maxn 100005

const double eps = 1e-8;
const double INF = 1e20;
const double pi = acos (-1.0);

int dcmp (double x) {
    if (fabs (x) < eps) return 0;
    return (x < 0 ? -1 : 1);
}
struct Point {
    double x, y, z;
    Point (double _x = 0, double _y = 0, double _z = 0) : x(_x), y(_y), z(_z) {}
    void input () {
        scanf ("%lf%lf%lf", &x, &y, &z);
    }
    void output () {
        printf ("%.2f %.2f %.2f\n", x, y, z);
    }
    bool operator == (const Point &b) const {
        return dcmp (x-b.x) == 0 && dcmp (y-b.y) == 0 && dcmp (z-b.z) == 0;
    }
    double len2 () {
        return x*x+y*y+z*z;
    }
    double len () {
        return sqrt (x*x + y*y + z*z);
    }
    Point operator - (const Point &b) const {
        return Point (x-b.x, y-b.y, z-b.z);
    }
    Point operator + (const Point &b) const {
        return Point (x+b.x, y+b.y, z+b.z);
    }
    Point operator * (const double &k) const {
        return Point (x*k, y*k, z*k);
    }
    Point operator / (const double &k) const {
        return Point (x/k, y/k, z/k);
    }
    double operator * (const Point &b) const {
        return x*b.x+y*b.y+z*b.z;
    }
    Point operator ^ (const Point &b) const {
        return Point (y*b.z-z*b.y, z*b.x-x*b.z, x*b.y-y*b.x);
    }
    double rad (Point a, Point b) {//两向量的夹角
        Point p = (*this);
        return acos (((a-p)*(b-p)) / (a.distance (p)*b.distance (p)));
    }
    Point trunc (double r) {
        double l = len ();
        if (!dcmp (l)) return *this;
        r /= l;
        return Point (x*r, y*r, z*r);
    }
    double distance (Point p) {
        return Point (x-p.x, y-p.y, z-p.z).len ();
    }
};

struct Line {
    Point s, e;
    Line () {}
    Line (Point _s, Point _e) {
        s = _s, e = _e;
    }
    void input () {
        s.input ();
        e.input ();
    }
    void output () {
        cout << "line:" << endl;
        s.output ();
        e.output ();
    }
    double len () {
        return (e-s).len ();
    }
    double point_to_line (Point p) {//点到直线的距离
        return ((e-s)^(p-s)).len () / s.distance (e);
    }
    Point prog (Point p) {//点在直线上的投影
        return s+(((e-s)*((e-s)*(p-s))) / ((e-s).len2 ()));
    }
    bool Point_on_line (Point p) {//点在直线上
        return dcmp (((s-p)^(e-p)).len ()) == 0 && dcmp ((s-p)*(e-p)) == 0;
    }
};

struct Plane {
    Point a, b, c, o;
    Plane () {}
    Plane (Point _a, Point _b, Point _c) {
        a = _a, b = _b, c = _c;
        o = pvec ();
    }
    Point pvec () {//法向量
        return (b-a)^(c-a);
    }
    Plane go (Point p, double x) {//平面沿p方向前进x距离
        o = o.trunc (x);
        double ang = a.rad (a+o, p);
        if (ang >= pi/2) o = o*(-1);
        return Plane (a+o, b+o, c+o);
    }
    bool point_on_plane (Point p) {//点在平面上
        return dcmp ((p-a)*o) == 0;
    }
    double plane_angle (Plane f) {//两平面夹角
        return acos (o*f.o)/(o.len ()*f.o.len ());
    }
    int line_cross_plane (Line u, Point &p) {//直线和平面的交点
        double x = o*(u.e-a);
        double y = o*(u.s-a);
        double d = x-y;
        if (dcmp (d) == 0) return 0;
        p = ((u.s*x) - (u.e*y))/d;
        return 1;
    }
    double area () {//三角形面积
        double x1 = (a-b).len (), x2 = (a-c).len (), x3 = (b-c).len ();
        double p = (x1+x2+x3)/2;
        return sqrt (p*(p-x1)*(p-x2)*(p-x3));
    }
    Point point_to_plane (Point p) {
        Line u = Line (p, p+o);
        line_cross_plane (u, p);
        return p;
    }
    int plane_cross_plane (Plane f, Line &u) {//平面交线
        Point oo = o^f.o;
        Point v = o^oo;
        double d = fabs (f.o*v);
        if (dcmp (d) == 0)
            return 0;
        Point q = a-(v*(f.o*(f.a-a))/d);
        u = Line (q, q+oo);
        return 1;
    }
};

Point A, B, C, D;
Plane p1, p2, p3, p4;

int main () {
    while (scanf ("%lf%lf%lf", &A.x, &A.y, &A.z) == 3) {
        B.input (), C.input (), D.input ();
        p1 = Plane (A, B, C);
        p2 = Plane (A, B, D);
        p3 = Plane (A, C, D);
        p4 = Plane (B, C, D);
        if (p1.point_on_plane (D)) {
            printf ("O O O O\n");
            continue;
        }
        Point tmp = p1.point_to_plane (D);
        double dis = (D-tmp).len ();
        double V = dis*p1.area ();
        double r = V/(p1.area () + p2.area () + p3.area () + p4.area ());
        p1 = p1.go (D, r);
        p2 = p2.go (C, r);
        p3 = p3.go (B, r);
        Line intersection;
        p2.plane_cross_plane (p3, intersection);
        Point ans;
        p1.line_cross_plane (intersection, ans);
        printf ("%.4f %.4f %.4f %.4f\n", ans.x+eps, ans.y+eps, ans.z+eps, r);
    }
    return 0;
}