ACM几何算法

//    叉乘
//    两个点的距离 

//    点到直线距离
//    返回直线 Ax + By + C =0  的系数 
//    线段
//    圆
//    两个圆的公共面积
//    矩形
//    根据下标返回多边形的边 
//    两个矩形的公共面积 
//    多边形  ,逆时针或顺时针给出x,y 
//    多边形顶点 
//    多边形的边 
//    多边形的周长 
//    判断点是否在线段上 
//    判断两条线断是否相交，端点重合算相交 
//    判断两条线断是否平行
//    判断两条直线断是否相交
//    直线相交的交点 
//    判断是否简单多边形 
//    求多边形面积 
//    判断是否在多边形上 
//    判断是否在多边形内部
//    点阵的凸包，返回一个多边形 

//  最近点对的距离
#include <cmath>
#include <cstdio>
#include <memory>
#include <algorithm>
#include <iostream>
using namespace std;

typedef double TYPE;  //把double 定义为TYPE

#define Abs(x) (((x)>0)?(x):(-(x)))  //用Abs()这个宏定义一个绝对值函数

#define Sgn(x) (((x)<0)?(-1):(1)) //取相反数
#define Max(a,b) (((a)>(b))?(a):(b)) //取大值
#define Min(a,b) (((a)<(b))?(a):(b)) //取小值

#define Epsilon 1e-10 
#define Infinity 1e+10
#define Pi 3.14159265358979323846 //定义几个常量
TYPE Deg2Rad(TYPE deg)
{return (deg * Pi / 180.0);}  //把角度制转化为弧度制

TYPE Rad2Deg(TYPE rad)
{return (rad * 180.0 / Pi);}  //把弧度制转化为角度制

TYPE Sin(TYPE deg)
{return sin(Deg2Rad(deg));} //对弧度制求正弦

TYPE Cos(TYPE deg)
{return cos(Deg2Rad(deg));} //对弧度制求余弦

TYPE ArcSin(TYPE val)
{return Rad2Deg(asin(val));}  //求反正弦

TYPE ArcCos(TYPE val)
{ return Rad2Deg(acos(val));} //求反余弦

TYPE Sqrt(TYPE val)
{ return sqrt(val);}

struct POINT
{
    TYPE x;
    TYPE y;
    TYPE z;
    POINT() : x(0), y(0), z(0) {};
    POINT(TYPE _x_, TYPE _y_, TYPE _z_ = 0) : x(_x_), y(_y_), z(_z_) {};
};

// cross product of (o->a) and (o->b)
// 叉乘 
TYPE Cross(const POINT & a, const POINT & b, const POINT & o)
{return (a.x - o.x) * (b.y - o.y) - (b.x - o.x) * (a.y - o.y);}  //叉积模板

// planar points' distance
//  两个点的距离 
TYPE Distance(const POINT & a, const POINT & b)
{return Sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y) + (a.z - b.z) * (a.z - b.z));}

struct LINE
{
    POINT a;
    POINT b;
    LINE() {};
    LINE(POINT _a_, POINT _b_) : a(_a_), b(_b_) {};
};   //直线由两点决定

//点到直线距离
double PointToLine(POINT p0 ,POINT p1 ,POINT p2 ,POINT &cp)
{ 
     double d = Distance(p1 ,p2);
     double s = Cross(p1 ,p2 ,p0) / d;
     cp.x = p0.x + s*( p2.y-p1.y) / d;
     cp.y = p0.y - s*( p2.x-p1.x) / d;
     return Abs(s);
}

// 返回直线 Ax + By + C =0  的系数 
void Coefficient(const LINE & L, TYPE & A, TYPE & B, TYPE & C)
{
    A = L.b.y - L.a.y;
    B = L.a.x - L.b.x;
    C = L.b.x * L.a.y - L.a.x * L.b.y;
}

void Coefficient(const POINT & p,const TYPE a,TYPE & A,TYPE & B,TYPE & C)
{
    A = Cos(a);
    B = Sin(a);
    C = - (p.y * B + p.x * A);
}

// 线段
struct SEG
{
    POINT a;
    POINT b;
    SEG() {};
    SEG(POINT _a_, POINT _b_):a(_a_),b(_b_) {};
};

// 圆 
struct CIRCLE
{
    TYPE x;
    TYPE y;
    TYPE r;
    CIRCLE() {}
    CIRCLE(TYPE _x_, TYPE _y_, TYPE _r_) : x(_x_), y(_y_), r(_r_) {}
};  //圆由圆心和半径组成

POINT Center(const CIRCLE & circle)
{ return POINT(circle.x, circle.y);}  //返回的是一个POINT类型的对象，他这里没有直接设出一个具体的对象，而是直接用类名来实现，也是可以的，而且他这里是把对象传递进来，返回的是他的圆心

 

//圆的面积
TYPE Area(const CIRCLE & circle)
{ return Pi * circle.r * circle.r;}

//两个圆的公共面积 
TYPE CommonArea(const CIRCLE & A, const CIRCLE & B)
{
    TYPE area = 0.0;

    const CIRCLE & M = (A.r > B.r) ? A : B;
    const CIRCLE & N = (A.r > B.r) ? B : A;  //按照圆半径的大小来把圆区分开来

    TYPE D = Distance(Center(M), Center(N));

    if ((D < M.r + N.r) && (D > M.r - N.r)) //判断出两个圆之间的关系是相交的
    {
        TYPE cosM = (M.r * M.r + D * D - N.r * N.r) / (2.0 * M.r * D);
        TYPE cosN = (N.r * N.r + D * D - M.r * M.r) / (2.0 * N.r * D);

        TYPE alpha = 2.0 * ArcCos(cosM);
        TYPE beta  = 2.0 * ArcCos(cosN);

        TYPE TM = 0.5 * M.r * M.r * Sin(alpha);
        TYPE TN = 0.5 * N.r * N.r * Sin(beta);

        TYPE FM = (alpha / 360.0) * Area(M);
        TYPE FN = (beta / 360.0) * Area(N);

        area = FM + FN - TM - TN;
  }//相交圆的公共部分的面积就是通过一个扇形的面积减去三角形积 ，然后把这两个部分相加得到

   else if (D <= M.r - N.r) //如果两圆是内含的关系，那么公共圆的面积就是一个圆的面积
    {
        area = Area(N);
    }
    return area;
}

//    矩形
//    矩形的线段 
//        2
//    --------------- b
//    |          | 
//   3 |          |  1
//   a --------------- 
//         0

struct RECT
{
    POINT a;                                   // 左下点 
    POINT b;                                   // 右上点 
    RECT() {};
    RECT(const POINT & _a_, const POINT & _b_)
    {a = _a_; b = _b_;}
};

//根据下标返回多边形的边 （没有看懂是干什么的，或者说是具体怎么操作的）
SEG Edge(const RECT & rect, int idx) //SEG是线段
{
    SEG edge;
    while (idx < 0) idx += 4;
    switch (idx % 4)
    {
    case 0:
        edge.a = rect.a;
        edge.b = POINT(rect.b.x, rect.a.y);
        break;
    case 1:
        edge.a = POINT(rect.b.x, rect.a.y);
        edge.b = rect.b;
        break;
    case 2:
        edge.a = rect.b;
        edge.b = POINT(rect.a.x, rect.b.y);
        break;
    case 3:
        edge.a = POINT(rect.a.x, rect.b.y);
        edge.b = rect.a;
        break;
    default:
        break;
    }
    return edge;
}

TYPE Area(const RECT & rect)
{return (rect.b.x - rect.a.x) * (rect.b.y - rect.a.y);}

//  两个矩形的公共面积 
TYPE CommonArea(const RECT & A, const RECT & B)
{
    TYPE area = 0.0;

    POINT LL(Max(A.a.x, B.a.x), Max(A.a.y, B.a.y)); 
    POINT UR(Min(A.b.x, B.b.x), Min(A.b.y, B.b.y)); //这里很妙，可以直接把相交部分的左下方的点与右上方的点的坐标求出来

    if ((LL.x <= UR.x) && (LL.y <= UR.y))  //判断是否两个矩形是否相交
    {
        area = Area(RECT(LL, UR));
    }
    return area;
}

// 多边形  ,逆时针或顺时针给出x,y 
struct POLY
{
    int n;        //n个点 
    TYPE * x;     //x,y为点的指针，首尾必须重合 
    TYPE * y;
    POLY() : n(0), x(NULL), y(NULL) {};
    POLY(int _n_, const TYPE * _x_, const TYPE * _y_)
    {     
        n = _n_;
        x = new TYPE[n + 1];
        memcpy(x, _x_, n*sizeof(TYPE));
        x[n] = _x_[0];

        y = new TYPE[n + 1];
        memcpy(y, _y_, n*sizeof(TYPE));
        y[n] = _y_[0];
    }
};

//多边形顶点 
POINT Vertex(const POLY & poly, int idx)
{
    idx %= poly.n;  //idx应该指的是点的个数
    return POINT(poly.x[idx], poly.y[idx]);  //由于POLY里面的x,y是指针，那么就相当于是关于多边形顶点的坐标的一个数组，用来记录下所有多边形顶点的坐标信息
}

//多边形的边 
SEG Edge(const POLY & poly, int idx)
{
    idx %= poly.n;
    return SEG(POINT(poly.x[idx], poly.y[idx]), 
        POINT(poly.x[idx + 1], poly.y[idx + 1]));
} 

//多边形的周长 
TYPE Perimeter(const POLY & poly)
{
    TYPE p = 0.0;
    for (int i = 0; i < poly.n; i++)
        p = p + Distance(Vertex(poly, i), Vertex(poly, i + 1));
    return p;
}

bool IsEqual(TYPE a, TYPE b)
{return (Abs(a - b) < Epsilon);} //基础的用于判断两个数是否相等

bool IsEqual(const POINT & a, const POINT & b)
{return (IsEqual(a.x, b.x) && IsEqual(a.y, b.y));}  //重载一下用于判断两个点是否相等

bool IsEqual(const LINE & A, const LINE & B)
{
    TYPE A1, B1, C1;
    TYPE A2, B2, C2;

    Coefficient(A, A1, B1, C1);  //把直线A的系数返回给A1,B1,C1
    Coefficient(B, A2, B2, C2);

    return IsEqual(A1 * B2, A2 * B1) && 
        IsEqual(A1 * C2, A2 * C1) &&
        IsEqual(B1 * C2, B2 * C1);
}  //判断两条直线是否相等

// 判断点是否在线段上 
bool IsOnSeg(const SEG & seg, const POINT & p)
{
    return (IsEqual(p, seg.a) || IsEqual(p, seg.b)) ||
        (((p.x - seg.a.x) * (p.x - seg.b.x) < 0 || 
        (p.y - seg.a.y) * (p.y - seg.b.y) < 0) &&
        (IsEqual(Cross(seg.b, p, seg.a), 0)));
}  //前面两个IsEqual是用来说明要判断的点是否就是直线的一个端点，中间两个用于判断这个点是否会在线的延长线上，最后一个用于判断三点是否共线。只有这样才能有效的保证点落在线段上

//判断两条线断是否相交，端点重合算相交 
bool IsIntersect(const SEG & u, const SEG & v)
{
    return (Cross(v.a, u.b, u.a) * Cross(u.b, v.b, u.a) >= 0) &&
        (Cross(u.a, v.b, v.a) * Cross(v.b, u.b, v.a) >= 0) &&
        (Max(u.a.x, u.b.x) >= Min(v.a.x, v.b.x)) && 
        (Max(v.a.x, v.b.x) >= Min(u.a.x, u.b.x)) && 
        (Max(u.a.y, u.b.y) >= Min(v.a.y, v.b.y)) && 
        (Max(v.a.y, v.b.y) >= Min(u.a.y, u.b.y));
} //后面的四个比较用于限定线段相交时的一些点坐标之间的关系，只有这样限定后才可以保证线段相交，否则如果没有这四个条件的话，只能保证的是直线的相交。。。

//判断两条线断是否平行
bool IsParallel(const LINE & A, const LINE & B)
{
    TYPE A1, B1, C1;
    TYPE A2, B2, C2;

    Coefficient(A, A1, B1, C1);
    Coefficient(B, A2, B2, C2);

    return (A1 * B2 == A2 * B1) && 
        ((A1 * C2 != A2 * C1) || (B1 * C2 != B2 * C1));
}  //只要系数成比例就行

//判断两条直线断是否相交
bool IsIntersect(const LINE & A, const LINE & B)
{return !IsParallel(A, B);} //不平行就相交

//直线相交的交点 
POINT Intersection(const LINE & A, const LINE & B)
{
    TYPE A1, B1, C1;
    TYPE A2, B2, C2;

    Coefficient(A, A1, B1, C1);
    Coefficient(B, A2, B2, C2);

    POINT I(0, 0);
    I.x = - (B2 * C1 - B1 * C2) / (A1 * B2 - A2 * B1);
    I.y =   (A2 * C1 - A1 * C2) / (A1 * B2 - A2 * B1);
    return I;
}

//判断矩形是否在圆内
bool IsInCircle(const CIRCLE & circle, const RECT & rect)
{
    return (circle.x - circle.r >= rect.a.x) &&
        (circle.x + circle.r <= rect.b.x) &&
        (circle.y - circle.r >= rect.a.y) &&
        (circle.y + circle.r <= rect.b.y);
}

//判断是否简单多边形 
bool IsSimple(const POLY & poly)
{
    if (poly.n < 3)
        return false;
    SEG L1, L2;
    for (int i = 0; i < poly.n - 1; i++)
    {
        L1 = Edge(poly, i); //用Edge取出多边形的边
        for (int j = i + 1; j < poly.n; j++) //用二重循环来对多边形上的边进行一一判断
        {
            L2 = Edge(poly, j);
            if (j == i + 1)  //相邻的两条边进行比较
            {
                if (IsOnSeg(L1, L2.b) || IsOnSeg(L2, L1.a))  return false;  //如果第二条直线的右边的点在第一条直线上或者第一条直线的坐边的点缀第二条直线上，说明这个多边形是凹的，不是简单的                  
            }
            else if (j == poly.n - i - 1)  //同上面的一样，不过这里取出的是左边的相邻直线
            {
                if (IsOnSeg(L1, L2.a) || IsOnSeg(L2, L1.b))  return false;                    
            }
            else
            {
                if (IsIntersect(L1, L2))  return false; //不相邻的直线但是相交，那么说明非简单
            }
        } // for j
    } // for i
    return true;
}

//求多边形面积 
TYPE Area(const POLY & poly)
{
    if (poly.n < 3) return TYPE(0); //如果边数小于3说明非多边形，没有面积可言
    double s = poly.y[0] * (poly.x[poly.n - 1] - poly.x[1]);
    for (int i = 1; i < poly.n; i++)
    {
        s += poly.y[i] * (poly.x[i - 1] - poly.x[(i + 1) % poly.n]);
    }
    return s/2;
}  //把多边形分割成三角形的和，不过这里的面积计算的方法没有看懂

//判断点是否在多边形上 
bool IsOnPoly(const POLY & poly, const POINT & p)
{
    for (int i = 0; i < poly.n; i++)
    {
        if (IsOnSeg(Edge(poly, i), p))  return true;
    }
    return false;
}

//判断点是否在多边形内部 
bool IsInPoly(const POLY & poly, const POINT & p) 
{
    SEG L(p, POINT(Infinity, p.y)); //Infinity==1e+10，L这条直线相当于是过p点且平行于x轴的直线
    int count = 0;
    for (int i = 0; i < poly.n; i++)
    {
        SEG S = Edge(poly, i);
        if (IsOnSeg(S, p))
        {
            return false;                        //如果想让 在poly上则返回 true,则改为true 
        }//点在多边形的边上，非内部，返回false
        if (!IsEqual(S.a.y, S.b.y))
        {
            POINT & q = (S.a.y > S.b.y)?(S.a):(S.b);  //q取这条多边形边上y坐标大的点
            if (IsOnSeg(L, q))
            {
                ++count;
            }
            else if (!IsOnSeg(L, S.a) && !IsOnSeg(L, S.b) && IsIntersect(S, L))
            {
                ++count;
            }
        }
    }
    return (count % 2 != 0);
}//这里用count这个计数器的奇偶来判断点是否在多边形内没有看懂

// 点阵的凸包，返回一个多边形 ，Graham-scan算法
POLY ConvexHull(const POINT * set, int n)            // 不适用于点少于三个的情况
{
    POINT * points = new POINT[n];
    memcpy(points, set, n * sizeof(POINT));

    TYPE * X = new TYPE[n];
    TYPE * Y = new TYPE[n];

    int i, j, k = 0, top = 2;
    for(i = 1; i < n; i++)
    {
        if ((points[i].y < points[k].y) ||
            ((points[i].y == points[k].y) &&
            (points[i].x < points[k].x)))
        {
            k = i;
        }
    } //找到p0点，一般取y坐标最小，如果y相同则取x坐标最小，即最左边的点

    std::swap(points[0], points[k]); 

    for (i = 1; i < n - 1; i++)
    {
        k = i;
        for (j = i + 1; j < n; j++)
        {
            if ((Cross(points[j], points[k], points[0]) > 0) ||
                ((Cross(points[j], points[k], points[0]) == 0) &&
                (Distance(points[0], points[j]) < Distance(points[0], points[k]))))
            {
                k = j;
            }
        }
        std::swap(points[i], points[k]);
    }

    X[0] = points[0].x; Y[0] = points[0].y;
    X[1] = points[1].x; Y[1] = points[1].y;
    X[2] = points[2].x; Y[2] = points[2].y;

    for (i = 3; i < n; i++)
    {
        while (Cross(points[i], POINT(X[top], Y[top]), 
            POINT(X[top - 1], Y[top - 1])) >= 0)
        {
            top--;
        }
        ++top;
        X[top] = points[i].x;
        Y[top] = points[i].y;
    }

    delete [] points;
    POLY poly(++top, X, Y);
    delete [] X;
    delete [] Y;
    return poly;
}

 

//最近点对的距离, Written By PrincessSnow
#define MAXN 100000
POINT pt[MAXN];

bool cmp(POINT n1, POINT n2)
{return (n1.x<n2.x || n1.x==n2.x && n1.y<n2.y);}

double Get(double dis, int mid, int start, int end)
{
    int s=mid, e=mid, i, j;
    double t;
    while(s > start && pt[mid].x - pt[s].x <= dis)     s--;
    while(e < end && pt[e].x - pt[mid].x <= dis)     e++;
    for(i=s; i <= e; i++)
        for(j=i+1; j <= e && j <= i+7; j++) {
            t = Distance(pt[i], pt[j]);
            if(t < dis)     dis=t;
        }
     return dis;
}

double ClosestPairDistance(int start, int end)
{
    int m = end-start+1, mid, i;
    double t1, t2, dis=-1, t;
    if(m <= 3) {
        for(i=start; i < end; i++) {
            t = Distance(pt[i] , pt[i+1]);
            if(t < dis || dis == -1)     dis = t;
        }
        t = Distance(pt[start] , pt[end]);
        if(t < dis) dis=t;
        return dis;
    }
     
    if(m%2 == 0)     mid = start + m/2 - 1;
    else               mid = start + m/2;
    if(m%2 == 0) {
        t1 = ClosestPairDistance(start, mid);
        t2 = ClosestPairDistance(mid+1, end);
    }
    else {
        t1 = ClosestPairDistance(start, mid);
        t2 = ClosestPairDistance(mid+1, end);
    }
    if(t1 < t2)     dis = t1;
    else          dis = t2;
    dis = Get(dis, mid, start, end);
    return dis;
}