网站域名的密码,WordPress文章可视化php代码,手机投资app平台,wordpress lms缘起众所周知#xff0c;二维凸包可以使用 Graham 扫描 内解决.所以本文来学习一下三维空间中凸包的一种直观算法——增量算法(increment algorithm)分析有一条叫 Willy 的苹果虫一直快乐的居住在一个苹果中#xff0c;直到有一天有一只仓鼠想吃这个苹果#xff0c;Willy 自… 缘起众所周知二维凸包可以使用 Graham 扫描 内解决.所以本文来学习一下三维空间中凸包的一种直观算法——增量算法(increment algorithm)分析有一条叫 Willy 的苹果虫一直快乐的居住在一个苹果中直到有一天有一只仓鼠想吃这个苹果Willy 自然不想被仓鼠吃掉所以 Willy 要开始逃亡为了简化问题我们将苹果视作一个3维的凸多面体我们将给出这个凸多面体的所有顶点的数据以及一系列的 Willy 可能位于的位置坐标. 请你设计一个程序计算 Willy 最少需要爬行的距离才能爬到苹果的表面. 【输入】多样例. 每个测试样例开始于一个 n, 4≤n≤1,000, 表示该凸多面体的顶点个数. 然后 n 行每行三个整数x y z , |x|,|y|,|z| 皆 10000, 整个苹果是这 n 个点的三维凸包. 这里保证输入不会出现四点共面的情况.然后一个 q, 1≤q≤100,000, 表示q次询问然后 q 行每行三个整数 x y z , |x|, |y|, |z| 10000输入保证 Willy 一定在苹果内部. 输入以 n 0 为结束并且你不需要处理 n 0 这组数据【输出】对每次询问输出答案. 精确到四位小数. 【样例输入】60 0 0100 0 00 100 00 0 10020 20 2030 20 1041 1 130 30 357 8 990 2 20【样例输出】1.00002.88687.00002.0000【限制】20s 64MB因为保证了没有四点共面所以该凸多面体的每个面都恰好是一个三角形. 本题的思路是显然的——首先计算出三维凸包然后计算虫子到凸包的各个三角面的距离然后这些距离取最小就是答案. 计算点到面的距离是很简单的. 只需要使用平行六面体的体积除以平行四边形底面的面积即可. 前者通过混合积即可得到. 后者通过叉积即可得到. 所以本题的关键是计算三维凸包.首先我们知道计算二维凸包的算法是 Graham 扫描. 它其实本质上讲是一种增量算法. 什么是增量算法? 这涉及到一些算法的基本思想. 这里清晰的阐述一下.首先我们知道 分治 就是将大问题不断的分解为(规模不等的)小的问题. 然后自底向上(bottom up)将小问题的答案合并为大问题的答案. 这种思想运用在 dp 和 递归中. 而增量法本质上是第一数学归纳法. 即假设问题规模为 n - 1 的时候已经计算得到了答案现在要解决问题规模为 n. 增量算法是自顶向下的(top-down). 打个不恰当的比方增量算法有点像贪心.分治的复杂度公式为而增量的复杂度公式为根据上面的学习我们就不难理解为什么说 Graham 扫描是一种增量算法了.那么放到三维怎么使用增量算法求三维凸包呢? 其实和 Graham 扫描是一样的. 就是伊始选定四个不共面的点组成初始的四面体这是待求解的凸包的初始状态. 然后不断的一个一个的往点集中加入点与此过程中不断的修改(或者说维护)凸包 (下面简记 CH Convex Hull) 的样子. 直至成功加入最后的点则凸包就构建完毕了. 那么加入点的过程中自然会遇到两种情况加入的点在 CH 内部那么这种点显然是不需要考虑的. 直接continue考察下一个点就行了.加入的点不在 CH 内部那么这种点是要考虑的. 还是画图举例比较好. 假定当前 CH 就是一个四面体 ABCD. 然后我们考虑加入一个点 P下面考虑这种加入导致 CH 的变化.那么显然我们要做的就是删除 BCD 这个三角面然后新增 PCD、PDB、PBC 三个三角面, 就像下图这样即凸包从原先的四面体变成了现在的六面体.那么我们想想看为什么要删除 BCD 而不需要删除 ABC、ADC、ADB 三个面呢? 原因很简单因为如果你想象一下 P 处放一个单点光源那么BCD 面将被照耀(白天)但是 ABC、ADC、ADB 三个面将处于一片漆黑因为不能被照到.所以我们必须想办法刻画能被照耀这件事情.自然的聪明如你想到了外法向量这个有力工具. 强调一遍这是自然的因为外法向量本身就是曲面定向的工具.如上图所示n1是三角面BCD 的外法向量、n2 是三角面ADC的外法向量、n3是三角面ADB的外法向量、n4是三角面ABC的外法向量. P的加入之所以引发了 BCD 面被删除而 ADB、ADC、ABC 三个三角面不需要被删除的根本原因在于 P 在 BCD 的侧刚好是 n1 指向的方向. 用数学的话来讲就是而 P 在另外三个面的侧和其相应的外法向量是反向的. 注意CH 的所有面中只要有一个面出现了 BCD 这种情况我们立马就可以肯定 P 一定在 CH 的外部.于是我们就知道了维护一个三角面的外法向是极为重要的.所以你可以开一个四维数组来保存每个三角面的外法向量. 但是这种存储并不够优美和内蕴(implicit). 更为内蕴的做法是——用混合积. (1)式的充要条件其实是因此我们只需要记录每个三角面的三个顶点的顺序即可. 例如对于三角面 BCD, 它的顺序是 D -- B -- C, 也就是我们如果用一个数据结构 Surface 来刻画三角面的话就长下面的样子class Surface { int a, b, c, flag; Surface(int a, int b, int c, int flag 0):a(a), b(b), c(c){}}则 a 就代表 D 点b 就代表 B 点 c 就代表 C 点(因为这样的话 就是三角面 BCD 的外法向量)flag 0 或者1因为 CH 在整个凸包构建过程中不断地变化所以有一些面会曾经被加入后来又被删除(例如 BCD 面)所以 flag 1 表示该面用于构成当前 CH flag 0 表示该面不用于构成当前 CH. 然后 P 在该面的外部(也就是(1))的充要条件就化为如下混合积大于0于是我们就很清楚 CH 的整个构建过程了就是每次加入 P 之后遍历检查每个当前参与构建 CH 的面 S(即flag 0 的Surface)然后检查 (2) 是否成立如果成立的话就表明S需要从 CH 中移除(即置 S.flag 0). 但是要注意P的引入可能不止引起一个面被移除所以我们还需要去考虑 P 会不会引起和 S 相邻的三角面的移除. 还是画图说明方便上图中 P 的引入已经导致了 abc 三角面的移除但是 P 的引入可能不仅仅导致 abc 面的移除还可能导致 abd、bce、acf 这三个和 abc 接壤的三角面的移除. 所以我们还需要考察 p 和 这三个三角面之间是否满足 (2) 如果满足就要移除并且再去考察接壤的三角面是否需要被移除. 例如 上图中如果我们发现 abd 也要被移除的话我们就会再去考察和 abd 接壤的诸如 bde 这个三角面是否需要被移除. 如果发现不需要移除了(即(2)不被满足)则我们就要新建一张三角面这张三角面将参与到 CH 的构建中(即 CH 的样子发生了变化). 就好比图1的例子我们发现 ABC、ADC、ADB 三张三角面不需要被移除所以就新增了 PDC、PDB、PCB 三张三角面.所以我们发现了本质上移除新增的过程可以使用一个 dfs 来进行. 递归的出口发生三角面的新增.于是就可以来切代码了.//#define LOCAL#pragma GCC optimize(2)#pragma G optimize(2)#pragma warning(disable:4996)#pragma comment(linker, /STACK:1024000000,1024000000)#include #include #include #include #include #include #include #include #include #include //#include #include #include #include #include #include #include #include #include using namespace std;//#define int unsigned long long//#define int long long#define re register int#define ci const int#define ui unsigned int typedef pairint, int P;#define FE(cur) for(re h head[cur], to; ~h; h g[h].nxt)#define ilv inline void#define ili inline int#define ilc inline char#define ild inline double#define ilp inline P#define LEN(cur) (hjt[cur].r - hjt[cur].l)#define MID(cur) (hjt[cur].l hjt[cur].r 1)#define SQUARE(x) ((x) * (x))typedef vectorint::iterator vit;typedef setint::iterator sit;typedef mapint, int::iterator mit;const int inf ~0u 1;//const int inf 0x3f3f3f3f;const double PI acos(-1.0), eps 1e-8;namespace fastio{ const int BUF 1 21; char fr[BUF], fw[BUF], * pr1 fr, * pr2 fr; int pw; ilc gc() { return pr1 pr2 (pr2 (pr1 fr) fread(fr, 1, BUF, stdin), *pr2 0, pr1 pr2) ? EOF : *pr1; } ilv flush() { fwrite(fw, 1, pw, stdout); pw 0; } ilv pc(char c) { if (pw BUF) flush(); fw[pw] c; } ili read(int x){ x 0; int f 1; char c gc(); if (!~c) return EOF; while (!isdigit(c)) { if (c -) f -1; c gc(); } while (isdigit(c)) x (x 3) (x 1) (c ^ 48), c gc(); x * f; return 1; } ili read(double x){ int xx 0; double f 1.0, fraction 1.0; char c gc(); if (!~c) return EOF; while (!isdigit(c)) { if (c -) f -1.0; c gc(); } while (isdigit(c)) { xx (xx 3) (xx 1) (c ^ 48), c gc(); } x xx; if (c ^ .) { x f * xx; return 1; } c gc(); while (isdigit(c)) x (c ^ 48) * (fraction / 10), c gc(); x * f; return 1; } ilv write(int x) { if (x 0) pc(-), x -x; if (x 9) write(x / 10); pc(x % 10 48); } ilv writeln(int x) { write(x); pc(10); } ili read(char* x){ char c gc(); if (!~c) return EOF; while (!isalpha(c) !isdigit(c)) c gc(); while (isalpha(c) || isdigit(c)) *x c, c gc(); *x 0; return 1; } ili readln(char* x){ char c gc(); if (!~c) return EOF; while (c 10) c gc(); while (c 32 c 126) *x c, c gc(); *x 0; return 1; } ilv write(char* x) { while (*x) pc(*x); } ilv write(const char* x) { while (*x) pc(*x); } ilv writeln(char* x) { write(x); pc(10); } ilv writeln(const char* x) { write(x); pc(10); } ilv write(char c) { pc(c); } ilv writeln(char c) { write(c); pc(10); }} using namespace fastio;#define cp const Point#define cs const Surfaceci maxn 1005;ili dbcmp(double x){ if (fabs(x) { return 0; } return x eps ? 1 : -1;}ilv mmin(double a, double b){ if (a b) { a b; }}struct CH3D // 三维凸包{ struct Point // 点 { double x, y, z; Point(double x 0, double y 0, double z 0) :x(x), y(y), z(z) {} Point operator (cp o) const { return Point(x o.x, y o.y, z o.z); } Point operator - (cp o) const { return Point(x - o.x, y - o.y, z - o.z); } double operator * (cp o) const { return x * o.x y * o.y z * o.z; } Point operator / (cp o) const { return Point(y * o.z - z * o.y, z * o.x - x * o.z, x * o.y - y * o.x); } double magnitude() const{ return sqrt(SQUARE(x) SQUARE(y) SQUARE(z)); } Point scalar(double lambda) const { return Point(x * lambda, y * lambda, z * lambda); } } ps[maxn]; int n; struct Surface { int a, b, c, flag; Surface(int a 0, int b 0, int c 0) :a(a), b(b), c(c) { flag 1; } } surfaces[maxn 3]; int num; int g[maxn][maxn]; ild area(cp a, cp b, cp c) { return ((b - a) / (c - a)).magnitude() / 2.0; } ild area(cs sur) { return area(ps[sur.a], ps[sur.b], ps[sur.c]); } ild area() // 三维凸包的表面积{ double ans 0; if (n 3) { return area(ps[0], ps[1], ps[2]); } for (re i 0; i { ans area(ps[surfaces[i].a], ps[surfaces[i].b], ps[surfaces[i].c]); } return ans; } ild volume(cp a, cp b, cp c, cp d) { return ((b - a) / (c - a)) * (d - a) / 6.0; } ild volume(cp p, cs sur) { return volume(ps[sur.a], ps[sur.b], ps[sur.c], p); } ild volume() // 三维凸包的体积{ double ans 0; Point o; for (re i 0; i { ans volume(o, ps[surfaces[i].a], ps[surfaces[i].b], ps[surfaces[i].c]); } return ans; } ili ck(cp p, cs sur) { return dbcmp(volume(ps[sur.a], ps[sur.b], ps[sur.c], p)) 1; } ili onsurface(cp p, cs sur) { return !dbcmp(volume(ps[sur.a], ps[sur.b], ps[sur.c], p)); } ili same(cs s, cs t) { return onsurface(ps[s.a], t) onsurface(ps[s.b], t) onsurface(ps[s.c], t); } Point gg() // 计算三维凸包的重心{ Point ans, o ps[0]; double v 0, t; for (re i 0; i { if (surfaces[i].flag) { cp a ps[surfaces[i].a]; cp b ps[surfaces[i].b]; cp c ps[surfaces[i].c]; t volume(o, surfaces[i]); ans ans (a b c o).scalar(t / 4); v t; } } return ans.scalar(1 / v); } ili countSurface() // 返回三维凸包有多少张表面{ int ans 0; for (re i 0, ok; i { ok 1; for (re j 0; j { if (same(surfaces[i], surfaces[j])) { ok 0; break; } } ans ok; } return ans; } ili prework(){ if (n 4) { return 0; } int ok 1; for (re i 1; i { if (dbcmp((ps[0] - ps[i]).magnitude())) { swap(ps[1], ps[i]); ok 0; break; } } if (ok) { return 0; } ok 1; for (re i 2; i { if (dbcmp(area(ps[0], ps[1], ps[i]))) { swap(ps[2], ps[i]); ok 0; break; } } if (ok) { return 0; } ok 1; for (re i 3; i { if (dbcmp(volume(ps[0], ps[1], ps[2], ps[i]))) { swap(ps[3], ps[i]); ok 0; break; } } return !ok; } ilv construct() // 构造三维凸包{ num 0; if (!prework()) { return; } Surface add; for (re i 0; i 4; i) { Surface add((i 1) % 4, (i 2) % 4, (i 3) % 4); if (ck(ps[i], add)) { swap(add.b, add.c); } g[add.a][add.b] g[add.b][add.c] g[add.c][add.a] num; surfaces[num] add; } for (re i 4; i { for (re j 0; j { if (surfaces[j].flag ck(ps[i], surfaces[j])) { dfs(i, j); break; } } } int tmp num; num 0; for (re i 0; i { if (surfaces[i].flag) { surfaces[num] surfaces[i]; } } } void dfs(int p, int j) { surfaces[j].flag 0; rmv(p, surfaces[j].b, surfaces[j].a); rmv(p, surfaces[j].c, surfaces[j].b); rmv(p, surfaces[j].a, surfaces[j].c); } void rmv(int p, int a, int b) { int f g[a][b]; if (surfaces[f].flag) { if (ck(ps[p], surfaces[f])) { dfs(p, f); } else { Surface add(b, a, p); g[b][a] g[a][p] g[p][b] num; surfaces[num] add; } } } ild dis(cp p) // 计算p到三维凸包的每个三角面的距离中的最小值{ double ans inf; for (re i 0; i { mmin(ans, dis(p, surfaces[i])); } return ans; } ild dis(cp p, cs sur) { return fabs(volume(ps[sur.a], ps[sur.b], ps[sur.c], p)) / area(sur) * 3.0; }}ch;signed main(){#ifdef LOCAL FILE* ALLIN freopen(d:\\data.in, r, stdin); // freopen(d:\\my.out, w, stdout);#endif while (read(ch.n), ch.n) { for (re i 0; i ch.construct(); int q; read(q); CH3D::Point tp; while (q--) { read(tp.x), read(tp.y), read(tp.z); printf(%.4lf\n, ch.dis(tp)); } } flush();#ifdef LOCAL fclose(ALLIN);#endif return 0;}ac情况Accepted 5678ms 7252kB C最后指出的是代码中封装了三维凸包的体积、表面积、面的个数等等有用函数都经过题目检验都是正确无误的, 大可放心使用.温馨提示如果你喜欢本文请分享到朋友圈想要获得更多信息请关注ACM算法日常。点赞的时候请宠溺一点
