计算几何模板 - 全

#include <cmath>
#include <queue>
#include <stack>
#include <vector>
#include <cstdio>
#include <string>
#include <cctype>
#include <climits>
#include <cstring>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <algorithm>
using namespace std;
#define LL long long
#define lson l, m, rt<<1
#define rson m+1, r, rt<<1|1
#define PI 3.1415926535897932626
#define EXIT exit(0);
#define DEBUG puts("Here is a BUG");
#define CLEAR(name, init) memset(name, init, sizeof(name))
const double eps = 1e-6;
const int MAXN = (int)1e9 + 5;


#define Vector Point

int dcmp(double x)
{
    return fabs(x) < eps ? 0 : (x < 0 ? -1 : 1);
}

struct Point
{
    double x, y;
    int flag;//点在哪条边上
    Point(const Point& rhs): x(rhs.x), y(rhs.y) { }                                 //拷贝构造函数
    Point(double x = 0.0, double y = 0.0, int flag = -1): x(x), y(y), flag(flag) {}                           //构造函数

     friend istream& operator >> (istream& in, Point& P) { return in >> P.x >> P.y; }
     friend ostream& operator << (ostream& out, const Point& P) { return out << P.x << ' ' << P.y; }

    friend Vector operator + (const Vector& A, const Vector& B)
    {
        return Vector(A.x + B.x, A.y + B.y);
    }
    friend Vector operator - (const Point& A, const Point& B)
    {
        return Vector(A.x - B.x, A.y - B.y);
    }
    friend Vector operator * (const Vector& A, const double& p)
    {
        return Vector(A.x * p, A.y * p);
    }
    friend Vector operator / (const Vector& A, const double& p)
    {
        return Vector(A.x / p, A.y / p);
    }
    friend bool operator == (const Point& A, const Point& B)
    {
        return dcmp(A.x - B.x) == 0 && dcmp(A.y - B.y) == 0;
    }
    friend bool operator < (const Point& A, const Point& B)
    {
        return A.x < B.x || (A.x == B.x && A.y < B.y);
    }

    void in(void)
    {
        scanf("%lf%lf", &x, &y);
    }
    void out(void)
    {
        printf("%lf %lf", x, y);
    }
};



struct Line
{
    Point P;                                                                        //直线上一点
    Vector dir;                                                                     //方向向量(半平面交中该向量左侧表示相应的半平面)
    double ang;                                                                     //极角，即从x正半轴旋转到向量dir所需要的角（弧度）

    Line() {}                                                                      //构造函数
    Line(const Line& L): P(L.P), dir(L.dir), ang(L.ang) { }
    Line(const Point& P, const Vector& dir): P(P), dir(dir)
    {
        ang = atan2(dir.y, dir.x);
    }

    bool operator < (const Line& L) const                                           //极角排序
    {
        return ang < L.ang;
    }

    Point point(double t)
    {
        return P + dir * t;
    }
};

typedef vector<Point> Polygon;

struct Circle
{
    Point c;                                                                        //圆心
    double r;                                                                       //半径
    Circle() { }
    Circle(const Circle& rhs): c(rhs.c), r(rhs.r) { }
    Circle(const Point& c, const double& r): c(c), r(r) { }

    Point point(double ang)
    {
        return Point(c.x + cos(ang) * r, c.y + sin(ang) * r); //圆心角所对应的
    }
};


double Dot(const Vector& A, const Vector& B)
{
    return A.x * B.x + A.y * B.y; //点积
}
double Length(const Vector& A)
{
    return sqrt(Dot(A, A) + eps);
}
double juli(Point A, Point B)
{
    return sqrt((A.x - B.x) * (A.x - B.x) + (A.y - B.y) * (A.y - B.y));
}

double Angle(const Vector& A, const Vector& B)
{
    return acos(Dot(A, B) / Length(A) / Length(B)); //向量夹角
}
double Cross(const Vector& A, const Vector& B)
{
    return A.x * B.y - A.y * B.x; //叉积
}
double Area(const Point& A, const Point& B, const Point& C)
{
    return fabs(Cross(B - A, C - A));
}

//三边构成三角形的判定
bool check_length(double a, double b, double c)
{
    return dcmp(a + b - c) > 0 && dcmp(fabs(a - b) - c) < 0;
}
bool isTriangle(double a, double b, double c)
{
    return check_length(a, b, c) && check_length(a, c, b) && check_length(b, c, a);
}

//向量绕起点旋转
Vector Rotate(const Vector& A, const double& rad)
{
    return Vector(A.x * cos(rad) - A.y * sin(rad), A.x * sin(rad) + A.y * cos(rad));
}

//向量的单位法线(调用前请确保A 不是零向量)
Vector Normal(const Vector& A)
{
    double len = Length(A);
    return Vector(-A.y / len, A.x / len);
}

//两直线交点(用前确保两直线有唯一交点，当且仅当Cross(A.dir, B.dir)非0)
Point GetLineIntersection(const Line& A, const Line& B)
{
    Vector u = A.P - B.P;
    double t = Cross(B.dir, u) / Cross(A.dir, B.dir);
    return A.P + A.dir * t;
}

//点到直线距离
double DistanceToLine(const Point& P, const Line& L)
{
    Vector v1 = L.dir, v2 = P - L.P;
    return fabs(Cross(v1, v2)) / Length(v1);
}

//点到线段距离
double DistanceToSegment(const Point& P, const Point& A, const Point& B)
{
    if (A == B) return Length(P - A);
    Vector v1 = B - A, v2 = P - A, v3 = P - B;
    if (dcmp(Dot(v1, v2)) < 0) return Length(v2);
    if (dcmp(Dot(v1, v3)) > 0) return Length(v3);
    return fabs(Cross(v1, v2)) / Length(v1);
}

//点在直线上的投影
Point GetLineProjection(const Point& P, const Line& L)
{
    return L.P + L.dir * (Dot(L.dir, P - L.P) / Dot(L.dir, L.dir));
}

//点在线段上的判定
bool isOnSegment(const Point& P, const Point& A, const Point& B)
{
    //若允许点与端点重合，可关闭下面的注释
    //if (P == A || P == B) return true;
    return dcmp(Cross(A - P, B - P)) == 0 && dcmp(Dot(A - P, B - P)) < 0;
}
bool specialOnSegment(const Point& P, const Point& A, const Point& B)
{
    //若允许点与端点重合，可关闭下面的注释
    if (P == A || P == B) return true;
    return dcmp(Cross(A - P, B - P)) == 0 && dcmp(Dot(A - P, B - P)) < 0;
}

//线段相交判定
bool SegmentProperIntersection(const Point& a1, const Point& a2, const Point& b1, const Point& b2)
{
    //若允许在端点处相交，可适当关闭下面的注释
    //if (isOnSegment(a1, b1, b2) || isOnSegment(a2, b1, b2) || isOnSegment(b1, a1, a2) || isOnSegment(b2, a1, a2)) return true;
    double c1 = Cross(a2 - a1, b1 - a1), c2 = Cross(a2 - a1, b2 - a1);
    double c3 = Cross(b2 - b1, a1 - b1), c4 = Cross(b2 - b1, a2 - b1);
    return dcmp(c1) * dcmp(c2) < 0 && dcmp(c3) * dcmp(c4) < 0;
}

//多边形的有向面积
double PolygonArea(Polygon po)
{
    int n = po.size();
    double area = 0.0;
    for (int i = 1; i < n - 1; i++)
    {
        area += Cross(po[i] - po[0], po[i + 1] - po[0]);
    }
    return area * 0.5;
}

//直线和圆的交点
int GetLineCircleIntersection(Line& L, const Circle& C)
{
    double t1, t2;
    double a = L.dir.x, b = L.P.x - C.c.x, c = L.dir.y, d = L.P.y - C.c.y;
    double e = a * a + c * c, f = 2.0 * (a * b + c * d), g = b * b + d * d - C.r * C.r;
    double delta = f * f - 4 * e * g;                                                   //判别式
    if (dcmp(delta) < 0) return 0;                                                      //相离
    if (dcmp(delta) == 0)                                                               //相切
    {
        t1 = t2 = -f / (2 * e);
        sol.push_back(L.point(t1));
        return 1;
    }
    // 相交
    t1 = (-f - sqrt(delta)) / (2.0 * e);
    sol.push_back(L.point(t1));
    t2 = (-f + sqrt(delta)) / (2.0 * e);
    sol.push_back(L.point(t2));
    return 2;
}

//凸包(Andrew算法)
//如果不希望在凸包的边上有输入点，把两个 <= 改成 <
//如果不介意点集被修改，可以改成传递引用
Polygon ConvexHull(vector<Point> p)
{
    //预处理，删除重复点
    sort(p.begin(), p.end());
    p.erase(unique(p.begin(), p.end()), p.end());
    int n = p.size(), m = 0;
    Polygon res(n + 1);
    for (int i = 0; i < n; i++)
    {
        while (m > 1 && Cross(res[m - 1] - res[m - 2], p[i] - res[m - 2]) <= 0) m--;
        res[m++] = p[i];
    }
    int k = m;
    for (int i = n - 2; i >= 0; i--)
    {
        while (m > k && Cross(res[m - 1] - res[m - 2], p[i] - res[m - 2]) <= 0) m--;
        res[m++] = p[i];
    }
    m -= n > 1;
    res.resize(m);
    return res;
}

//点P在有向直线L左边的判定(线上不算)
bool isOnLeft(const Line& L, const Point& P)
{
    return Cross(L.dir, P - L.P) > 0;
}

//半平面交主过程
//如果不介意点集被修改，可以改成传递引用
Polygon HalfPlaneIntersection(vector<Line> L)
{
    int n = L.size();
    int head, rear;                                                                 //双端队列的第一个元素和最后一个元素的下标
    vector<Point> p(n);                                                             //p[i]为q[i]和q[i+1]的交点
    vector<Line> q(n);                                                              //双端队列
    Polygon ans;

    sort(L.begin(), L.end());                                                       //按极角排序
    q[head = rear = 0] = L[0];                                                      //双端队列初始化为只有一个半平面L[0]
    for (int i = 1; i < n; i++)
    {
        while (head < rear && !isOnLeft(L[i], p[rear - 1])) rear--;
        while (head < rear && !isOnLeft(L[i], p[head])) head++;
        q[++rear] = L[i];
        if (fabs(Cross(q[rear].dir, q[rear - 1].dir)) < eps)                        //两向量平行且同向，取内侧的一个
        {
            rear--;
            if (isOnLeft(q[rear], L[i].P)) q[rear] = L[i];
        }
        if (head < rear) p[rear - 1] = GetLineIntersection(q[rear - 1], q[rear]);
    }
    while (head < rear && !isOnLeft(q[head], p[rear - 1])) rear--;                  //删除无用平面
    if (rear - head <= 1) 
        return ans;                                               //空集
    p[rear] = GetLineIntersection(q[rear], q[head]);                                //计算首尾两个半平面的交点

    for (int i = head; i <= rear; i++)                                              //从deque复制到输出中
    {
        ans.push_back(p[i]);
    }
    return ans;
}
//圆的面积交
double calc(double r1, double r2, double d) {
    double alpha = acos((r2*r2+d*d-r1*r1)/2/r2/d)*2;
    double sector = alpha * r2 * r2 / 2;
    double triangle = r2 * r2 * sin(alpha) / 2;
    return sector - triangle;
}
double AreaofCircle(double x1, double y1, double r1, double x2, double y2, double r2)
{
     double d = sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
     if(d <= max(r1, r2) - min(r1, r2) + eps)
     return acos(-1.0) * min(r1, r2) * min(r1, r2);  
     else if (d + eps >= r1 + r2)
         return 0.000;
     else
     return calc(r1, r2, d) + calc(r2, r1, d);  
     //cout << setprecision(20) << AreaofCircle(x1, y1, r1, x2, y2, r2) << endl ;  
}

int main(int argc, char const *argv[])
{

    return 0;
}