2025.8.12 计算几何
点与向量
点与向量二者类似，在平面直角坐标系中都可以用二元组 \((x,y)\) 来精确刻画. 通过一些代数运算，我们可以刻画向量之间的位置、方向关系. 这是计算几何的核心.
点击查看代码
struct pt{
double x, y;
pt(double _x = 0, double _y = 0) {x = _x, y = _y;}
inline void read() {cin >> x >> y; return;}
};
基本运算
- 向量的模长就是到原点的距离.
点击查看代码
inline double Len(const vec &a) {return sqrt(a.x * a.x + a.y * a.y);}
- 向量的辐角与复数类似，是由 \(x\) 正半轴逆时针旋转所成的角度.
点击查看代码
inline double angle(const vec &a) {return atan2(a.y, a.x);}
- 向量的加法与减法是对应坐标的加减，向量的数乘与数除是对应坐标乘或除以一个数.
点击查看代码
inline vec operator + (const vec &a, const vec &b) {return vec(a.x + b.x, a.y + b.y);}
inline vec operator - (const vec &a, const vec &b) {return vec(a.x - b.x, a.y - b.y);}
inline vec operator * (const vec &a, double b) {return vec(a.x * b, a.y * b);}
inline vec operator / (const vec &a, double b) {return vec(a.x / b, a.y / b);}
点积与叉积
对于向量 \(\vec a=(x_1,y_1),\vec b=(x_2,y_2)\)，定义点积：
\[\vec a\cdot \vec b=|\vec a||\vec b|
\cos\langle \vec a,\vec b\rangle=x_1x_2+y_1y_2
\]
- 点积具有交换律，其结果与 \(\vec a\cdot \vec b\) 的顺序无关. 而点积与 \(0\) 的大小关系可以用于判断两向量的夹角大小：\(\vec a\cdot \vec b\) 为正，夹角为锐角；\(\vec a\cdot \vec b\) 为 \(0\)，相互垂直；\(\vec a\cdot \vec b\) 为负，夹角为钝角.
点击查看代码
inline double operator * (const vec &a, const vec &b) {return a.x * b.x + a.y * b.y;}
定义叉积：
\[\vec a\times \vec b=|\vec a||\vec b|
\sin\langle \vec a,\vec b\rangle=x_1y_2-x_2y_1
\]
-
叉积具有反交换律：\(\vec a\times \vec b=-\vec b\times \vec a\)，也就是说叉积运算的两者有先后的区分. 具体来讲对于 \(\vec a\times \vec b\)，当 \(\vec b\) 在 \(\vec a\) 的逆时针方向时叉积为正，顺时针方向时叉积为负. 这个性质在求凸包、旋转卡壳、半平面交都有应用.
-
叉积绝对值的一半就是两向量围成三角形的面积.
点击查看代码
inline double operator ^ (const vec &a, const vec &b) {return a.x * b.y - a.y * b.x;}
inline double S(const pt &x, const pt &y, const pt &z) {return Abs((y - x) ^ (z - x));}
旋转
根据 \(\sin\) 与 \(\cos\) 的和角公式，可以推得向量 \(\vec a=(x,y)\) 逆时针旋转 \(\theta\) 的坐标变换公式：
\[\vec {a'}=(x\cos\theta-y\sin\theta,y\cos\theta+x\sin\theta)
\]
向量绕点旋转可以先将参考系换为原点，逆时针旋转之后再换回来.
点击查看代码
inline vec rotate(vec a, double theta) {
double x = a.x * cos(theta) - a.y * sin(theta);
double y = a.y * cos(theta) + a.x * sin(theta);
return vec(x, y);
}
inline pt rotate_P(pt a, pt b, double theta) {return rotate(a - b, theta) + b;}
直线与线段
直线与线段都由两个有序的点对 \((s,t)\) 确定，方向 \(s\rightarrow t\). 这是为了方便确定半平面的位置. 两者的表示没有任何区别，各种操作需要手动分讨.
点击查看代码
struct line {
pt s, t;
line(pt _s = pt(0, 0), pt _t = pt(0, 0)) {s = _s, t = _t;}
};
位置关系
- 点是否在直线/线段上：前者可以用叉积是否为 \(0\) 实现；后者可以用点积实现.
点击查看代码
inline bool judge_PL(pt a, line b) {return !dcmp((a - b.s) ^ (b.t - b.s));}
inline bool judge_PS(pt a, line b) {return !dcmp((a - b.s) ^ (b.t - b.s)) && dcmp((a - b.s) * (a - b.t)) <= 0;}
- 点到直线/线段的距离：前者可以用三角形面积除以底；后者先用点积判断夹角，分讨求解.
点击查看代码
inline double dis_PL(pt a, line b) {return Abs((a - b.s) ^ (a - b.t)) / Len(b.t - b.s);}
inline double dis_PS(pt a, line b) {
pt x = a - b.s, y = a - b.t, z = b.t - b.s;
if(dcmp(x * z) < 0) return Len(x);
if(dcmp(y * z) > 0) return Len(y);
return dis_PL(a, b);
}
- 点到直线的垂足与点到线段距离最短的点：前者可以用点积直接求出投影向量；后者与点到线段距离的分讨类似.
点击查看代码
inline pt Footprint(pt a, line b) {
pt x = a - b.s, y = a - b.t, z = b.t - b.s;
double s1 = x * z, s2 = -1.0 * (y * z);
return b.s + z * (s1 / (s1 + s2));
}
inline pt point_PS(pt a, line b) {
pt x = a - b.s, y = a - b.t, z = b.t - b.s;
if(dcmp(x * z) < 0) return b.s;
if(dcmp(y * z) > 0) return b.t;
return Footprint(a, b);
}
- 判断线段与直线/线段是否相交：前者可以先求出线线交点再判断位置关系；后者先判断值域是否有交，再用叉积判断两端点是否在同侧.
点击查看代码
inline pt cross_LL(line a, line b) {
pt x = a.t - a.s, y = b.t - b.s, z = a.s - b.s;
return a.s + x * ((y ^ z) / (x ^ y));
}
inline bool judge_cross_SL(line a, line b) {
if(!dcmp((a.t - a.s) ^ (b.t - a.s))) return false;
return judge_PS((cross_LL(a, b)), a);
}
inline double maxx(const line L) {return max(L.s.x, L.t.x);}
inline double maxy(const line L) {return max(L.s.y, L.t.y);}
inline double minx(const line L) {return min(L.s.x, L.t.x);}
inline double miny(const line L) {return min(L.s.y, L.t.y);}
inline bool judge_cross_SS(line a, line b) {
if(maxx(a) < minx(b) || maxx(b) < minx(a)) return false;
if(maxy(a) < miny(b) || maxy(b) < miny(a)) return false;
double s1 = (a.t - a.s) ^ (b.s - a.s), s2 = (a.t - a.s) ^ (b.t - a.s);
double s3 = (b.t - b.s) ^ (a.s - b.s), s4 = (b.t - b.s) ^ (a.t - b.s);
return dcmp(s1) * dcmp(s2) <= 0 && dcmp(s3) * dcmp(s4) <= 0;
}
多边形
多边形是计算几何题中最常研究的对象，其中的点通常逆时针给出.
struct polygon{
vector<pt> pts;
inline pt& operator [](int x) {return pts[x];}
inline void read(int n) {
pt cur;
for(int i = 0; i < n; i++) cur.read(), pts.push_back(cur);
}
inline void Clear() {pts.clear(); return;}
inline void insert(pt p) {pts.push_back(p); return;}
};
周长与面积
周长是简单的，因为
板子
以下板子封装的内容包括：
- 计算几何常用基本函数
- 多边形类
- Andrew 求凸包、旋转卡壳、半平面交
- 自适应辛普森积分法
#include<bits/stdc++.h>
using namespace std;
const double eps = 1e-8;
inline int dcmp(double x) {return (x < -eps) ? -1 : (x > eps ? 1 : 0);}
inline double Abs(double x) {return x * dcmp(x);}
namespace CG{
struct pt{
double x, y;
pt(double _x = 0, double _y = 0) {x = _x, y = _y;}
inline void read() {cin >> x >> y; return;}
};
inline bool cmpx(const pt &a, const pt &b) {return (a.x != b.x) ? a.x < b.x : a.y < b.y;}
typedef pt vec;
inline double Len(const vec &a) {return sqrt(a.x * a.x + a.y * a.y);}
inline double angle(const vec &a) {return atan2(a.y, a.x);}
inline vec operator + (const vec &a, const vec &b) {return vec(a.x + b.x, a.y + b.y);}
inline vec operator - (const vec &a, const vec &b) {return vec(a.x - b.x, a.y - b.y);}
inline vec operator * (const vec &a, double b) {return vec(a.x * b, a.y * b);}
inline vec operator / (const vec &a, double b) {return vec(a.x / b, a.y / b);}
inline double operator * (const vec &a, const vec &b) {return a.x * b.x + a.y * b.y;}
inline double operator ^ (const vec &a, const vec &b) {return a.x * b.y - a.y * b.x;}
inline vec rotate(vec a, double theta) {
double x = a.x * cos(theta) - a.y * sin(theta);
double y = a.y * cos(theta) + a.x * sin(theta);
return vec(x, y);
}
inline pt rotate_P(pt a, pt b, double theta) {return rotate(a - b, theta) + b;}
struct line {
pt s, t;
line(pt _s = pt(0, 0), pt _t = pt(0, 0)) {s = _s, t = _t;}
};
inline double maxx(const line L) {return max(L.s.x, L.t.x);}
inline double maxy(const line L) {return max(L.s.y, L.t.y);}
inline double minx(const line L) {return min(L.s.x, L.t.x);}
inline double miny(const line L) {return min(L.s.y, L.t.y);}
inline double ang(const line &L) {return angle(L.t - L.s);}
inline bool operator < (const line &a, const line &b) { //半平面交的极角排序
double a1 = angle(a.t - a.s), a2 = angle(b.t - b.s);
if(dcmp(a2 - a1) != 0) return dcmp(a2 - a1) > 0;
return dcmp((b.t - a.s) ^ (a.t - a.s)) > 0;
}
inline bool operator == (pt a, pt b) {return !dcmp(a.x - b.x) && !dcmp(a.y - b.y);}
inline double dis_PP(pt a, pt b) {return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));}
inline bool judge_PL(pt a, line b) {return !dcmp((a - b.s) ^ (b.t - b.s));}
inline bool judge_PS(pt a, line b) {return !dcmp((a - b.s) ^ (b.t - b.s)) && dcmp((a - b.s) * (a - b.t)) <= 0;}
inline pt Footprint(pt a, line b) {
pt x = a - b.s, y = a - b.t, z = b.t - b.s;
double s1 = x * z, s2 = -1.0 * (y * z);
return b.s + z * (s1 / (s1 + s2));
}
inline pt mirror(pt a, line b) {return 2.0 * Footprint(a, b) - a;}
inline double dis_PL(pt a, line b) {return Abs((a - b.s) ^ (a - b.t)) / Len(b.t - b.s);}
inline double dis_PS(pt a, line b) {
pt x = a - b.s, y = a - b.t, z = b.t - b.s;
if(dcmp(x * z) < 0) return Len(x);
if(dcmp(y * z) > 0) return Len(y);
return dis_PL(a, b);
}
inline pt point_PS(pt a, line b) {
pt x = a - b.s, y = a - b.t, z = b.t - b.s;
if(dcmp(x * z) < 0) return b.s;
if(dcmp(y * z) > 0) return b.t;
return Footprint(a, b);
}
inline pt cross_LL(line a, line b) {
pt x = a.t - a.s, y = b.t - b.s, z = a.s - b.s;
return a.s + x * ((y ^ z) / (x ^ y));
}
inline bool judge_cross_SL(line a, line b) {
if(!dcmp((a.t - a.s) ^ (b.t - a.s))) return false;
return judge_PS((cross_LL(a, b)), a);
}
inline bool judge_cross_SS(line a, line b) {
if(maxx(a) < minx(b) || maxx(b) < minx(a)) return false;
if(maxy(a) < miny(b) || maxy(b) < miny(a)) return false;
double s1 = (a.t - a.s) ^ (b.s - a.s), s2 = (a.t - a.s) ^ (b.t - a.s);
double s3 = (b.t - b.s) ^ (a.s - b.s), s4 = (b.t - b.s) ^ (a.t - b.s);
return dcmp(s1) * dcmp(s2) <= 0 && dcmp(s3) * dcmp(s4) <= 0;
}
} using namespace CG;
namespace Polygon{
struct polygon{
vector<pt> pts;
inline pt& operator [](int x) {return pts[x];}
inline void read(int n) {
pt cur;
for(int i = 0; i < n; i++) cur.read(), pts.push_back(cur);
}
inline void Clear() {pts.clear(); return;}
inline int nxt(int x) {return x < pts.size() - 1 ? x + 1 : 0;}
inline int include(pt p) {
int cnt = 0;
for(int i = 0; i < pts.size(); i++) {
pt s = pts[i], t = pts[nxt(i)];
line L = line(s, t);
if(judge_PS(p, L)) return 2;
if(dcmp(p.y - miny(L)) >= 0 && dcmp(maxy(L) - p.y) > 0) {//竖直向上射线
double nx = s.x + ((p.y - s.y) / (t.y - s.y) * (t.x - s.x));
if(dcmp(nx - p.x) > 0) cnt++;
}
}
return cnt & 1;
}
inline double area() {
double ans = 0;
for(int i = 1; i < pts.size() - 1; i++) ans += (pts[i] - pts[0]) ^ (pts[nxt(i)] - pts[0]);
return Abs(ans) / 2;
}
inline double peri() {
double ans = 0;
for(int i = 0; i < pts.size(); i++) ans += dis_PP(pts[i], pts[nxt(i)]);
return ans;
}
inline bool Left(pt x, pt l, pt r) {return dcmp((l - x) ^ (r - x)) <= 0;}
inline int include_bs(pt p) {
int n = pts.size();
if(!Left(pts[0], p, pts[1])) return 0;
if(!Left(pts[0], pts[n - 1], p)) return 0;
if(judge_PS(p, line(pts[0], pts[1]))) return 2;
if(judge_PS(p, line(pts[0], pts[n - 1]))) return 2;
int l = 1, r = n - 2, ans = 1;
while(l <= r) {
int mid = (l + r) >> 1;
if(!Left(pts[0], pts[mid], p)) l = mid + 1, ans = mid;
else r = mid - 1;
}
if(!Left(pts[ans], p, pts[nxt(ans)])) return 0;
if(judge_PS(p, line(pts[ans], pts[nxt(ans)]))) return 2;
return 1;
}
inline void insert(pt p) {pts.push_back(p); return;}
};
} using namespace Polygon;
namespace Hull{
inline void Andrew(polygon &p) {
vector<pt> q(p.pts.size() << 1);
sort(p.pts.begin(), p.pts.end(), cmpx);
int tp = 0, n = p.pts.size();
q[0] = p[0];
for(int i = 1; i < n; i++) {
while(tp && dcmp((q[tp - 1] - q[tp]) ^ (p[i] - q[tp])) >= 0) tp--;
q[++tp] = p[i];
}
int nw = tp;
for(int i = n - 2; i >= 0; i--) {
while(tp > nw && dcmp((q[tp - 1] - q[tp]) ^ (p[i] - q[tp])) >= 0) tp--;
q[++tp] = p[i];
}
if(n > 1) --tp;
p.Clear(); for(int i = 0; i <= tp; i++) p.insert(q[i]);
return;
}
inline double S(const pt &x, const pt &y, const pt &z) {return Abs((y - x) ^ (z - x));}
inline double diam(polygon &p) {
int n = p.pts.size(); double ans = 0;
for(int i = 0, j = 1; i < n; i++) {
for(; ; j = p.nxt(j)) if(dcmp(S(p[j], p[i], p[p.nxt(i)]) - S(p[p.nxt(j)], p[i], p[p.nxt(i)])) >= 0) break;
ans = max(ans, max(dis_PP(p[i], p[j]), dis_PP(p[p.nxt(i)], p[j])));
}
return ans;
}
inline polygon SI(vector<line> q) {
int n = q.size();
sort(q.begin(), q.end());
vector<line> li(n + 1), L(n + 1); vector<pt> p(n + 1);
int len = 0;
for(int i = 0; i < n; i++) {
if(i > 0 && !dcmp(ang(q[i]) - ang(q[i - 1]))) continue;
L[++len] = q[i];
}
int l = 1, r = 0;
for(int i = 1; i <= len; i++) {
while(l < r && dcmp((L[i].t - p[r]) ^ (L[i].s - p[r])) > 0) --r;
while(l < r && dcmp((L[i].t - p[l + 1]) ^ (L[i].s - p[l + 1])) > 0) ++l;
li[++r] = L[i]; if(r > l) p[r] = cross_LL(li[r], li[r - 1]);
}
while(l < r && dcmp((L[l].t - p[r]) ^ (li[l].s - p[r])) > 0) --r;
while(l < r && dcmp((li[r].t - p[l + 1]) ^ (li[r].s - p[l + 1])) > 0) l++;
p[r + 1] = cross_LL(li[r], li[l]), ++r;
polygon P;
for(int i = l + 1; i <= r; i++) P.insert(p[i]);
return P;
}
}
using Hull :: Andrew;
using Hull :: diam;
using Hull :: SI;
namespace Simpson{
inline double f(double x) {return ;} //欲求函数
inline double simpson(double l, double r) {return (f(l) + 4 * f((l + r) / 2) + f(r)) * (r - l) / 6;}
inline double asr(double l, double r, double Eps, double ans) {
double mid = (l + r) / 2;
double ansl = simpson(l, mid), ansr = simpson(mid, r);
if(Abs(ansl + ansr - ans) <= 15 * Eps) return ansl + ansr + (ansl + ansr - ans) / 15;
return asr(l, mid, Eps / 2, ansl) + asr(mid, r, Eps / 2, ansr);
}
inline double Asr(double l, double r) {return asr(l, r, eps, simpson(l, r));}
}
using Simpson :: Asr;
浙公网安备 33010602011771号