计算几何

旋转卡壳

其实旋转卡壳就是枚举每一条凸包上的边，然后找一找离这条边最远的点是哪个。确定一条边，找那个点用的就是双指针法。因此，旋转卡壳大多是枚举的边，在边上找特征。

2938. 周游世界

题意：给定一个二维平面，平面上有 N 个点。每个点的位置可由一对整数坐标 (x,y) 来表示，不同的点位置不同。请你求出平面上距离最远的点对之间的距离是多少。

#include<iostream>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;#define x first
#define y secondtypedef pair<int, int> P;
const int maxn = 50010;
int N;
P q[maxn];
int stk[maxn], top;
bool used[maxn];P operator-(P a, P b) {return { a.x - b.x, a.y - b.y };
}int operator*(P a, P b) {return a.x * b.y - a.y * b.x;
}int area(P a, P b, P c) {return (b - a) * (c - a);
}int get_dist(P a, P b) {int dx = a.x - b.x;int dy = a.y - b.y;return dx * dx + dy * dy;
}void get_convex() {sort(q, q + N);for (int i = 0; i < N; i++) {while (top >= 2 && area(q[stk[top - 2]], q[stk[top - 1]], q[i]) <= 0) {if (area(q[stk[top - 2]], q[stk[top - 1]], q[i]) < 0) {used[stk[--top]] = false;}else top--;}stk[top++] = i;used[i] = true;}used[0] = false;for (int i = N - 1; i >= 0; i--) {if (used[i]) continue;while (top >= 2 && area(q[stk[top - 2]], q[stk[top - 1]], q[i]) <= 0)top--;stk[top++] = i;}//因为0号点在起点和终点都插入了一次，因此删掉一个。top--;
}int rotating_calipers() {if (top <= 2) return get_dist(q[0], q[N - 1]);int res = 0;for (int i = 0, j = 2; i < top; i++) {//这里 i + 1 不用取模，因为根据上面那个凸包的函数，此时的 stk[top] 就是第一个点的下标auto d = q[stk[i]], e = q[stk[i + 1]];while (area(d, e, q[stk[j]]) < area(d, e, q[stk[j + 1]])) j = (j + 1) % top;res = max(res, max(get_dist(q[stk[j]], d), get_dist(q[stk[j]], e)));}return res;
}int main() {scanf("%d", &N);for (int i = 0; i < N; i++) scanf("%d%d", &q[i].x, &q[i].y);get_convex();printf("%d\n", rotating_calipers());return 0;
}

2142. 最小矩形覆盖

已知平面上不共线的一组点的坐标，求覆盖这组点的面积最小的矩形。输出矩形的面积和四个顶点的坐标。
上面那个证明出来，矩形至少有一条边和凸包的一条边重合。那么，我们其实就是枚举凸包的每一条边，然后找到离这条边最远的边，即投影最大的。然后，找两边的点时，就是找投影最小的。这样矩形的三个点就找到了，最后一个点可以根据矩形的边的关系算出来。

#include<cstring>
#include<algorithm>
#include<iostream>
#include<cmath>
using namespace std;#define x first
#define y secondtypedef pair<double, double> P;
const int maxn = 50010;
const double eps = 1e-12, PI = acos(-1), INF = 1e20;int N;
P q[maxn], ans[10];
double min_area = INF;int stk[maxn], top;
bool used[maxn];int sign(double x) {if (fabs(x) < eps) return 0;if (x < 0) return -1;return 1;
}
int dcmp(double x, double y) {if (fabs(x - y) < eps) return 0;if (x < y) return -1;return 1;
}
P operator+(P a, P b) {return { a.x + b.x, a.y + b.y };
}
P operator-(P a, P b) {return{ a.x - b.x, a.y - b.y };
}
double operator*(P a, P b) {return a.x * b.y - a.y * b.x;
}
P operator*(P a, double t) {return { a.x * t, a.y * t };
}
P operator/(P a, double t) {return { a.x / t, a.y / t };
}
double operator&(P a, P b) {return a.x * b.x + a.y * b.y;
}
double area(P a, P b, P c) {return (b - a) * (c - a);
}
double get_len(P a) {return sqrt(a & a);
}
double project(P a, P b, P c) {//AC 在 AB 上面的投影。return ((b - a) & (c - a)) / get_len(b - a);
}P unit(P a) {//单位向量return a / get_len(a);
}
P rotate(P a, double b) {return { a.x * cos(b) + a.y * sin(b), -a.x * sin(b) + a.y * cos(b) };
}void get_convex() {sort(q, q + N);for (int i = 0; i < N; i++) {while (top >= 2 && area(q[stk[top - 2]], q[stk[top - 1]], q[i]) <= 0) {if (area(q[stk[top - 2]], q[stk[top - 1]], q[i]) < 0)used[stk[--top]] = false;else top--;}stk[top++] = i;used[i] = true;}used[0] = false;for (int i = N - 1; i >= 0; i--) {if (used[i]) continue;while (top >= 2 && area(q[stk[top - 2]], q[stk[top - 1]], q[i]) <= 0)top--;stk[top++] = i;}top--;
}void rotating_calipers() {//由题意知不用判断是否共线for (int i = 0, a = 2, b = 2, c = 2; i < top; i++) {auto d = q[stk[i]], e = q[stk[i + 1]];while (dcmp(area(d, e, q[stk[a]]), area(d, e, q[stk[a + 1]])) < 0)a = (a + 1) % top;while (dcmp(project(d, e, q[stk[b]]), project(d, e, q[stk[b + 1]])) < 0)b = (b + 1) % top;//c这个点一开始必须从a开始往左走，这样才能保证c是在往极小值点移动。if (!i) c = a;//注意这里是 > 0 了，因为在左边投影是负的。//一定要小心这个 > 0 写在括号外面，不然会陷入死循环。while (dcmp(project(d, e, q[stk[c]]), project(d, e, q[stk[c + 1]])) > 0) c = (c + 1) % top;auto x = q[stk[a]], y = q[stk[b]], z = q[stk[c]];//矩形的高和宽auto h = area(d, e, x) / get_len(e - d);auto w = ((y - z) & (e - d)) / get_len(e - d);if (h * w < min_area) {min_area = h * w;ans[0] = d + unit(e - d) * project(d, e, y);ans[3] = d + unit(e - d) * project(d, e, z);auto u = unit(rotate(e - d, -PI / 2));ans[1] = ans[0] + u * h;ans[2] = ans[3] + u * h;}}
}int main() {scanf("%d", &N);for (int i = 0; i < N; i++) scanf("%lf%lf", &q[i].x, &q[i].y);get_convex();rotating_calipers();printf("%.5f\n", min_area);int k = 0;for (int i = 1; i < 4; i++) {if (dcmp(ans[i].y, ans[k].y) < 0 ||!dcmp(ans[i].y, ans[k].y) && dcmp(ans[i].x, ans[k].x) < 0)k = i;}for (int i = 0; i < 4; i++, k++) {auto x = ans[k % 4].x, y = ans[k % 4].y;//直接输出的话，会有-0.00000if (!sign(x)) x = 0;if (!sign(y)) y = 0;printf("%.5f %.5f\n", x, y);}return 0;
}

三角剖分

3034. 望远镜

讨论三角形与圆的交点

#include<iostream>
#include<cstring>
#include<algorithm>
#include<cmath>#define x first
#define y secondusing namespace std;typedef pair<double, double> P;
const int maxn = 55;
const double eps = 1e-8;
//千万别写成 cos(-1)
const double PI = acos(-1);double R;
int N;
P q[maxn], o;int sign(double x) {if (fabs(x) < eps) return 0;if (x < 0) return -1;return 1;
}//小心，这里 < eps 千万别写成 < 0.
int dcmp(double x, double y) {if (fabs(x - y) < eps) return 0;if (x < y) return -1;return 1;
}P operator+(P a, P b) {return { a.x + b.x, a.y + b.y };
}
P operator-(P a, P b) {return { a.x - b.x, a.y - b.y };
}
double operator*(P a, P b) {return a.x * b.y - a.y * b.x;
}
P operator*(P a, double t) {return { a.x * t, a.y * t };
}
P operator/(P a, double t) {return { a.x / t, a.y / t };
}
double operator&(P a, P b) {return a.x * b.x + a.y * b.y;
}double area(P a, P b, P c) {return (b - a) * (c - a);
}
double get_len(P a) {return sqrt(a & a);
}
double get_dist(P a, P b) {return get_len(b - a);
}
double project(P a, P b, P c) {//AC 在 AB 上面的投影return ((b - a) & (c - a)) / get_len(b - a);
}
P rotate(P a, double b){return { a.x * cos(b) + a.y * sin(b), -a.x * sin(b) + a.y * cos(b) };
}
P unit(P a) {return a / get_len(a);
}
bool on_segment(P p, P a, P b) {return !sign((p - a) * (p - b)) && ((p - a) & (p - b)) <= 0;
}
P get_line_intersection(P p, P v, P q, P w) {P u = p - q;double t = (w * u) / (v * w);return p + v * t;
}double get_sector(P a, P b) {//a 和 b 形成的扇形的面积，很容易看出圆心角一定小于 PI.auto angle = acos((a & b) / get_len(a) / get_len(b));//如果 oab 形成有向面积是负的，那么返回的扇形的面积也要是负的。if (sign(a * b) < 0) angle = -angle;return R * R * angle / 2;
}double get_circle_line_intersection(P a, P b, P& pa, P& pb) {auto e = get_line_intersection(a, b - a, o, rotate(b - a, PI / 2));//mind 表示 o 到 线段 ab 的距离auto mind = get_dist(o, e);//判断 e 是否在 ab上面if (!on_segment(e, a, b)) mind = min(get_dist(o, a), get_dist(o, b));//线段ab与圆无交点if (dcmp(R, mind) <= 0) return mind;//求垂足到线段与圆交点的距离auto len = sqrt(R * R - get_dist(o, e) * get_dist(o, e));pa = e + unit(a - b) * len;pb = e + unit(b - a) * len;return mind;
}double get_circle_triangle_area(P a, P b) {//一定要注意，是按照从 a 到 b 逆时针的顺序计算叉积和扇形面积。//a 和 b 共线if (!sign(a * b)) return 0;auto da = get_dist(o, a), db = get_dist(o, b);//第一种情况if (dcmp(da, R) <= 0 && dcmp(db, R) <= 0) return a * b / 2;P pa, pb;  //与圆的交点//mind 表示 o 到线段 ab 的距离auto mind = get_circle_line_intersection(a, b, pa, pb);//第二种情况if (dcmp(R, mind) <= 0) return get_sector(a, b);//第三种情况if (dcmp(da, R) <= 0) return a * pb / 2 + get_sector(pb, b);//第四种情况if (dcmp(db, R) <= 0) return pa * b / 2 + get_sector(a, pa);//第五种情况return get_sector(a, pa) + get_sector(pb, b) + pa * pb / 2;
}double work() {//题目没说是顺时针还是逆时针。//如果给的是逆时针，那么求得面积是正的。//如果给的是顺时针，那么求的面积是负的。double res = 0;for (int i = 0; i < N; i++) {res += get_circle_triangle_area(q[i], q[(i + 1) % N]);}return fabs(res);
}int main() {while (cin >> R >> N) {for (int i = 0; i < N; i++) scanf("%lf%lf", &q[i].x, &q[i].y);printf("%.2f\n", work());}return 0;
}

扫描线

主要是用来求面积或者周长

3068. 扫描线

在二维平面中给定 n 个两条边分别与 x 轴和 y 轴平行的矩形，请你求出它们的面积并。
先把举行切割成一片一片的，再把每一片的矩形做一个区间和并（排序），复杂度 O ( n 2 l o g n ) O(n^2logn) O(n2logn).

#include<bits/stdc+
using namespace std;#define x first
#define y secondconst int maxn = 1010;
typedef long long ll;
typedef pair<ll, ll> P;
P l[maxn], r[maxn];
vector<ll> xs;
int N;
P q[maxn];  //存储高度的区间ll range_area(ll a, ll b) {int cnt = 0;for (int i = 0; i < N; i++) {if (l[i].x <= a && b <= r[i].x) {q[cnt++] = { l[i].y, r[i].y };}}if (!cnt) return 0;sort(q, q + cnt);ll st = q[0].first, ed = q[0].second;ll res = 0;for (int i = 1; i < cnt; i++) {if (q[i].first > ed) {res += ed - st;st = q[i].first, ed = q[i].second;}else ed = max(ed, q[i].second);}res += ed - st;return res * (b - a);
}int main() {scanf("%d", &N);for (int i = 0; i < N; i++) {scanf("%lld%lld%lld%lld", &l[i].x, &l[i].y, &r[i].x, &r[i].y);xs.push_back(l[i].x), xs.push_back(r[i].x);}sort(xs.begin(), xs.end());ll ans = 0;for (int i = 0; i + 1 < xs.size(); i++) {if (xs[i] != xs[i + 1]) {ans += range_area(xs[i], xs[i + 1]);}}printf("%lld\n", ans);return 0;
}

2801. 三角形面积并

#include<iostream>
#include<cstring>
#include<algorithm>
#include<vector>
#include<cmath>
using namespace std;#define x first
#define y secondtypedef pair<double, double> P;
const int maxn = 110;
const double eps = 1e-8, INF = 1e6;
int N;
P tr[maxn][3];
P q[maxn];
vector<double> xs;int sign(double x) {if (fabs(x) < eps) return 0;if (x < 0) return -1;return 1;
}int dcmp(double x, double y) {if (fabs(x - y) < eps) return 0;if (x < y) return -1;return 1;
}P operator+(P a, P b) {return { a.x + b.x, a.y + b.y };
}P operator-(P a, P b) {return { a.x - b.x, a.y - b.y };
}P operator*(P a, double t) {return { a.x * t, a.y * t };
}double operator*(P a, P b) {return { a.x * b.y - b.x * a.y };
}double operator&(P a, P b) {return a.x * b.x + a.y * b.y;
}
//这里的 on_segment 需要和下面的 get_line_intersection，才能求线段交点。
bool on_segment(P p, P a, P b) {return sign((p - a) & (p - b)) <= 0;
}P get_line_intersection(P p, P v, P q, P w) {//注: p, q 是线段起点，v, w 是方向向量。要是想得到线段终点需要 p + v, q + w.//两条线段平行。if (!sign(v * w)) return { INF, INF };auto u = p - q;auto t = w * u / (v * w);auto o = p + v * t;//需要判断 o 是否在两条线段上。if (!on_segment(o, p, p + v) || !on_segment(o, q, q + w)) {return { INF, INF };}return o;
}double line_area(double a, int side) {int cnt = 0;for (int i = 0; i < N; i++) {auto t = tr[i];//只有 a 在 t[0].x 到 t[2].x 之间，三角形才与 x = a 有交点。 if (dcmp(t[0].x, a) > 0 || dcmp(t[2].x, a) < 0) continue;//需要特判三角形的一条边与 x = a 重合if (!dcmp(t[0].x, a) && !dcmp(t[1].x, a)) {if (side) q[cnt++] = { t[0].y, t[1].y };}else if (!dcmp(t[1].x, a) && !dcmp(t[2].x, a)) {if (!side) q[cnt++] = { t[1].y, t[2].y };}else {//一般情况，三角形与 x = a 要么是两个交点，要么是三个交点（顶点在x = a上）double d[3];int u = 0;for (int j = 0; j < 3; j++) {auto o = get_line_intersection(t[j], t[(j + 1) % 3] - t[j], { a, -INF }, { 0, INF * 2 });if(dcmp(INF, o.x)) d[u++] = o.y;}if (u) {sort(d, d + u);q[cnt++] = { d[0], d[u - 1] };}}}if (!cnt) return 0;for (int i = 0; i < cnt; i++) {if (q[i].first > q[i].second)swap(q[i].first, q[i].second);}sort(q, q + cnt);double res = 0, st = q[0].x, ed = q[0].y;for (int i = 1; i < cnt; i++) {if (q[i].x <= ed) ed = max(q[i].y, ed);else {res += ed - st;st = q[i].x, ed = q[i].y;}}res += ed - st;return res;
}double range_area(double a, double b) {return (line_area(a, 1) + line_area(b, 0)) * (b - a) / 2;
}int main() {scanf("%d", &N);for (int i = 0; i < N; i++) {for (int j = 0; j < 3; j++) {scanf("%lf%lf", &tr[i][j].x, &tr[i][j].y);xs.push_back(tr[i][j].x);}sort(tr[i], tr[i] + 3);}for (int i = 0; i < N; i++) {for (int j = i + 1; j < N; j++) {for (int x = 0; x < 3; x++) {for (int y = 0; y < 3; y++) {auto o = get_line_intersection(tr[i][x], tr[i][(x + 1) % 3] - tr[i][x],tr[j][y], tr[j][(y + 1) % 3] - tr[j][y]);if(dcmp(INF, o.x)) xs.push_back(o.x);}}}}sort(xs.begin(), xs.end());double res = 0;for (int i = 0; i + 1 < xs.size(); i++) {if (dcmp(xs[i], xs[i + 1])) {res += range_area(xs[i], xs[i + 1]);}}printf("%.2f\n", res);return 0;
}

自适应辛普森积分

当图形过于复杂时，可以用这个自适应辛普森积分求面积。
辛普森积分法是用二次函数来模拟这个函数，是在数值分析上，根据经验发现效果非常好的一个方法。
那么，可以这样近似

∫ a b f ( x ) d x ≈ b − a 6 ( f ( a ) + 4 f ( a + b 2 ) + f ( b ) ) \int_a^bf(x)dx \approx \frac{b-a}{6}\biggl(f(a) + 4f(\frac{a+b}{2})+f(b)\biggr) ∫abf(x)dx≈6b−a(f(a)+4f(2a+b)+f(b))

如果原函数本身的次数不超过2，那么可以求出很精确的值。

3074. 自适应辛普森积分

题意：给定两个整数 a,b，请计算如下积分：

∫ a b S i n ( x ) x d x \int_a^b \frac{Sin(x)}{x}dx ∫abxSin(x)dx

#include<cstring>
#include<algorithm>
#include<iostream>
#include<cmath>
using namespace std;const double eps = 1e-12;double f(double x) {return sin(x) / x;
}double simpson(double l, double r) {auto mid = (l + r) / 2;return (r - l) * (f(l) + 4 * f(mid) + f(r)) / 6;
}double asr(double l, double r, double s) {auto mid = (l + r) / 2;auto left = simpson(l, mid), right = simpson(mid, r);if (fabs(left + right - s) < eps) return left + right;return asr(l, mid, left) + asr(mid, r, right);
}int main() {double l, r;scanf("%lf%lf", &l, &r);printf("%.6f\n", asr(l, r, simpson(l, r)));return 0;
}

3069. 圆的面积并

#include<cstring>
#include<algorithm>
#include<iostream>
#include<cmath>#define x first
#define y secondusing namespace std;const double eps = 1e-8;
typedef pair<double, double> P;
const int maxn = 1010;int N;
struct Circle {P o;double R;
}c[maxn];P q[maxn];int dcmp(double x, double y) {if (fabs(x - y) < eps) return 0;if (x < y) return -1;return 1;
}double f(double x) {int cnt = 0;for (int i = 0; i < N; i++) {auto X = fabs(x - c[i].o.x), R = c[i].R;if (dcmp(X, R) < 0) {auto Y = sqrt(R * R - X * X);q[cnt++] = { c[i].o.y - Y, c[i].o.y + Y };}}if (!cnt) return 0;sort(q, q + cnt);double res = 0, st = q[0].x, ed = q[0].y;for (int i = 1; i < cnt; i++) {if (q[i].x <= ed) ed = max(q[i].y, ed);else {res += ed - st;st = q[i].x, ed = q[i].y;}}return res + ed - st;
}double simpson(double l, double r) {auto mid = (l + r) / 2;return (r - l) * (f(l) + 4 * f(mid) + f(r)) / 6;
}double asr(double l, double r, double s) {auto mid = (l + r) / 2;auto left = simpson(l, mid), right = simpson(mid, r);if (fabs(left + right - s) < eps) return left + right;return asr(l, mid, left) + asr(mid, r, right);
}int main() {scanf("%d", &N);for (int i = 0; i < N; i++) {scanf("%lf%lf%lf", &c[i].o.x, &c[i].o.y, &c[i].R);}double l = -2000, r = 2000;printf("%.3f\n", asr(l, r, simpson(l, r)));return 0;
}

算法模板（7）：计算几何（2）相关推荐

【蓝桥杯算法模板题--蓝桥题库Java】
PDF下载地址:点击即可文章目录 ==算法模板== 1 排序(ArrayList,sort) 题目描述输入描述输出描述输入输出样例示例 1 运行限制 2 小明的彩灯(差分) 输入输出样例示 ...
GitHub超4.4k星：程序员求职，一个算法模板就够了
来源:新智元本文约1800字,建议阅读5分钟本文为你介绍程序员求职中科学的.高效的刷题方式. [ 导读 ] 近日,GitHub上一个名为"算法模板"的项目引发热议,获得了超4. ...
hdu 2255 奔小康赚大钱--KM算法模板
题目链接:http://acm.hdu.edu.cn/showproblem.php?pid=2255 题意:有N个人跟N个房子,每个人跟房子都有一定的距离,现在要让这N个人全部回到N个房子里面去,要 ...
CSP认证201509-4 高速公路[C++题解]：强连通分量、tarjan算法模板题
题目分析来源:acwing 分析: 所求即为强连通分量的个数,然后计算每个强连通分量中点的个数,相加即可. 所谓强连通分量,它是一个子图,其中任意两点可以相互到达,并且再加一个点,就不能满足任意两点 ...
最大流ISAP算法模板
这两天学习了最大流,下面是ISAP算法模板: const int inf = 0x3fffffff; template <int N, int M> struct Isap {int to ...
九十五、二叉树的递归和非递归的遍历算法模板
@Author:Runsen 刷Leetcode,需要知道一定的算法模板,本次先总结下二叉树的递归和非递归的遍历算法模板. 二叉树的四种遍历方式,前中后加上层序遍历.对于二叉树的前中后层序遍历,每种遍 ...
Python 算法模板库，Pythonista 找工作利器
来源:Github-dashidhy https://github.com/dashidhy/algorithm-pattern-python [导语]程序员找工作,刷算法题是必不可少的一步,这里给广 ...
1470: 区间求最值（RMQ问题，ST算法模板）
1470: 区间求最值 Time Limit: 1 Sec Memory Limit: 128 MB [Submit][Status][Web Board] Description 给定一个长度为N ...
算法模板-广度优先遍历
本文转载自Maple博客的算法模板-01背包问题,转载请注明出处. 简介广度优先遍历(Breadth-First-Search,BFS),其特点便体现在广度,它是按"层"&quo ...
关于算法模板和规范的说明
这里记录的算法模板和规范仅供个人使用,可能会存在BUG.由于使用本博客模板导致的BUG及其负面影响,本人概不负责. 转载于:https://www.cnblogs.com/hzs2000/p/6682 ...

算法模板（7）：计算几何（2）

计算几何

旋转卡壳

2938. 周游世界

2142. 最小矩形覆盖

三角剖分

3034. 望远镜

扫描线

3068. 扫描线

2801. 三角形面积并

自适应辛普森积分

3074. 自适应辛普森积分

3069. 圆的面积并

算法模板（7）：计算几何（2）相关推荐

最新文章

热门文章