缘起
众所周知，二维凸包可以使用 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 这组数据
【输出】
对每次询问，输出答案. 精确到四位小数.
【样例输入】
6
0 0 0
100 0 0
0 100 0
0 0 100
20 20 20
30 20 10
4
1 1 1
30 30 35
7 8 9
90 2 2
0
【样例输出】
1.0000
2.8868
7.0000
2.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 Surface
ci 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算法日常。
![Spring [CVE-2022-22965]漏洞处理](http://pic.xiahunao.cn/Spring [CVE-2022-22965]漏洞处理)
![[初级]Java中的switch对整型、字符型、字符串的具体实现细节](http://pic.xiahunao.cn/[初级]Java中的switch对整型、字符型、字符串的具体实现细节)
)