计算几何入门
计算几何入门，一窍不通
参考资料：
一个有趣的入门做题网站：CGL_all < Library < Courses | Aizu Online Judge。
实际上还有挺多的，一时间到处学板子感觉不够理智，学一点是一点吧。
精度
众所周知，计算几何常用浮点数，根据计算机储存和运算浮点数的原理，精度问题在所难免。
为减小精度误差，常不直接使用 > = < 等比较符号，而是设置一个 \(\text{eps}\)（根据题目要求设定，太大太小都是问题）。
对应的运算：
- \(a=b\)：\(|a-b|<\text{eps}\)。
- \(a<b\)：\(a<b-\text{eps}\)；\(a\leq b\)：\(a<b+\text{eps}\)。
- \(a>b\)：\(a>b+\text{eps}\)；\(a\geq b\)：\(a>b-\text{eps}\)。
点和向量
使用一个 struct 封装，点和向量都可以用坐标表示。
const double eps = 1e-10;
typedef double Db;
struct Pt {Db x, y; Pt(Db x_ = 0, Db y_ = 0) {x = x_, y = y_;}};
简单学习了一下向量外积，也即叉乘。下文令 \(\theta =\left \langle \vec a,\vec b\right \rangle\)，\(\vec a=(x_1,y_1),\vec b=(x_2,y_2)\)。
外积得到的是一个向量，方向垂直于两个向量所在平面，具体可用右手定则判定正负。模长为 \(|\vec a||\vec b|\sin \theta\)。
集合意义是两个向量围成的平行四边形的有向面积，应用广泛，坐标运算为 \(|\vec a\times \vec b|=x_1y_2-x_2y_1\)。
因为向量的方向性不满足交换律，而且有 \(\vec a\times \vec b=-\vec b\times \vec a\)。
然后就有亿些基本运算。
inline Pt operator + (const Pt &a, const Pt &b) {return Pt(a.x + b.x, a.y + b.y);}
inline Pt operator - (const Pt &a, const Pt &b) {return Pt(a.x - b.x, a.y - b.y);}
inline Pt operator * (const Pt &a, const Db &v) {return Pt(a.x * v, a.y * v);}
inline Pt operator * (const Db &v, const Pt &a) {return Pt(a.x * v, a.y * v);}
inline Pt operator / (const Pt &a, const Db &v) {return Pt(a.x / v, a.y / v);}
inline bool operator == (const Pt &a, const Pt &b) {return fabs(a.x - b.x) < eps && fabs(a.y - b.y) < eps;}
inline bool operator != (const Pt &a, const Pt &b) {return ! (a == b);}
inline double operator * (const Pt &a, const Pt &b) {return a.x * b.x + a.y * b.y;}
inline double operator ^ (const Pt &a, const Pt &b) {return a.x * b.y - a.y * b.x;}
inline double len(const Pt &a) {return sqrt(a.x * a.x + a.y * a.y);}
据此能够得到 向量夹角 和 旋转 的方法：
inline double angle(const Pt &a, const Pt &b) {return acos(a * b / len(a) / len(b));}
inline Pt rotate(const Pt &a, const double &t) {
double s = sin(t), c = cos(t);
return Pt(a.x * c - a.y * s, a.x * s + a.y * c);
}
直线与线段
直线常用 方向向量 + 直线上任意一点 表示，线段则是熟悉的两个端点。
其中方向向量默认为单位方向向量，直接 \(\vec a\gets \dfrac{\vec a}{|\vec a|}\) 即可。
两者在多方面是类似的，所以干脆直接用三个 Pt 一起封装。（\(\vec d\) 表示单位方向向量）
struct Ln {Pt a, b, d;};
引入四种点与直线的运算：
- 点到直线的垂足：\(\vec {OQ}=\vec {OA}+(\vec {AP}\cdot \vec{d})\times \vec d\)。
- 点关于直线对称点：求出垂足后用中点坐标。
- 点与直线的位置关系：利用叉积。共线时叉积 \(=0\)，不共线时根据叉积正负看顺逆时针。线段可以直接利用模长。
- 点到直线距离：利用叉积求面积，即：\(\text{dis}=\left|\dfrac{\vec{AP}\times \vec{BP}}{|\vec{AB}|}\right|\)。到线段需要判断垂足是否在线段上，如果不在就取两端点距离 \(\min\)。
判断两线段是否相交，利用快速排斥实验与跨立实验。
- 快速排斥实验：定义一条线段所在矩形为完全包含该线段的极小矩形，如果矩形无交则不可能相交。
- 跨立实验：如果 \(l_1\) 的两端点在 \(l_2\) 的两侧，且 \(l_2\) 的端点在 \(l_1\) 的两侧，就判断有交，否则无交。
快速排斥实验可以通过坐标 \(\min,\max\) 的比较实现，跨立实验则可以用先前的叉积求点与线的关系。
inline bool inter(const Ln &u, const Ln &v) {
if(min(u.a.x, u.b.x) > max(v.a.x, v.b.x) + eps || max(u.a.x, u.b.x) < min(v.a.x, v.b.x) - eps) return false;
if(min(u.a.y, u.b.y) > max(v.a.y, v.b.y) + eps || max(u.a.y, u.b.y) < min(v.a.y, v.b.y) - eps) return false;
return ( ((u.a - v.a) ^ v.d) * ((u.b - v.a) ^ v.d) < eps
&& ((v.a - u.a) ^ u.d) * ((v.b - u.a) ^ u.d) < eps );
}
判断有交之后，求交点就可以把线段当直线处理。
懒得写 vec 了，设 \(p=a_1+kd_1\)，则 \((p-a_2)\times d_2=0\)，叉积满足分配律，解方程得到 \(k=\dfrac{(a_2-a_1)\times d_2}{d_1\times d_2}\)。
带回原式即可。
inline Pt intes(const Ln &u, const Ln &v) {
double k = ((v.a - u.a) ^ v.d) / (u.d ^ v.d);
return u.a + (k * u.d);
}
三角剖分求多边形面积
有了以上基础就能干些有用的活。对于多边形上点有序给出，那么：
理性感性都能理解。
二维凸包
核心性质是，逆时针方向看，对于凸包上相邻两个向量，叉积一定 \(>0\)（一般共线不用保留，所以不取等）。
将点按照 \(x\) 为第一关键字，\(y\) 为第二关键字升序排序。
那么最小和最大的点一定在凸包中，据此划分上凸包和下凸包，分别用单调栈维护即可。
const double eps = 1e-10;
typedef double Db;
const int N = 1e5 + 10;
int n, stk[N], top; bool used[N];
struct Pt {Db x, y; Pt(Db x_ = 0, Db y_ = 0) {x = x_, y = y_;}} a[N];
inline Pt operator - (const Pt &a, const Pt &b) {return Pt(a.x - b.x, a.y - b.y);}
inline double operator ^ (const Pt &a, const Pt &b) {return a.x * b.y - a.y * b.x;}
inline double len(const Pt &a) {return sqrt(a.x * a.x + a.y * a.y);}
int main() {
n = read();
rep(i, 1, n) scanf("%lf %lf", &a[i].x, &a[i].y);
sort(a + 1, a + n + 1, [](Pt c, Pt d) {return fabs(c.x - d.x) < eps ? c.y < d.y - eps : c.x < d.x - eps;});
stk[++ top] = 1;
rep(i, 2, n) {
while(top > 1 && ((a[stk[top]] - a[stk[top - 1]]) ^ (a[i] - a[stk[top]])) < eps)
used[stk[top --]] = false;
used[stk[++ top] = i] = true;
}
int ups = top;
per(i, n - 1, 1) if(! used[i]) {
while(top > ups && ((a[stk[top]] - a[stk[top - 1]]) ^ (a[i] - a[stk[top]])) < eps)
used[stk[top --]] = false;
used[stk[++ top] = i] = true;
}
double sum = 0;
rep(i, 2, top) sum += len(a[stk[i]] - a[stk[i - 1]]);
printf("%.2lf\n", sum);
return 0;
}
注意这里的凸包求出来是首尾相连的，\(\text{stk}(1)=\text{stk}(n)\)。
旋转卡壳
时刻维护两条平行线，其中一条与凸多边形的边重合，另一条则是保持与多边形有交的极远直线。
双指针顺时针旋转维护，发现最远点也是单调移动的，利用叉积算点到直线距离，就可以 \(O(n)\) 维护了。
先要求一遍凸包，目的是使点有序（或者题目本身就是求凸包直径）。
注意函数是单峰的但是可能有连续相等的情况，所以双指针一定是 \(\leq\) 就移动。也因此要特判凸包只有两个点的情况。
rep(i, 1, top) b[i] = a[stk[i]];
int t = top, p = 1; double ans = 0;
if(t == 3) ans = len(b[1] - b[2]);
else {
rep(i, 1, t - 1) {
while(dis(b[i], b[i + 1], b[p]) < dis(b[i], b[i + 1], b[p + 1]) + eps)
p = p % (t - 1) + 1;
ans = max(ans, max(len(b[p] - b[i]), len(b[p] - b[i + 1])));
}
}
printf("%.0lf\n", ans);
半平面交
一条直线能把平面分成两个半平面，对于多边形来说，它等价于逆时针得到的相邻向量的左侧半平面交（或者顺时针右侧）。
半平面交需要先进行极角排序，对直线的方向向量使用 atan2(y, x) 就能得到与 \(x\) 轴的夹角 arg。
默认求左侧半平面交，那么就利用一个单调队列，时刻维护 \(p_i\) 表示队列中 \(i\) 与 \(i-1\) 两个向量的交点。
每次加入向量时，只会影响队首和队尾，判断 \(p_r\) 与 \(p_{l+1}\) 即可，它们不能出现在当前向量的右侧。
注意一定先判断队尾，再判断队首，因为队尾要限制队首，而队尾首先得自我限制（
判断无交需要注意，根据 P4196 半平面交 - suxxsfe 的博客 的介绍得到一种方法。
发现无非是有两个互相反向平行且互相在对方右侧的向量，如果出现这种情况，对队尾的限制一定会把所有向量全部 pop。
也就是只需要判断 \(r\) 与 \(r-1\)，如果平行，若互相在对方右侧（叉乘 \(<0\)）且反向平行（点乘 \(<0\)）就无交，否则选取较左的一个向量。
巧妙巧妙。
const int N = 10 * 55;
int n, m, t, tot; Pt p[N], ans[N]; Ln a[N], que[N];
inline Pt intes(const Ln &u, const Ln &v) {
double k = ((v.a - u.a) ^ v.d) / (u.d ^ v.d);
return u.a + (u.d * k);
}
inline bool onRight(const Ln &l, const Pt &p) {return (l.d ^ (p - l.a)) < - eps;}
bool cmp(Ln c, Ln d) {return c.ang < d.ang - eps;}
bool halfPlane() {
sort(a + 1, a + t + 1, cmp);
int l, r; que[l = r = 1] = a[1];
rep(i, 2, t) {
while(l < r && onRight(a[i], p[r])) r --;
while(l < r && onRight(a[i], p[l + 1])) l ++;
que[++ r] = a[i];
if(fabs(que[r].d ^ que[r - 1].d) < eps) {
if(onRight(que[r], que[r - 1].a) && que[r].d * que[r - 1].d < - eps)
return false;
r --; if(! onRight(que[r], a[i].a)) que[r] = a[i];
}
if(l < r) p[r] = intes(que[r - 1], que[r]);
}
while(l < r && onRight(que[l], p[r])) r --;
if(r - l <= 1) return false;
p[l] = intes(que[l], que[r]); rep(i, l, r) ans[++ tot] = p[i];
return true;
}
int main() {
n = read();
rep(s, 1, n) { m = read();
rep(i, 1, m) p[i].x = read(), p[i].y = read();
rep(i, 1, m) {
int j = (i == m) ? 1 : i + 1; t ++;
a[t].a = p[i];
a[t].b = p[j];
a[t].d = p[j] - p[i]; a[t].d = a[t].d / len(a[t].d);
a[t].ang = atan2(a[t].d.y, a[t].d.x);
}
}
if(! halfPlane()) return puts("0.00"), 0;
double sum = 0;
rep(i, 1, tot) sum += (ans[i] ^ ans[i % tot + 1]);
printf("%.3lf\n", fabs(sum) / 2.0);
return 0;
}
