计算几何相关总结
一、基础知识
-
\(\pi = \arccos(-1)\)
-
余弦定理 \(c^2=a^2+b^2-2ab\cos\theta\)
-
浮点数的比较
const double eps=1e-8; int sign(double x) { if (fabs(x) < eps) return 0; if (x < 0) return -1; return 1; }int cmp(double x, double y) { if (fabs(x - y) < eps) return 0; if (x < y) return -1; return 1; } -
向量
-
内积、点积
double dot(Point a, Point b) { return a.x * b.x +a.y * b.y; } -
外积、叉积
double cross(Point a, Point b) { return a.x * b.y - b.x * a.y; } -
向量长度
sqrt(dot(a, a)); -
向量夹角
acos(dot(a, b) / get_length(a) / get_length(b)); -
平行四边形面积
cross(b - a, c - a); -
向量旋转
顺时针
\[\begin{bmatrix}x&y\end{bmatrix}\times\begin{bmatrix}\cos\theta&-\sin\theta\\\sin\theta&\cos\theta\end{bmatrix} \]逆时针
\[\begin{bmatrix}x&y\end{bmatrix}\times\begin{bmatrix}\cos \theta & \sin\theta\\-\sin\theta&\cos\theta\end{bmatrix} \]
-
-
直线
-
点向式 \(\boldsymbol{v}t+P_0\)。
-
判断点在直线上 \(\boldsymbol{A}\times\boldsymbol{B}=0\)。
-
相交
Point get_line_intersection(Point p, Vector v, Point q, Vector w) { Vector u = p - q; double t = cross(w, u) / cross(v, w); return p + v * t; } -
点到直线距离
double distance_to_line(Point p, Point a, Point b) { Vector v1 = b - a, v2 = p - a; return fabs(cross(v1, v2) / get_length(v1)); } -
点到线段距离
double distance_to_segment(Point p, Point a, Point b) { if (a == b) return get_length(p - a); Vector v1 = b - a, v2 = p - a, v3 = p - b; if (sign(dot(v1, v2)) < 0) return get_length(v2); if (sign(dot(v1, v3)) > 0) return get_length(v3); return distance_to_line(p, a, b); } -
点在直线上投影
double get_line_projection(Point p, Point a, Point b) { Vector v = b - a; return a + v * (dot(v, p - a) / dot(v, v)); } -
点是否在线段上
bool on_segment(Point p, Point a, Point b) { return sign(cross(p - a, p - b)) == 0 && sign(dot(p - a, p - b)) <= 0; } -
判断两线段是否相交
bool segment_intersection(Point a1, Point a2, Point b1, Point b2) { double c1 = cross(a2 - a1, b1 - a1), c2 = cross(a2 - a1, b2 - a1); double c3 = cross(b2 - b1, a2 - b1), c4 = cross(b2 - b1, a1 - b1); return sign(c1) * sign(c2) <= 0 && sign(c3) * sign(c4) <= 0; }
-
-
三角形
-
面积：叉积、海伦公式：
\[p = (a+b+c) / 2\\ S = \sqrt{p(p-a)(p-b)(p-c)} \] -
外心、内心、垂心、重心
重心是到三角形三顶点距离平方和最小的点，三角形内到三边距离之积最大的点。
-
-
普通多边形（通常按逆时针方向存储所有顶点）
-
求多边形面积：
double polygon_area(Point p[], int n) { double s = 0; for (int i = 1; i < n - 1; ++ i) s += cross(p[i] - p[0], p[i + 1] - p[0]); return s / 2; } -
判断点是否在多边形里
射线法
从该点出发任意做一条和所有边都不平行的射线。交点个数为偶数，则在多边形外；为奇数，则在多边形内。
-
判断点是否在凸多边形内
只需判断点是否在所有边左边
-
皮克定理（顶点在格点上）
\[S=a+b/2+1 \]\(a\) 表示多边形内部的点数，\(b\) 表示多边形边上的点数
-
二、平面凸包
Andrew 算法。将凸包分成上面和下面两部分。
先将点按照横坐标从小到大排序。对于上凸壳来说，从左向右扫描每一个点；对于下凸壳，从右向左扫描每一个点。用一个栈维护当前凸壳的轮廓。
inline double andrew(void) { // 求平面凸包周长
sort(a + 1, a + n + 1);
stk[++ top] = 1; used[1] = true;
for (int i = 1; i <= n; ++ i) {
while (top > 1 && sign(cross(a[stk[top]] - a[stk[top - 1]], a[i] - a[stk[top - 1]])) >= 0) {
if (sign(cross(a[stk[top]] - a[stk[top - 1]], a[i] - a[stk[top - 1]])) > 0)
used[stk[top --]] = false;
else -- top;
}
stk[++ top] = i; used[i] = true;
}
used[1] = false;
for (int i = n - 1; i >= 1; -- i) {
if (used[i]) continue;
while (top > 1 && sign(cross(a[stk[top]] - a[stk[top - 1]], a[i] - a[stk[top - 1]])) >= 0)
-- top;
stk[++ top] = i;
}
double sum = 0;
for (int i = top - 1; i >= 1; -- i) {
sum += get_dist(a[stk[i]], a[stk[i + 1]]);
}
return sum;
}
三、半平面交
先将所有直线按照倾斜角从小到大排序。然后，依次扫描每一条直线，用双端队列维护当前凸壳的轮廓。
inline void half_plane_intersection(void) { // 求半平面交面积
sort(a + 1, a + tot + 1, cmp);
hh = 0; tt = -1;
q[++ tt] = 1;
for (int i = 2; i <= tot; ++ i) {
while (hh <= tt - 1 && on_right(a[i], get_line_intersection(a[q[tt]], a[q[tt - 1]])))
-- tt;
while (hh <= tt - 1 && on_right(a[i], get_line_intersection(a[q[hh]], a[q[hh + 1]])))
++ hh;
q[++ tt] = i;
}
while (hh <= tt - 1 && on_right(a[q[tt]], get_line_intersection(a[q[hh]], a[q[hh + 1]])))
++ hh;
while (hh <= tt - 1 && on_right(a[q[hh]], get_line_intersection(a[q[tt]], a[q[tt - 1]])))
-- tt;
q[++ tt] = q[hh];
tot = 0;
for (int i = hh; i < tt; ++ i)
ver[++ tot] = get_line_intersection(a[q[i]], a[q[i + 1]]);
double res = 0;
for (int i = 2; i < tot; ++ i)
res += area(ver[1], ver[i], ver[i + 1]);
printf("%.3lf\n", res / 2);
}
四、最小圆覆盖
最小覆盖圆一定是唯一的。
依次确定最小覆盖圆上的三个点。
inline void solve(void) {
random_shuffle(a + 1, a + n + 1);
if (n == 2) {
get_circle(a[1], a[2]).print();
return;
}
Circle c = get_circle(a[1], a[2], a[3]);
for (int i = 4; i <= n; ++ i) {
if (dcmp(get_dist(a[i], c.o), c.r) <= 0) continue;
c = get_circle(a[i]);
for (int j = 1; j <= i - 1; ++ j) {
if (dcmp(get_dist(a[j], c.o), c.r) <= 0) continue;
c = get_circle(a[i], a[j]);
for (int k = 1; k <= j - 1; ++ k) {
if (dcmp(get_dist(a[k], c.o), c.r) <= 0) continue;
c = get_circle(a[i], a[j], a[k]);
}
}
}
c.print();
}
五、三维凸包
思路比二维凸包简单。
将平面分割成三角形，每一个三角形是一个平面。每次加进来一个点的时候，把新点想象成一个光源，当前的凸包是一个障碍物，那么被光源能照到的平面都要舍弃，照不到的平面不能舍弃。
double solve(void) {
plane[++ m] = Plane(1, 2, 3);
plane[++ m] = Plane(1, 3, 2);
for (int i = 4; i <= n; ++ i) {
int cnt = 0;
for (int j = 1; j <= m; ++ j) {
bool t = above(plane[j], q[i]);
if (!t) np[++ cnt] = plane[j];
for (int k = 0; k < 3; ++ k) {
g[plane[j].v[k]][plane[j].v[(k + 1) % 3]] = t;
}
}
for (int j = 1; j <= m; ++ j) {
for (int k = 0; k < 3; ++ k) {
int a = plane[j].v[k], b = plane[j].v[(k + 1) % 3];
if (g[a][b] && !g[b][a]) {
np[++ cnt] = Plane(a, b, i);
}
}
}
m = cnt;
swap(plane, np);
}
double res = 0;
for (int i = 1; i <= m; ++ i) {
res += area(plane[i]);
}
return res;
}
六、旋转卡壳
对踵点的最远距离。
决策单调性。
inline int solve(void) {
sort(a + 1, a + n + 1);
for (int i = 1; i <= n; ++ i) {
while (top > 1 && area(a[i], a[stk[top]], a[stk[top - 1]]) >= 0) {
if (area(a[i], a[stk[top]], a[stk[top - 1]]) > 0)
used[stk[top --]] = false;
else -- top;
}
stk[++ top] = i;
used[i] = true;
}
used[1] = false;
for (int i = n - 1; i >= 1; -- i) {
if (used[i]) continue;
while (top > 1 && area(a[i], a[stk[top]], a[stk[top - 1]]) >= 0)
-- top;
stk[++ top] = i;
}
-- top;
if (top <= 1)
return get_dist(a[1], a[n]);
int res = 0;
for (int i = 1, j = 3; i <= top; ++ i) {
while (area(a[stk[i]], a[stk[i + 1]], a[stk[j == top ? 1 : j + 1]]) > area(a[stk[i]], a[stk[i + 1]], a[stk[j]]))
j = j == top ? 1 : j + 1;
res = max(res, max(get_dist(a[stk[i]], a[stk[j]]), get_dist(a[stk[i + 1]], a[stk[j]])));
}
return res;
}
七、平面最近点对
这是一个分治问题，和计算几何没多大关系，但我还是想放在这里。
double solve(int l, int r) {
if (l == r) return inf;
int mid = (l + r) >> 1;
double mid_x = a[mid].x;
double res = min(solve(l, mid), solve(mid + 1, r));
{
int i = l, j = mid + 1, k = 0;
while (i <= mid && j <= r) {
if (a[i].y < a[j].y) temp[++ k] = a[i ++];
else temp[++ k] = a[j ++];
}
while (i <= mid) temp[++ k] = a[i ++];
while (j <= r) temp[++ k] = a[j ++];
for (int i = l, j = 1; i <= r; ++ i, ++ j) a[i] = temp[j];
}
int k = 0;
for (int i = l; i <= r; ++ i)
if (mid_x - res <= a[i].x && a[i].x <= mid_x + res)
temp[++ k] = a[i];
for (int i = 1; i <= k; ++ i)
for (int j = i - 1; j >= 1 && temp[i].y - temp[j].y <= res; -- j)
res = min(res, get_dist(temp[i], temp[j]));
return res;
}
八、扫描线
扫描线是一种思想，有一定的局限性，体现在它只能适用于图形被扫描线分割后，每个部分的面积比较容易求的情况。
矩形的面积并事实上是一道线段树的题目，就不在这里放代码了。
三角形面积并：
inline bool on_segment(Point p, Point a, Point b) {
return sign((p - a) & (p - b)) <= 0;
}
inline Point get_line_intersection(Point p, Point v, Point q, Point w) {
if (sign(v * w) == 0) return Point(inf, inf);
Point u = p - q;
double t = (u * w) / (w * v);
Point res = p + v * t;
if (!on_segment(res, p, v + p) || !on_segment(res, q, w + q)) return Point(inf, inf);
return res;
}
inline double calc(double X, bool side) {
int cnt = 0;
for (int i = 1; i <= n; ++i) {
Tri now = a[i];
if (dcmp(now.v[0].x, X) > 0 || dcmp(now.v[2].x, X) < 0) continue;
if (!dcmp(now.v[0].x, X) && !dcmp(now.v[1].x, X)) {
if (!side)
np[++cnt] = Point(now.v[0].y, now.v[1].y);
continue;
}
if (!dcmp(now.v[1].x, X) && !dcmp(now.v[2].x, X)) {
if (side)
np[++cnt] = Point(now.v[1].y, now.v[2].y);
continue;
}
double u[3]; int idx = 0;
for (int k = 0; k < 3; ++k) {
Point res = get_line_intersection(now.v[k], now.v[(k + 1) % 3] - now.v[k], Point(X, -inf), Point(0, 2 * inf));
if (res.x == inf || res.y == inf) continue;
u[idx++] = res.y;
}
sort(u, u + idx);
np[++cnt] = Point(u[0], u[idx - 1]);
}
for (int i = 1; i <= cnt; ++i)
if (np[i].x > np[i].y) swap(np[i].x, np[i].y);
sort(np + 1, np + cnt + 1);
double st = -1e9, ed = -1e9, res = 0;
for (int i = 1; i <= cnt; ++i)
if (np[i].x <= ed) ed = max(ed, np[i].y);
else {
res += ed - st;
st = np[i].x;
ed = np[i].y;
}
return res + ed - st;
}
int main() {
n = read();
for (int i = 1; i <= n; ++i) {
for (int k = 0; k < 3; ++k) {
a[i].v[k].input();
vec.push_back(a[i].v[k].x);
}
sort(a[i].v, a[i].v + 3);
}
for (int i = 1; i <= n; ++i)
for (int j = i + 1; j <= n; ++j)
for (int k = 0; k < 3; ++k)
for (int l = 0; l < 3; ++l) {
Point now = get_line_intersection(a[i].v[k], a[i].v[(k + 1) % 3] - a[i].v[k], a[j].v[l], a[j].v[(l + 1) % 3] - a[j].v[l]);
if (now.x == inf || now.y == inf) continue;
vec.push_back(now.x);
}
sort(vec.begin(), vec.end());
double res = 0.0;
for (int i = 0; i + 1 < (int)vec.size(); ++i)
if (dcmp(vec[i + 1], vec[i]))
res += (vec[i + 1] - vec[i]) * (calc(vec[i], 0) + calc(vec[i + 1], 1));
res /= 2;
printf("%.2lf\n", res);
return 0;
}
