计算几何模板(入门级)

尽管才刚入门 然而最近各方面能力都在下降 还是先把这个计算几何入门级模板放上来 弃坑做点其他的

----------------------------------------------------------------------------------------------------------------------

UPD0 16.3.26

UPD1 16.4.3 新增 向量旋转 两圆求交点 余弦定理 以及一些hint 

UPD2 16.4.13 新增 轴对称点

UPD3 16.7.18 新增 三角形面积 凸包上最大三角形 

UPD4 16.7.27 新增 三点确定平面一般方程 点到平面距离

#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;
const double eps = 1e-8, inf = 1e8, pi = acos(-1);
const int N = 110;
struct point
{
    double x, y;
    point(){}
    point(double _x, double _y)
    {
        x = _x;
        y = _y;
    }
    point operator - (const point &p) const
    {
        return point(x - p.x, y - p.y);
    }
    point operator + (const point &p) const
    {
        return point(x + p.x, y + p.y);
    }
    double operator * (const point &p) const
    {
        return x * p.y - y * p.x;
    }
    double operator / (const point &p) const
    {
        return x * p.x + y * p.y;
    }
}a[N], b[N << 1];

struct line
{
    point p1, p2;
    double ang;
    line(){}
    line(point p01, point p02)
    {
        p1 = p01;
        p2 = p02;
    }
    void getang()
    {
        ang = atan2(p2.y - p1.y, p2.x - p1.x);
    }
}c[N], q[N];

struct circle
{
    double x, y, r;
}cc;

struct point3
{
    double x, y, z;
    point3(){}
    point3(double _x, double _y, double _z)
    {
        x = _x;
        y = _y;
        z = _z;
    }
    point3 operator - (const point3 &p) const
    {
        return point3(x - p.x, y - p.y, z - p.z);
    }
    point3 operator + (const point3 &p) const
    {
        return point3(x + p.x, y + p.y, z - p.z);
    }
    point3 operator * (const point3 &p)const
    {
        return point3(y * p.z - z * p.y, z * p.x - x * p.z, x * p.y - y * p.x);
    }
    double operator / (const point3 &p)const
    {
        return x * p.x + y * p.y + z * p.z;
    }
};

int n;

bool cmp(const point &aa, const point &bb)
{
    return aa.x < bb.x || (aa.x == bb.x && aa.y < bb.y);
}

bool cmpl(const line &aa, const line &bb)
{
    if(abs(aa.ang - bb.ang) > eps)
        return aa.ang < bb.ang;
    return (bb.p2 - aa.p1) * (bb.p1 - bb.p2) > 0;
}

int sgn(double x)
//正负判断
{
    if(x > eps)
        return 1;
    if(x < -eps)
        return -1;
    return 0;
}

double getdist2(const point &aa)
{
    return (aa.x * aa.x + aa.y * aa.y);
}

bool online(const point &aa, const point &bb, const point &cc)
//三点共线判断
{
    if(!(min(bb.x, cc.x) <= aa.x && aa.x <= max(bb.x, cc.x)))
        return 0;
    if(!(min(bb.y, cc.y) <= aa.y && aa.y <= max(bb.y, cc.y)))
        return 0;
    return !sgn((bb - aa) * (cc - aa));
}

bool checkline(const point &aa, const point &bb, const point &cc, const point &dd)
//线段是否相交
{
    int c1 = sgn((bb - aa) * (cc - aa)), c2 = sgn((bb - aa) * (dd - aa)), 
    c3 = sgn((dd - cc) * (aa - cc)), c4 = sgn((dd - cc) * (bb - cc));
    if(c1 * c2 < 0 && c3 * c4 < 0)
        return 1;//规范相交
    if(c1 == 0 && c2 == 0)
    {
        if(min(aa.x, bb.x) > max(cc.x, dd.x) || min(cc.x, dd.x) > max(aa.x, bb.x)
        || min(aa.y, bb.y) > max(cc.y, dd.y) || min(cc.y, dd.y) > max(aa.y, bb.y))
            return 0;
        return 1;//共线有重叠
    }
    if(online(cc, aa, bb) || online(dd, aa, bb) 
    || online(aa, cc, dd) || online(bb, cc, dd))
        return 1;//端点相交
    return 0;
}

point getpoint(const point &aa, const point &bb,
 const point &cc, const point &dd)
//两直线交点
{
    double a1, b1, c1, a2, b2, c2;
    point re;
    a1 = aa.y - bb.y;
    b1 = bb.x - aa.x;
    c1 = aa * bb;
    a2 = cc.y - dd.y;
    b2 = dd.x - cc.x;
    c2 = cc * dd;
    //以下为交点横纵坐标
    re.x = (c1 * b2 - c2 * b1) / (a2 * b1 - a1 * b2);
    re.y = (a2 * c1 - a1 * c2) / (a1 * b2 - a2 * b1);
    return re;
}

bool protrusion(point aa[])
//判断多边形凹(凸)
{
    int flag = sgn((aa[1] - aa[0]) * (aa[2] - aa[0]));
    for(int i = 1; i < n; ++i)
        if(sgn((aa[(i + 1) % n] - aa[i]) * (aa[(i + 2) % n] - aa[i])) != flag)
            return 1;
    return 0;
}

bool inpolygon(const point &aa)
//点在凸多边形内(外)
{
    int flag = sgn((a[0] - aa) * (a[1] - aa));
    for(int i = 1; i < n; ++i)
        if(sgn((a[i] - aa) * (a[(i + 1) % n] - aa)) != flag)
            return 0;
    return 1;
}

void transline(line &aa, double dist)
//逆时针方向平移线段
{
    double d = sqrt((aa.p1.x - aa.p2.x) * (aa.p1.x - aa.p2.x) + 
    (aa.p1.y - aa.p2.y) * (aa.p1.y - aa.p2.y));
    point ta;
    ta.x = aa.p1.x + dist / d * (aa.p1.y - aa.p2.y);
    ta.y = aa.p1.y - dist / d * (aa.p1.x - aa.p2.x);
    aa.p2 = ta + aa.p2 - aa.p1;
    aa.p1 = ta;
}

void convexhull(point aa[], point bb[])
//凸包
{
    int len = 0;
    sort(aa, aa + n, cmp);
    bb[len++] = aa[0];//bb数组要开2倍(防止出现直线)
    bb[len++] = aa[1];
    for(int i = 2; i < n ;++i)
    {
        while(len > 1 && (aa[i] - bb[len - 2]) * (bb[len - 1] - bb[len - 2]) > 0)
        //若严格则加上等于
            --len;
        bb[len++] = aa[i];
    }
    int t = len;
    for(int i = n - 2; i >= 0; --i)
    {
        while(len > t && (aa[i] - bb[len - 2]) * (bb[len - 1] - bb[len - 2]) > 0)
        //同上
            --len;
        bb[len++] = aa[i];
    }
    --len;
    n = len;
}

bool checkout(const line &aa, const line &bb, const line &cc)
//检查交点是否在向量顺时针侧
{
    point p = getpoint(aa.p1, aa.p2, bb.p1, bb.p2);
    return (cc.p1 - p) * (cc.p2 - p) < 0;//如果不允许共线或算面积 则此处不取等
}

double halfplane(point aa[], line bb[])
//半平面交 
{
    sort(bb, bb + n, cmpl);
    int n2 = 1;
    for(int i = 1; i < n; ++i)
    {
        if(bb[i].ang - bb[i - 1].ang > eps)
            ++n2;
        bb[n2 - 1]= bb[i];
    }
    n = n2;
    int front = 0, tail = 0;
    q[tail++] = bb[0], q[tail++] = bb[1];
    for(int i = 2; i < n; ++i)
    {
        while(front + 1 < tail && checkout(q[tail - 2], q[tail - 1], bb[i]))
            --tail;
        while(front + 1 < tail && checkout(q[front], q[front + 1], bb[i]))
            ++front;
        q[tail++] = bb[i];
    }
    while(front + 2 < tail && checkout(q[tail - 2], q[tail - 1], q[front]))
        --tail;
    while(front + 2 < tail && checkout(q[front], q[front + 1], q[tail - 1]))
        ++front;
    if(front + 2 >= tail)
        return 0;
    int j = 0;
    for(int i = front; i < tail; ++i, ++j)
    {
        aa[j] = getpoint(q[i].p1, q[i].p2, q[(i != tail - 1 ? i + 1 : front)].p1,
         q[(i != tail - 1 ? i + 1 : front)].p2);
    }
    double re = 0;
    for(int i = 1; i < j - 1; ++i)
        re += (aa[i] - aa[0]) * (aa[i + 1] - aa[0]);
    return abs(re * 0.5);
}

double calipers(point aa[])
//旋转卡壳
{
    double re = 0;
    aa[n] = aa[0];
    int now = 1;
    for(int i = 0; i < n; ++i)
    {
        while((aa[i + 1] - aa[i]) * (aa[now + 1] - aa[i]) >
        (aa[i + 1] - aa[i]) * (aa[now] - aa[i]))
            now = (now == n - 1 ? 0 : now + 1);
        re = max(re, getdist2(aa[now] - aa[i]));
    }
    return re;
}

line bisector(const point &aa, const point &bb)
//中垂线(正方形)
{
    double mx, my;
    mx = (aa.x + bb.x) / 2;
    my = (aa.y + bb.y) / 2;
    line cc;
    cc.p1.x = mx - (aa.y - bb.y) / 2;
    cc.p1.y = my + (aa.x - bb.x) / 2;
    cc.p2 = aa + bb - cc.p1;
    cc.getang();
    return cc;
}

point rotate(const point &p, double cost, double sint)
//逆时针向量旋转
{
    double x = p.x, y = p.y;
    return point(x * cost - y * sint, x * sint + y * cost);
}

void getpoint(circle c1, circle c2)
//已确保两圆有交点时求出两圆交点
{
    long double dab = sqrt(getdist2(point(c1.x, c1.y) - 
    point(c2.x, c2.y)));
    if(c1.r > c2.r)
        swap(c1, c2);
    long double cost = (c1.r * c1.r + dab * dab - c2.r * c2.r) / 
    (c1.r * dab * 2);
    long double sint = sqrt(1 - cost * cost);
    point re = rotate(point(c2.x, c2.y) - point(c1.x, c1.y), cost, sint);
    re.x = c1.x + re.x * (c1.r / dab);
    re.y = c1.y + re.y * (c1.r / dab);
    point re2 = rotate(point(c2.x, c2.y) - point(c1.x, c1.y), cost, -sint);
    re2.x = c1.x + re2.x * (c1.r / dab);
    re2.y = c1.y + re2.y * (c1.r / dab);
}

double mycos(double B, double C, double A)
//余弦定理 给定边长
{
    return (B * B + C * C - A * A) / (B * C * 2);
}

double mycos2(const point &aa, const point &bb, const point &cc)
//余弦定理 给定点坐标
{
    double C2 = getdist2(aa - bb), A2 = getdist2(bb - cc), B2 = getdist2(cc - aa);
    return (B2 + C2 - A2) / (sqrt(B2 * C2) * 2);
}

point mirror_point(const point &aa, const point &bb, const point &cc)
//轴对称点 也可用来求垂足
{
    double cost, sint;
    cost = mycos2(bb, cc, aa);
    sint = sqrt(1 - cost * cost);
    if((cc - bb) * (aa - bb) > 0)
        sint = -sint;
    point re;
    re = rotate(aa - bb, cost, sint) + bb;
    re = rotate(re - bb, cost, sint) + bb;
    return re;
}

double area(const point &aa, const point &bb, const point &cc)
//求三角形面积
{
    return abs((bb - aa) * (cc - aa) / 2); 
}

void max_triangle(const point aa[])
//求凸包上最大三角形
{
    double ans = 0;
    int p1, p2, p3;
    for(int i = 0; i < n; ++i)
    {
        int j = (i + 1) % n;
        int k = (j + 1) % n;
        while(k != i && area(aa[i], aa[j], aa[k]) < area(aa[i], aa[j], aa[(k + 1) % n]))
            k = (k + 1) % n;
        if(k == i)
            continue;
        int kk = (k + 1) % n;
        while(j != kk && k != i)
        {
            if(ans < area(aa[i], aa[j], aa[k]))
            {
                ans = area(aa[i], aa[j], aa[k]);//三角形面积
                p1 = i;
                p2 = j;
                p3 = k;
            }
            while(k != i && area(aa[i], aa[j], aa[k]) < area(aa[i], aa[j], aa[(k + 1) % n]))
                k = (k + 1) % n;
            j = (j + 1) % n;
        }
    }
    point q1, q2, q3;
    q1 = b[p2] + b[p3] - b[p1];
    q2 = b[p1] + b[p3] - b[p2];
    q3 = b[p1] + b[p2] - b[p3];
    //三角形上三个点
}

void get_panel(const point3 &p1, const point3 &p2, const point3 &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 = (-(A * p1.x + B * p1.y + C * p1.z));
}

double dist_panel(const point3 &pt, double A, double B, double C, double D)
//点到平面距离
{
    return abs(A * pt.x + B * pt.y + C * pt.z + D) / 
    sqrt(A * A + B * B + C * C);
}

void hints()
{
    //皮克定理 ans = s - point / 2 + 1
    //atan2一定要先作差 再取模 再取劣弧
    //精度不够时尝试用余弦定理替代三角函数
    //算出sin cos tan中的一种后 其他的直接公式算出
    /*
    四面体重心
    x1=(s1*x1+s2*x2+s3*x3+s4*x4)/(s1+s2+s3+s4);
    y1=(s1*y1+s2*y2+s3*y3+s4*y4)/(s1+s2+s3+s4);
    z1=(s1*z1+s2*z2+s3*z3+s4*z4)/(s1+s2+s3+s4);
    */
}

int main()
{
    return 0;
}