3维计算几何模板

#include <iostream>
#include <cmath>
#include <stdio.h>
#include <vector>
#include <string.h>
#include <stdlib.h>
#include <algorithm>
using namespace std;

#define MAX_N 110

/*------------------常量区-------------------*/

const double INF        = 1e10;      // 无穷大
const double EPS        = 1e-5;      // 计算精度
const double PI         = acos(-1.0);// PI
const int PINXING       = 0;         // 平行
const int XIANGJIAO     = 1;         // 相交
const int XIANGLI       = 0;         // 相离
const int GONGXIAN      = 2;         // 共线
const int CHONGDIE      = -1;        // 重叠
const int INSIDE        = 1;         // 点在图形内部
const int OUTSIDE       = 0;         // 点在图形外部
const int BORDER        = 2;         // 点在图形边界

/*-----------------类型定义区----------------*/

struct Point {              // 二维点或矢量
    double x, y;
    //double angle, dis;
    Point() {}
    Point(double x0, double y0): x(x0), y(y0) {}
    void read()
    {
        scanf("%lf%lf",&x,&y);
    }
};

struct Line {               // 二维的直线或线段
    Point p1, p2;
    Line() {}
    Line(Point p10, Point p20): p1(p10), p2(p20) {}
    void read()
    {
        scanf("%lf%lf%lf%lf",&p1.x,&p1.y,&p2.x,&p2.y);
    }
};

struct Rect {              // 用长宽表示矩形的方法 w, h分别表示宽度和高度
    double w, h;
    Rect() {}
    Rect(double _w,double _h) : w(_w),h(_h) {}
};
struct Rect_2 {             // 表示矩形，左下角坐标是(xl, yl)，右上角坐标是(xh, yh)
    double xl, yl, xh, yh;
    Rect_2() {}
    Rect_2(double _xl,double _yl,double _xh,double _yh) : xl(_xl),yl(_yl),xh(_xh),yh(_yh) {}
};
struct Circle {            //圆
    Point c;
    double r;
    Circle() {}
    Circle(Point _c,double _r) :c(_c),r(_r) {}
};

typedef vector<Point> Polygon;      // 二维多边形
typedef vector<Point> Points;       // 二维点集

/*-------------------基本函数区---------------------*/

inline double max(double x,double y)
{
    return x > y ? x : y;
}
inline double min(double x, double y)
{
    return x > y ? y : x;
}
inline bool ZERO(double x)              // x == 0
{
    return (fabs(x) < EPS);
}
inline bool ZERO(Point p)               // p == 0
{
    return (ZERO(p.x) && ZERO(p.y));
}

inline bool EQ(double x, double y)      // eqaul, x == y
{
    return (fabs(x - y) < EPS);
}
inline bool NEQ(double x, double y)     // not equal, x != y
{
    return (fabs(x - y) >= EPS);
}
inline bool LT(double x, double y)     // less than, x < y
{
    return ( NEQ(x, y) && (x < y) );
}
inline bool GT(double x, double y)     // greater than, x > y
{
    return ( NEQ(x, y) && (x > y) );
}
inline bool LEQ(double x, double y)     // less equal, x <= y
{
    return ( EQ(x, y) || (x < y) );
}
inline bool GEQ(double x, double y)     // greater equal, x >= y
{
    return ( EQ(x, y) || (x > y) );
}

// 输出浮点数前，防止输出-0.00调用该函数进行修正！
inline double FIX(double x)
{
    return (fabs(x) < EPS) ? 0 : x;
}



/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
//-------------------3D 区域----------------------------//

struct Point3D {            //三维点或矢量
    double x, y, z;
    Point3D() {}
    Point3D(double x0, double y0, double z0): x(x0), y(y0), z(z0) {}
    void read()
    {
        scanf("%lf%lf%lf",&x,&y,&z);
    }
};

struct Line3D {             // 三维的直线或线段
    Point3D p1, p2;
    Line3D() {}
    Line3D(Point3D p10, Point3D p20): p1(p10), p2(p20) {}
    void read()
    {
        scanf("%lf%lf%lf%lf%lf%lf",&p1.x,&p1.y,&p1.z,&p2.x,&p2.y,&p2.z);
    }
};

struct Area3D{
    Point3D p1,p2,p3;
    Area3D(){}
    Area3D(Point3D p10, Point3D p20,Point3D p30): p1(p10), p2(p20), p3(p30){}
    void read()
    {
        p1.read(); p2.read(); p3.read();
        //scanf("%lf%lf%lf%lf%lf%lf%lf%lf%lf",&p1.x,&p1.y,&p1.z,&p2.x,&p2.y,&p2.z,&p3.x,&p3.y,&p3.z);
    }
};

inline bool ZERO(Point3D p)              // p == 0
{
    return (ZERO(p.x) && ZERO(p.y) && ZERO(p.z));
}

//////////////////////////////////////////////////////////////////////////////////////
//三维矢量运算
bool operator==(Point3D p1, Point3D p2)
{
    return ( EQ(p1.x, p2.x) && EQ(p1.y, p2.y) && EQ(p1.z, p2.z) );
}
bool operator<(Point3D p1, Point3D p2)
{
    if (NEQ(p1.x, p2.x)) {
        return (p1.x < p2.x);
    } else if (NEQ(p1.y, p2.y)) {
        return (p1.y < p2.y);
    } else {
        return (p1.z < p2.z);
    }
}
Point3D operator+(Point3D p1, Point3D p2)
{
    return Point3D(p1.x + p2.x, p1.y + p2.y, p1.z + p2.z);
}
Point3D operator-(Point3D p1, Point3D p2)
{
    return Point3D(p1.x - p2.x, p1.y - p2.y, p1.z - p2.z);
}
Point3D operator * (const Point3D& A, double p) { return Point3D(A.x*p, A.y*p, A.z*p); }
Point3D operator / (const Point3D& A, double p) { return Point3D(A.x/p, A.y/p, A.z/p); }

Point3D operator*(Point3D p1, Point3D p2) // 计算叉乘 p1 x p2
{
    return Point3D(p1.y * p2.z - p1.z * p2.y,
                   p1.z * p2.x - p1.x * p2.z,
                   p1.x * p2.y - p1.y * p2.x );
}
double operator&(Point3D p1, Point3D p2) { // 计算点积 p1·p2
    return (p1.x * p2.x + p1.y * p2.y + p1.z * p2.z);
}
double Norm(Point3D p) // 计算矢量p的模
{
    return sqrt(p.x * p.x + p.y * p.y + p.z * p.z);
}

//取平面法向量
Point3D GetV(Area3D s){
    return (s.p1-s.p2)*(s.p2-s.p3);
}

//判三点共线
int PointOnLine(Point3D p1,Point3D p2,Point3D p3){
    return ZERO( (p1-p2)*(p2-p3) );
}


//判四点共面
int PointOnArea(Point3D a,Point3D b,Point3D c,Point3D d){
    return ZERO( GetV(Area3D(a, b, c))&(d-a) );
}

//求三维空间中两点间的距离
double Dis(Point3D p1, Point3D 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 Dis(Point3D p, Line3D L)
{
    if(L.p1==L.p2) return Dis(p, L.p1);
    return Norm((p - L.p1) * (L.p2 - L.p1)) / Norm(L.p2 - L.p1);
}

bool OnLine(Point3D p, Line3D L) // 判断三维空间中点p是否在直线L上
{
    if(L.p1==L.p2 && p==L.p1) return true;//共点时，返回true
    return ZERO( (p - L.p1) * (L.p2 - L.p1) );
}
bool OnLineSeg(Point3D p, Line3D L) // 判断三维空间中点p是否在线段l上
{
    return ( ZERO((L.p1 - p) * (L.p2 - p)) &&
            EQ( Norm(p - L.p1) + Norm(p - L.p2), Norm(L.p2 - L.p1)) );
}

// 点p到平面Ap-Al的距离。
double Dis(Point3D p, Point3D Ap, Point3D Al) {
    return fabs((p-Ap)&Al)/Norm(Al); // 如果不取绝对值，得到的是有向距离
}

// 点p在平面Ap-Al上的投影。
Point3D PointToArea(Point3D p,Point3D Ap, Point3D Al) {
    Al=Al/(Norm(Al));//把Al变成法向量。
    return p-Al*((p-Ap)&Al);
}
//得到点p到直线L的距离，并返回p到直直线L的最近点rep
double PointToLine(Point3D p,Line3D L,Point3D& rep)
{
    if(L.p1==L.p2)
    {
        rep=L.p1;
        return Dis(p,L.p1);
    }
    Point3D a,b;
    a = L.p2-L.p1;
    b = p-L.p1;
    double dis12 = Dis(L.p1,L.p2);
    double dis = ( Norm(a*b) )/dis12;
    
    double k = (a&b)/(Norm(a)*dis12) ;
    rep = L.p1+(L.p2-L.p1)*k;
    return dis;
}
//求两条直线之间的关系（三维）
//输入：两条不为点的直线
//输出：相交返回XIANGJIAO和交点p,平行返回PINGXING,共线返回GONGXIAN
int LineAndLine(Line3D L1,Line3D L2,Point3D &p)
{
    Point3D px,py;
    px = L1.p1 - L1.p2;
    py = L2.p1 - L2.p2;
    
    if( ZERO(px*py) )//平行或者共线
    {
        if( ZERO( (L2.p1-L1.p1)*py ) ) //共线
        {
            return GONGXIAN;
        }
        return PINXING;
    }
    //判断是否共面
    Point3D tp=(L1.p1-L2.p1)*py;
    if( !ZERO(tp&px) ) return XIANGLI;//XIANGLI与平行相同
    
    p = L1.p1;
    Point3D tp1=(L2.p1-L1.p1)*(L2.p1-L2.p2);
    Point3D tp2=(L1.p2-L1.p1)*(L2.p1-L2.p2);
    double _t = Norm(tp1)/Norm(tp2);
    //tp1和tp2肯定是共线的，如果反向则_t 为负
    if( LT( (tp1&tp2),0 ) ) _t*=-1;
    p.x += (L1.p2.x-L1.p1.x)*_t;
    p.y += (L1.p2.y-L1.p1.y)*_t;
    p.z += (L1.p2.z-L1.p1.z)*_t;
    return XIANGJIAO;
}

//空间两直线最近点对。直线不能平行，直线不能为点.
//ans1为直线a1,b1上的最近点
Point3D ans1,ans2;
double SegSegDistance(Point3D a1, Point3D b1, Point3D a2, Point3D b2)
{
    Point3D v1 = (a1-b1), v2 = (a2-b2);
    Point3D N = v1*v2;
    Point3D ab = (a1-a2);
    double ans = (N&ab) / Norm(N);
    Point3D p1 = a1, p2 = a2;
    Point3D d1 = b1-a1, d2 = b2-a2;
    double t1, t2;
    t1 = ((p2-p1)*d2 )&(d1*d2);
    t2 = ((p2-p1)*d1 )&(d1*d2);
    double dd = Norm( (d1*d2) );
    t1 /= dd*dd;
    t2 /= dd*dd;
    ans1=a1+(b1-a1)*t1;
    ans2=a2+(b2-a2)*t2;
    return fabs(ans);
}

//直线与平面交
int LineAndArea(Line3D l1,Point3D Ap,Point3D Al,Point3D &p)
{
    //输入直线，和平面的点法式
    //第一步，判断直线与平面是否平行。
    if( ZERO((l1.p2-l1.p1)&Al)  ) return 0;//直线与平面平行。
    double _t =( (Ap-l1.p1)&Al ) / ((l1.p1-l1.p2)&Al);
    p = l1.p1+(l1.p1-l1.p2)*_t;
    return 1;
}

void dfs(int x,double &len)
{
    len++;
    dfs(x-1,len);
    dfs(x-2,len);
}

//空间两直线最近点对
//注意：直线不能平行
double LineAndLine(Line3D l1,Line3D l2,Point3D &p1,Point3D &p2)
{
    //先求出法向量
    Point3D v1,v2;
    v1 = l1.p2-l1.p1;
    v2 = l2.p2-l2.p1;
    Point3D vt=v1*v2;
    //然后先把l2投影到 l1所在的平面上
    double len = ((l2.p1-l1.p1)&vt)/Norm(vt);
    double normvt = -len/Norm(vt);
    
    vt.x = vt.x*normvt;
    vt.y = vt.y*normvt;
    vt.z = vt.z*normvt;
    
    Line3D tl2;
    tl2.p1 = l2.p1+vt;
    tl2.p2 = l2.p2+vt;
    
    int sign=LineAndLine(l1, tl2, p1);
    /*
     //测试用
     if(sign!=XIANGJIAO)
     {
     int x=0;
     printf("%lf
",len/x);
     dfs(100000000,len);
     }
     */
    return fabs(len);
}

//已知空间四面体6条边，求体积
double P( double a,double b,double c,double d,double e ){ return a*(b*c-d*e); }
double Get4V(int OA,int OB,int OC,int AB,int CA,int BC)
{
    OA*=OA;OB*=OB;OC*=OC;AB*=AB;CA*=CA;BC*=BC;
    double ans=0;
    ans+=P( OA,OB,OC,(OB+OC-BC)/2.0,(OB+OC-BC)/2.0 );
    ans-=P( (OA+OB-AB)/2.0,(OA+OB-AB)/2.0,OC,(OA+OC-CA)/2.0,(OB+OC-BC)/2.0 );
    ans+=P( (OA+OC-CA)/2.0,(OA+OB-AB)/2.0,(OB+OC-BC)/2.0,OB,(OA+OC-CA)/2.0);
    return sqrt(ans/36);
}

//求两面相交,平行或共面返回PINGXING,否则返回XIANGJIAO和直线rel
int AreaAndArea(Area3D a1,Area3D a2,Line3D &rel)
{
    Point3D va1 = GetV(a1),va2 = GetV(a2);
    Point3D lv = va1*va2;//相交直线的方向向量
    if( ZERO(lv) )//平行
    {
        return PINXING;
    }
    //然后得到某一个交点
    Point3D p1;
    if(LineAndArea(Line3D(a1.p1,a1.p2), a2.p1, va2, p1) == 0)
       if(LineAndArea(Line3D(a1.p1,a1.p3), a2.p1, va2, p1) == 0)
           LineAndArea(Line3D(a1.p2,a1.p3), a2.p1, va2, p1);
    rel.p1 = p1; rel.p2 = p1 + (lv*10);
    return XIANGJIAO;
}

//已知3点坐标，求平面ax+by+cz+d=0;
void GetAreaABCD(Point3D p1,Point3D p2,Point3D p3,double &a,double &b,double &c,double &d)
{
    a = ( (p2.y-p1.y)*(p3.z-p1.z)-(p2.z-p1.z)*(p3.y-p1.y) );
    b = ( (p2.z-p1.z)*(p3.x-p1.x)-(p2.x-p1.x)*(p3.z-p1.z) );
    c = ( (p2.x-p1.x)*(p3.y-p1.y)-(p2.y-p1.y)*(p3.x-p1.x) );
    d = ( 0-(a*p1.x+b*p1.y+c*p1.z) );
}

//////////////////////////////////////////////////////////////////////////////////////


/*---------------------代码区---------------------------*/