计算几何知识

注意

• 浮点数0，有可能输出-0。需要x = sgn(x) == 0 ? 0 : x;处理一下。
• 精度要尽量开高一点，一般1e-10比较合适。
• double atan2(double y，double x)反正切函数，能判断象限输出$(-\pi ,\pi ]$范围的角度，且避免除零问题。

向量叉乘的意义

方向：右手定则
大小：两向量组成的平行四边形的面积

直线和线段是否相交

跨立实验，judge_point其中一个为0，或者异号。

struct Point{
    double x, y;
    double operator* (const Point &p) const {
        return x*p.y - p.x*y;
    }
    Point operator- (const Point &p) const {
        return Point{x-p.x, y-p.y};
    }
    double len(){
        return sqrt(x*x+y*y);
    }

};
struct Line{
    Point a,b;
    Point dir() const{
        return a - b;
    }
    double judge_point(const Point &p) const{
        return this->dir() * (a - p);
    }
};

两线段是否相

bool iscross(Line s1,Line s2){
    return min(s1.a.x,s1.b.x) <= max(s2.a.x,s2.b.x)
        && min(s1.a.y,s1.b.y) <= max(s2.a.y,s2.b.y)
        && min(s2.a.x,s2.b.x) <= max(s1.a.x,s1.b.x)
        && min(s2.a.y,s2.b.y) <= max(s1.a.y,s1.b.y)
        && s1.judge_point(s2.a) * s1.judge_point(s2.b) < 1e-8
        && s2.judge_point(s1.a) * s2.judge_point(s1.b) < 1e-8;
}

求两直线交点

先判断两直线是否平行（方向向量叉积为0）或共线（并有一个点）。再用以下公式求

x0 = ((x3-x4) * (x2*y1 - x1*y2) - (x1-x2) * (x4*y3 - x3*y4)) / ((x3-x4) * (y1-y2) - (x1-x2) * (y3-y4));
y0 = ((y3-y4) * (y2*x1 - y1*x2) - (y1-y2) * (y4*x3 - y3*x4)) / ((y3-y4) * (x1-x2) - (y1-y2) * (x3-x4));

坐标系旋转

在平面坐标上，任意点P(x1,y1)，绕一个坐标点Q(x2,y2)旋转$\theta$角度后,新的坐标设为(x, y)的计算公式：

x= (x1 - x2)*cos(θ) - (y1 - y2)*sin(θ) + x2 ;
y= (x1 - x2)*sin(θ) + (y1 - y2)*cos(θ) + y2 ;

(原理:先把坐标系中心平移到Q点, 把P点在极坐标下旋转$\theta$, 再把坐标轴平移回去)

极角排序

通过叉乘和距离判断

bool cmp(const Point &a, const Point &b){
    Point da = a - mark;
    Point db = b - mark;
    switch(sgn(da ^ db)){
    case 1:
        return 1;
    case -1:
        return 0;
    default:
        return da.len() < db.len();
    }
}

判断一个点是否在多边形内

回转数法

double get_angle(Point a, Point b, Point c){
    a = a - c;
    b = b - c;
    return acos((a^b) / (a.len() * b.len()));  // ^是内积
}
bool is_in_ploygn(Point p){  // 回转数法 提前把v[0]放到最后面
    double tot = 0;
    for(int i=0;i<n;++i){
        if((v[i]-p)*(v[i+1]-p) > 0) tot += get_angle(v[i],v[i+1],p);
        else tot -= get_angle(v[i],v[i+1],p);
    }
    tot = fabs(tot);
    if(sgn(tot) == 0)return false; // out
    if(sgn(tot - 2 * PI) == 0) return true; // in
    if(sgn(tot - PI) == 0) return true; // 在多边形边上 不在拐角点上
    else return true; // 在多边形拐点上
}

判断一个多边形是否为凸（已考虑顺逆时针）

bool is_convex(){   // 提前先把前两个点加到v后面
    int f=0;
    for(int tmp, i=0;i<n;++i){
        tmp = sgn( (v[i] - v[i+1]) * (v[i+1] - v[i+2]) );
        if(!f) f = tmp;
        else if(tmp * f < 0) return false;
    }
    return true;
}

Pick定理

以整点坐标组成的多边形面积为

内部格点数目$i$、边上格点数目$b$。

多边形面积

任意一个多边形的面积等于按顺序求相邻两个点与原点组成的向量的叉积之和除2的绝对值(注意正负)。

二维凸包

Graham 扫描法

1. 找出最左下角的点，作为凸包第一个点
2. 以这个点为key，对剩下的点进行极角排序（画图可知，凸包中含的点序列必然在此时升序，且第一个点和最后一个点必然为凸包上点）
3. 依次枚举，判断当前点和凸包前两个点是否为右拐，右拐则前一个点不是凸包上的点，退栈。直到可以左拐或栈只有一个点为止。加入该点。
   

int id = 0;
for(int i=1;i<n;++i){ // 选出最下，最左的点
    if(v[i].y < v[id].y || v[i].y == v[id].y && v[i].x < v[id].x)id = i;
}
swap(v[0], v[id]);
key = v[0];
sort(v.begin()+1, v.end(), cmp); // 极角排序
int top = 0;
sta[0] = v[0];
for(int i=1;i<n;++i){
    while(top > 0 && sgn( (sta[top] - sta[top-1]) * (v[i] - sta[top]) ) < 0)
        top--; // 踢掉右转的
    sta[++top] = v[i];
}

上凸包，下凸包

有时候，需要单求上凸包，或下凸包。求凸包方法，与上面类似，但是一开始不是用极角排序，而是按第一关键字x，第二关键字y排序。

bool cmp(const Point &a, const Point &b){
    return a.x < b.x || a.x == b.x && a.y < b.y;
}

下凸包，即保证叉乘大于零

sta[top = 1] = p[1];
for(int i = 2; i <= ln; ++i){ // 下凸包
    while(top > 1 && ((sta[top] - sta[top-1]) ^ (p[i] - sta[top])) <= 0)--top;
    sta[++top] = p[i];
}

上凸包，即保证叉乘小于零

sta[top = 1] = p[1];
for(int i = 2; i <= ln; ++i){ // 上凸包
    while(top > 1 && ((sta[top] - sta[top-1]) ^ (p[i] - sta[top])) >= 0)--top;
    sta[++top] = p[i];
}

旋转卡壳

被两条线夹逼着整个凸包，就产生了一对对踵点。可以证明对踵点的个数不超过$3N/2$个 也就是说对踵点的个数是$O(N)$的。旋转卡壳就是用来求对踵点对的$O(n)$算法。由于，固定一条边后，求离这条边最远的点，是一个单峰函数。并且容易看出，当边逆时针转动时，最远的点必然也在逆时针转动。就可得到以下代码，通过叉积得到三角形的面积，间接比较边到点之间的距离。

int area(int i, int j, int k){
    int res = (sta[j] - sta[i]) ^ (sta[k] - sta[i]);
    return res;
}

	key = a[0]; //前面已保证a[0]是y最小，x最小的点
    sort(a + 1, a + n, cmp); // 极角排序
    int top = 0;
    sta[0] = a[0]; // 求凸包
    for(int i = 1; i < n; ++i){
        while(top > 0 && ((sta[top] - sta[top - 1]) ^ (a[i] - sta[top])) <= 0) top--;
        sta[++top] = a[i];
    }
    int ans = 0;
    sta[++top] = sta[0];
    int j = 1; // 旋转卡壳
    for(int i = 0; i < top; ++i){
        while(area(i, i + 1, j) < area(i, i + 1, j + 1) ){
            j = (j + 1) % top;
        }
        ans = max(ans, (sta[i] - sta[j]).len());
    }

扫描线

求矩形组成的面积或周长。用线段树统计目前每一段的计数，计数大于0则统计该段长度。

#include <bits/stdc++.h>
#define endl '\n'
#define IOS ios::sync_with_stdio(false); cin.tie(nullptr), cout.tie(nullptr);
using namespace std;
const int N = 2e2 + 5;
double valx[N];
struct Point{
    double x1, x2, y;
    int flag;
    bool operator< (const Point &ot) const &{
        return y < ot.y;
    }
}p[N];
#define ls rt << 1
#define rs rt << 1 | 1
double sum[N << 3];
int lazy[N << 3];
void pushup(int rt, int l, int r){
    if(lazy[rt]){
        sum[rt] = valx[r - 1] - valx[l - 1];
    }else if(r + 1 != l)sum[rt] = sum[ls] + sum[rs];
    else sum[rt] = 0;
}
void update(int rt, int l, int r, int lx, int rx, int f){
    if(lx <= l && r <= rx){
        lazy[rt] += f;
        pushup(rt, l, r);
        return;
    }
    int mid = l + r >> 1;
    if(lx < mid)update(ls, l, mid, lx, rx, f);
    if(rx > mid)update(rs, mid, r, lx, rx, f);
    pushup(rt, l, r);
}
int getpos(double x, int n){
    return lower_bound(valx, valx + n, x) - valx + 1;
}
int main(){
    IOS;
    int n, cas = 1;
    while(cin >> n && n){
        for(int i = 0; i < n; ++i){
            double x1, x2, y1, y2;
            cin >> x1 >> y1 >> x2 >> y2;
            valx[i] = x1, valx[i + n] = x2;
            p[i].x1 = x1, p[i].x2 = x2, p[i].y = y1, p[i].flag = 1;
            p[n + i].x1 = x1, p[n + i].x2 = x2, p[n + i].y = y2, p[n + i].flag = -1;
        }
        memset(sum, 0, sizeof(sum));
        memset(lazy, 0, sizeof(lazy));
        sort(valx, valx + 2 * n);
        int m = unique(valx, valx + 2 * n) - valx;
        sort(p, p + 2 * n);
        double ans = 0;
        for(int i = 0; i < 2 * n; ++i){
            //cout << p[i].y << " " << p[i].x1 << " " << p[i].x2 << " " << p[i].flag << " sum=" << sum[1] << endl;
            if(i)ans += (p[i].y - p[i - 1].y) * sum[1];
            update(1, 1, m, getpos(p[i].x1, m), getpos(p[i].x2, m), p[i].flag);
        }
        cout << "Test case #" << cas++ << "\nTotal explored area: " << fixed << setprecision(2) << ans << endl << endl;
    }
}