计算几何篇

admin1

const double eps = 1e-6;
struct Point{
    double x, y;
    Point(double x = 0, double y = 0):x(x),y(y){}
};
typedef Point Vector;
Vector operator + (Vector A, Vector B){
    return Vector(A.x+B.x, A.y+B.y);
}
Vector operator - (Point A, Point B){
    return Vector(A.x-B.x, A.y-B.y);
}
Vector operator * (Vector A, double p){
    return Vector(A.x*p, A.y*p);
}
int sgn(double x){
    if(fabs(x) < eps)
        return 0;
    if(x < 0)
        return -1;
    return 1;
}
double Dot(Vector A, Vector B){
    return A.x*B.x + A.y*B.y;
}
double Cross(Vector A, Vector B){
    return A.x*B.y-A.y*B.x;
}
double Length(Vector A){
    return sqrt(Dot(A, A));
}
Vector Normal(Vector A){//向量A左转90°的单位法向量
    double L = Length(A);
    return Vector(-A.y/L, A.x/L);
}
struct Line{
    Point p;//直线上任意一点
    Vector v;//方向向量，它的左边就是对应的半平面
    double ang;//极角，即从x轴正半轴旋转到向量v所需要的角（弧度）
    Line(){}
    Line(Point p, Vector v) : p(p), v(v){
        ang = atan2(v.y, v.x);
    }
    bool operator < (const Line& L) const {//排序用的比较运算符
        return ang < L.ang;
    }
};
//点p在有向直线L的左侧
bool OnLeft(Line L, Point p){
    return Cross(L.v, p - L.p) > 0;
}
//两直线交点。假定交点唯一存在
Point GetIntersection(Line a, Line b){
    Vector u = a.p - b.p;
    double t = Cross(b.v, u)/Cross(a.v, b.v);
    return a.p + a.v*t;
}
//半平面交的主过程
int HalfplaneIntersection(Line* L, int n, Point* poly){
    sort(L, L + n);//按照极角排序
    int fst = 0, lst = 0;//双端队列的第一个元素和最后一个元素
    Point *P = new Point[n];//p[i] 为 q[i]与q[i + 1]的交点
    Line *q = new Line[n];//双端队列
    q[fst = lst = 0] = L[0];//初始化为只有一个半平面L[0]
    for(int i = 1; i < n; ++i){
        while(fst < lst && !OnLeft(L[i], P[lst - 1])) --lst;
        while(fst < lst && !OnLeft(L[i], P[fst])) ++fst;
        q[++lst] = L[i];
        if(sgn(Cross(q[lst].v, q[lst - 1].v)) == 0){
            //两向量平行且同向，取内侧一个
            --lst;
            if(OnLeft(q[lst], L[i].p)) q[lst] = L[i];
        }
        if(fst < lst)
            P[lst - 1] = GetIntersection(q[lst - 1], q[lst]);
    }
    while(fst < lst && !OnLeft(q[fst], P[lst - 1])) --lst;
    //删除无用平面
    if(lst - fst <= 1) return 0;//空集
    P[lst] = GetIntersection(q[lst], q[fst]);//计算首尾两个半平面的交点
    //从deque复制到输出中
    int m = 0;
    for(int i = fst; i <= lst; ++i) poly[m++] = P[i];
    return m;
}