计算几何入门
计算几何基础知识
目录
格式
写在最前
做计算几何题目，一般是有固定的模板和形式，以方便处理几何对象。一般规定为一下几点：
- 写全局函数，而非类方法，结构体只存数据；
- 每个函数标注依赖于哪些函数，且尽量减少依赖；
- 用较简略的名字，同时传值而非传 \(\texttt {const}\) 引用。
几何对象
struct Point { double x, y; }; // 点
using Vec = Point; // 向量
struct Line { Point P; Vec v; }; // 直线（点向式）
struct Seg { Point A, B; }; // 线段（存两个端点）
struct Circle { Point O; double r; }; // 圆（存圆心和半径）
当然，向量和点有时候是可以混用的，学过高中数学的同学都知道。而这样处理的目的是为了使很多函数变得简洁方便，如果后面有特殊情况，再补档。（省流：特殊在对称部分）
常用常量
const Point O = {0, 0}; // 原点
const Line Ox = {O, {1, 0}}, Oy = {O, {0, 1}}; // 坐标轴
const double PI = acos(-1), EPS = 1e-9; // pi 和 精度
关于 \(\pi\) 的求法，会在下文提到。
浮点比较
bool eq(double a, double b) { return abs(a - b) < EPS; } // ==
bool gt(double a, double b) { return a - b > EPS; } // >
bool lt(double a, double b) { return a - b < -EPS; } // <
bool ge(double a, double b) { return a - b > -EPS; } // >=
bool le(double a, double b) { return a - b < EPS; } // <=
由于很多计算几何题目，涉及到答案比较的问题，且对于答案的精度要求高，所以我们会做上述的工作。
求 \(\pi\) 的方法
考虑 \(\cos(\pi) = -1\) 所以易得 \(\pi = \arccos(-1)\)，代码中，我们一般写成 PI = arccos(-1)。
向量相关
Vec r90a(Vec v) { return {-v.y, v.x}; } // 逆时针旋转90度的向量
Vec r90c(Vec v) { return {v.y, -v.x}; } // 顺时针旋转90度的向量
Vec operator + (Vec u, Vec v) { return {u.x + v.x, u.y + v.y}; } // 向量加向量
Vec operator - (Vec u, Vec v) { return {u.x - v.x, u.y - v.y}; } // 向量减向量
Vec operator * (double k, Vec v) { return {k * v.x, k * v.y}; } // 数乘
double operator * (Vec u, Vec v) { return u.x * v.x + u.y * v.y; } // 点乘
double operator ^ (Vec u, Vec v) { return u.x * v.y - u.y * v.x; } // 叉乘
double len(Vec v) { return sqrt(v.x * v.x + v.y * v.y); } // 向量长度
double slope(Vec v) { return v.y / v.x; } // 斜率 // NOTE 不要用isinf判断斜率不存在，用后面的paral_y
// 两向量的夹角余弦
// DEPENDS len, V * V
double cos_t(Vec u, Vec v) { return u * v / len(u) / len(v); }
// 归一化向量（与原向量方向相同的单位向量）
// DEPENDS len
Vec norm(Vec v) { return {v.x / len(v), v.y / len(v)}; }
// 与原向量平行且横坐标大于等于0的单位向量
// DEPENDS d * V, len
Vec pnorm(Vec v) { return (v.x < 0 ? -1 : 1) / len(v) * v; }
// 线段的方向向量
// DEPENDS V - V
// NOTE 直线的方向向量直接访问属性v
Vec dvec(Seg l) { return l.B - l.A; }
直线相关
// 两点式直线
// DEPENDS V - V
Line line(Point A, Point B) { return {A, B - A}; }
// 斜截式直线
Line line(double k, double b) { return {{0, b}, {1, k}}; }
// 点斜式直线
Line line(Point P, double k) { return {P, {1, k}}; }
// 线段所在直线
// DEPENDS V - V
Line line(Seg l) { return {l.A, l.B - l.A}; }
// 给定直线的横坐标求纵坐标
// NOTE 请确保直线不与 y 轴平行
double at_x(Line l, double x) { return l.P.y + (x - l.P.x) * l.v.y / l.v.x; }
// 给定直线的纵坐标求横坐标
// NOTE 请确保直线不与 x 轴平行
double at_y(Line l, double y) { return l.P.x - (y + l.P.y) * l.v.x / l.v.y; }
// 点到直线的垂足
// DEPENDS V - V, V * V, d * V
Point pedal(Point P, Line l) { return l.P - (l.P - P) * l.v / (l.v * l.v) * l.v; }
// 过某点作直线的垂线
// DEPENDS r90c
Line perp(Line l, Point P) { return {P, r90c(l.v)}; }
// 角平分线
// DEPENDS V + V, len, norm
Line bisec(Point P, Vec u, Vec v) { return {P, norm(u) + norm(v)}; }
线段相关
// 线段的方向向量
// DEPENDS V - V
// NOTE 直线的方向向量直接访问属性 v
Vec dvec(Seg l) { return l.B - l.A; }
// 线段中点
Point midp(Seg l) { return {(l.A.x + l.B.x) / 2, (l.A.y + l.B.y) / 2}; }
// 线段中垂线
// DEPENDS r90c, V - V, midp
Line perp(Seg l) { return {midp(l), r90c(l.B - l.A)}; }
几何对象之间的关系
// 向量是否互相垂直
// DEPENDS eq, V * V
bool verti(Vec u, Vec v) { return eq(u * v, 0); }
// 向量是否互相平行
// DEPENDS eq, V ^ V
bool paral(Vec u, Vec v) { return eq(u ^ v, 0); }
// 向量是否与x轴平行
// DEPENDS eq
bool paral_x(Vec v) { return eq(v.y, 0); }
// 向量是否与y轴平行
// DEPENDS eq
bool paral_y(Vec v) { return eq(v.x, 0); }
// 点是否在直线上
// DEPENDS eq
bool on(Point P, Line l) { return eq((P.x - l.P.x) * l.v.y, (P.y - l.P.y) * l.v.x); }
// 点是否在线段上
// DEPENDS eq, len, V - V
bool on(Point P, Seg l) { return eq(len(P - l.A) + len(P - l.B), len(l.A - l.B)); }
// 两个点是否重合
// DEPENDS eq
bool operator==(Point A, Point B) { return eq(A.x, B.x) && eq(A.y, B.y); }
// 两条直线是否重合
// DEPENDS eq, on(L)
bool operator==(Line a, Line b) { return on(a.P, b) && on(a.P + a.v, b); }
// 两条线段是否重合
// DEPENDS eq, P == P
bool operator==(Seg a, Seg b) { return (a.A == b.A && a.B == b.B) || (a.A == b.B && a.B == b.A); }
// 以横坐标为第一关键词、纵坐标为第二关键词比较两个点
// DEPENDS eq, lt
bool operator<(Point A, Point B) { return lt(A.x, B.x) || (eq(A.x, B.x) && lt(A.y, B.y)); }
// 直线与圆是否相切
// DEPENDS eq, V ^ V, len
bool tangency(Line l, Circle C) { return eq(abs((C.O ^ l.v) - (l.P ^ l.v)), C.r * len(l.v)); }
// 圆与圆是否相切
// DEPENDS eq, V - V, len
bool tangency(Circle C1, Circle C2) { return eq(len(C1.O - C2.O), C1.r + C2.r); }
距离
// 两点间的距离
// DEPENDS len, V - V
double dis(Point A, Point B) { return len(A - B); }
// 点到直线的距离
// DEPENDS V ^ V, len
double dis(Point P, Line l) { return abs((P ^ l.v) - (l.P ^ l.v)) / len(l.v); }
// 平行直线间的距离
// DEPENDS d * V, V ^ V, len, pnorm
// NOTE 请确保两直线是平行的
double dis(Line a, Line b) { return abs((a.P ^ pnorm(a.v)) - (b.P ^ pnorm(b.v))); }
平移和旋转
// 平移
// DEPENDS V + V
Line operator + (Line l, Vec v) { return {l.P + v, l.v}; }
Seg operator + (Seg l, Vec v) { return {l.A + v, l.B + v}; }
// 旋转
// DEPENDS V + V, V - V
Point rotate(Point P, double rad) { return {cos(rad) * P.x - sin(rad) * P.y, sin(rad) * P.x + cos(rad) * P.y}; }
Point rotate(Point P, double rad, Point C) { return C + rotate(P - C, rad); } // DEPENDS ^1
Line rotate(Line l, double rad, Point C = O) { return {rotate(l.P, rad, C), rotate(l.v, rad)}; } // DEPENDS ^1, ^2
Seg rotate(Seg l, double rad, Point C = O) { return {rotate(l.A, rad, C), rotate(l.B, rad, C)}; } // DEPENDS ^1, ^2
对称
这里是本文中唯一一处点和向量表现不同的地方。
// 对称
// 关于点对称
Point reflect(Point A, Point P) { return {P.x * 2 - A.x, P.y * 2 - A.y}; }
Line reflect(Line l, Point P) { return {reflect(l.P, P), l.v}; } // DEPENDS ^1
Seg reflect(Seg l, Point P) { return {reflect(l.A, P), reflect(l.B, P)}; } // DEPENDS ^1
// 关于直线对称
// DEPENDS V - V, V * V, d * V, pedal
// NOTE 向量和点在这里的表现不同，求向量关于某直线的对称向量需要用 reflect_v
Point reflect(Point A, Line ax) { return reflect(A, pedal(A, ax)); } // DEPENDS ^1
Vec reflect_v(Vec v, Line ax) { return reflect(v, ax) - reflect(O, ax); } // DEPENDS ^1, ^4
Line reflect(Line l, Line ax) { return {reflect(l.P, ax), reflect_v(l.v, ax)}; } // DEPENDS ^1, ^4, ^5
Seg reflect(Seg l, Line ax) { return {reflect(l.A, ax), reflect(l.B, ax)}; } // DEPENDS ^1, ^4
交点
// 直线与直线交点
// DEPENDS eq, d * V, V * V, V + V, V ^ V
vector <Point> inter(Line a, Line b) {
double c = a.v ^ b.v;
if (eq(c, 0)) return {};
Vec v = 1 / c * Vec{a.P ^ (a.P + a.v), b.P ^ (b.P + b.v)};
return {{v * Vec{-b.v.x, a.v.x}, v * Vec{-b.v.y, a.v.y}}};
}
// 直线与圆交点
// DEPENDS eq, gt, V + V, V - V, V * V, d * V, len, pedal
vector <Point> inter(Line l, Circle C) {
Point P = pedal(C.O, l);
double h = len(P - C.O);
if (gt(h, C.r)) return {};
if (eq(h, C.r)) return {P};
double d = sqrt(C.r * C.r - h * h);
Vec vec = d / len(l.v) * l.v;
return {P + vec, P - vec};
}
// 圆与圆的交点
// DEPENDS eq, gt, V + V, V - V, d * V, len, r90c
vector <Point> inter(Circle C1, Circle C2) {
Vec v1 = C2.O - C1.O, v2 = r90c(v1);
double d = len(v1);
if (gt(d, C1.r + C2.r) || gt(abs(C1.r - C2.r), d)) return {};
if (eq(d, C1.r + C2.r) || eq(d, abs(C1.r - C2.r))) return {C1.O + C1.r / d * v1};
double a = ((C1.r * C1.r - C2.r * C2.r) / d + d) / 2;
double h = sqrt(C1.r * C1.r - a * a);
Vec av = a / len(v1) * v1, hv = h / len(v2) * v2;
return {C1.O + av + hv, C1.O + av - hv};
}
多边形
三角形四心
// 三角形的重心
Point barycenter(Point A, Point B, Point C) {return {(A.x + B.x + C.x) / 3, (A.y + B.y + C.y) / 3};}
// 三角形的外心
// DEPENDS r90c, V * V, d * V, V - V, V + V
// NOTE 给定圆上三点求圆，要先判断是否三点共线
Point circumcenter(Point A, Point B, Point C) {
double a = A * A, b = B * B, c = C * C;
double d = 2 * (A.x * (B.y - C.y) + B.x * (C.y - A.y) + C.x * (A.y - B.y));
return 1 / d * r90c(a * (B - C) + b * (C - A) + c * (A - B));
}
// 三角形的内心
// DEPENDS len, d * V, V - V, V + V
Point incenter(Point A, Point B, Point C) {
double a = len(B - C), b = len(A - C), c = len(A - B);
double d = a + b + c;
return 1 / d * (a * A + b * B + c * C);
}
// 三角形的垂心
// DEPENDS V * V, d * V, V - V, V ^ V, r90c
Point orthocenter(Point A, Point B, Point C) {
double n = B * (A - C), m = A * (B - C);
double d = (B - C) ^ (A - C);
return 1 / d * r90c(n * (C - B) - m * (C - A));
}
空间几何
凸包
半平面交
最小覆盖圆
最小覆盖圆。详见随机增量法。
#include <bits/stdc++.h>
using namespace std;
const double eps = 1e-8;
struct node {
double x, y;
friend node operator + (node a, node b) {return node{a.x + b.x, a.y + b.y};}
friend node operator - (node a, node b) {return node{a.x - b.x, a.y - b.y};}
friend node operator * (node a, double t) {return node{a.x * t, a.y * t};}
friend node operator / (node a, double t) {return node{a.x / t, a.y / t};}
} p[100005], ans, now;
int n;
double ansr, nowr;
node rotate_90(node a) {return node{a.y, -a.x};}
node mid(node a, node b) {return node{(a.x + b.x) / 2, (a.y + b.y) / 2};}
double dist(node a, node b) {return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) *(a.y - b.y));}
bool check(node a, double r, node b) {
if (r - dist(a, b) >= -eps) return true;
return false;
}
double cross(node a, node b) {return a.x * b.y - a.y * b.x;}
node get_intersection(node p1, node v1, node p2, node v2) {
double t = cross((p2 - p1), v2) / (cross(v1, v2));
return p1 + v1 * t;
}
node get_circle(node a, node b, node c) {return get_intersection(mid(a, b), rotate_90(b - a), mid(a, c), rotate_90(c - a));}
mt19937 rnd(chrono::system_clock::now().time_since_epoch().count());
signed main() {
scanf("%d", &n);
for (int i = 1; i <= n; i++) scanf("%lf%lf", &p[i].x, &p[i].y);
shuffle(p + 1, p + 1 + n, rnd);
ans = p[1];
ansr = 0;
for (int i = 2; i <= n; i++) {
if (check(ans, ansr, p[i])) continue;
now = p[i];
nowr = 0;
for (int j = 1; j <= i - 1; j++) {
if (check(now, nowr, p[j])) continue;
now = mid(p[i], p[j]);
nowr = dist(p[i], p[j]) / 2;
for (int k = 1; k <= j - 1; k++) {
if (check(now, nowr, p[k])) continue;
now = get_circle(p[i], p[j], p[k]);
nowr = dist(now, p[i]);
}
}
ans = now;
ansr = nowr;
}
printf("%.10lf\n", ansr);
printf("%.10lf %.10lf\n", ans.x, ans.y);
return 0;
}
