c++ 凸包 分治算法_三维凸包
缘起
众所周知,二维凸包可以使用 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 pair<int, 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 vector<int>::iterator vit;typedef set<int>::iterator sit;typedef map<int, 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算法日常。
点赞的时候,请宠溺一点
c++ 凸包 分治算法_三维凸包相关推荐
- matlab 凸包质心算法,求多边形凸包(线性算法)--陈氏凸包算法--Computing the convex hull of a simple polygon(源码)...
陈氏凸包算法-算法参考:Computing the convex hull of a simple polygon 作者:Chern-Lin Chen 陈氏算法提供了一个线性效率求凸包的算法,本文使用 ...
- python分治算法_分治法及其python实现例子
在前面的排序算法学习中,归并排序和快速排序就是用的分治法,分治法作为三大算法之一的,有非常多的应用例子. 分治法概念 将一个复杂的问题分成两个或更多的相同或相似的子问题,再把子问题分成更小的子问题-- ...
- python分治算法_黄哥Python:分治算法(Divide-and-Conquer)
分治算法(Divide-and-Conquer) 在计算机科学中,分而治之(简称分治法)是基于多分支递归的算法设计范例.分而治之算法的工作原理是将问题递归分解为两个或多个相同或相关类型的子问题,直到这 ...
- python分治算法_算法-分治
分治算法的基本思想是将一个规模为N的问题分解为K个规模较小的子问题,这些子问题相互独立且与原问题性质相同.求出子问题的解,就可得到原问题的解,是一种分目标完成程序算法,简单的问题可用二分法完成. 1. ...
- python矩阵乘法分治算法_矩阵乘法的Strassen算法详解 --(算法导论分治法求矩阵)...
1 题目描述 2 思路分析 3 解法 4 小结 1 题目描述 请编程实现矩阵乘法,并考虑当矩阵规模较大时的优化方法. 2 思路分析 根据wikipedia上的介绍:两个矩阵的乘法仅当第一个矩阵B的列数 ...
- 基于C++实现的几何学与计算机的交叉应用(四色定理、三维凸包)
浅谈几何学与计算机的交叉应用 一.摘要 就计算机在几何学领域的应用,本文介绍了四色定理(Four Color Theorem)及其历史与机器证明,以及四色问题的求解.就几何学在计算机领域的应用,本文介 ...
- python 三维凸包_浅尝则止 - SciPy科学计算 in Python
本文节选自作者的<Python编程基础及应用>视频教程.Python编程基础及应用_哔哩哔哩 (゜-゜)つロ 干杯~-bilibiliwww.bilibili.com 21. 浅尝则止 ...
- HDUOJ 1392凸包graham算法
#include <iostream> #include <cstdio> #include <cstring> #include <cmath> #i ...
- sgu290:Defend the Milky Way(三维凸包)
题目大意: ~~~~~~给出空间中 n(1≤n≤100) ~n(1\leq n \leq 100)~个点,求出其凸包.如果有点在凸包的面或棱上,也要将其算进凸包中,将答案按字典序输出. 分 ...
最新文章
- httpclient base64 文件上传_文件上传下载
- linux如何卸载virtualbox,如何在Mac上卸载VirtualBox | MOS86
- kernel mtd 分区与UBOOT 分区的理解
- 1136 A Delayed Palindrome (20 分)
- Springboot是什么?Springboot详解!入门介绍
- 这个小伙因WannaCry勒索软件一夜成名,获得一年免费披萨
- android个推快速集成,个推用户画像产品(个像)Android集成实践
- 【路径规划】基于matlab帝国企鹅算法求解机器人栅格地图避障路径规划问题【含Matlab源码 784期】
- Python读写txt文件
- Hybird App 应用开发中5个必备知识点复习
- uniapp编译支付宝小程序图片图标显示问题
- 利用自定义注解,统计方法执行时间
- 修改电脑IP地址和MAC地址
- googletest 学习笔记
- html表格左右布局,css table布局大法,解决你大部分居中、多列等高、左右布局的问题...
- windows 程序设计 第一章
- oracle10G安装与配置
- 周期性行业是什么意思_周期性行业
- 计算机技术 高中教案,高中信息技术 计算机软件教案
- EasyIcon的图标还是挺全的
热门文章
- P4062 [Code+#1]Yazid 的新生舞会(分治做法)
- 2020牛客国庆集训派对day2 AKU NEGARAKU
- 迎开学水题狂欢赛(舞踏会[dp+三叉树],HH去散步[矩阵快速幂],排序[模拟],铁路旅行[线段树])
- P6773-[NOI2020]命运【线段树合并,树形dp】
- P3586-[POI2015]LOG【线段树】
- 5、play中的json数据处理
- JavaFX UI控件教程(二十一)之Tooltip
- 这些保护Spring Boot 应用的方法,你都用了吗?
- Spring BeanFactory 容器
- JavaScript学习总结(九)——Javascript面向(基于)对象编程