三维凸包

ACM算法日常

发布于 2020-08-27 10:27:03

1.9K0

发布于 2020-08-27 10:27:03

文章被收录于专栏：ACM算法日常

缘起

众所周知，二维凸包可以使用 Graham 扫描

O(n\log n)

内解决.

所以本文来学习一下三维空间中凸包的一种直观算法——增量算法（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）. 打个不恰当的比方，增量算法有点像贪心.

分治的复杂度公式为

T(n) = \Sigma_{i=1}^{n-1} T(i)+f(n)

而增量的复杂度公式为

T(n) = T(n-1) + f(n)

根据上面的学习，我们就不难理解为什么说 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 指向的方向. 用数学的话来讲就是

\overrightarrow{n_1} \cdot \overrightarrow{DP} \gt 0 \ \ \ \ \ (1)

而 P 在另外三个面的侧和其相应的外法向量是反向的. 注意，CH 的所有面中，只要有一个面出现了 BCD 这种情况，我们立马就可以肯定 P 一定在 CH 的外部.

于是，我们就知道了，维护一个三角面的外法向是极为重要的.

所以，你可以开一个四维数组，来保存每个三角面的外法向量. 但是这种存储并不够优美和内蕴（implicit）. 更为内蕴的做法是——用混合积. （1）式的充要条件其实是

(\overrightarrow{DB} \times \overrightarrow{DC}) \cdot \overrightarrow{DP} \gt 0

因此，我们只需要记录每个三角面的三个顶点的顺序即可. 例如对于三角面 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 点（因为这样的话，

\overrightarrow{ab}\times \overrightarrow{ac}

就是三角面 BCD 的外法向量），flag = 0 或者1，因为 CH 在整个凸包构建过程中不断地变化，所以有一些面会曾经被加入，后来又被删除（例如 BCD 面），所以 flag = 1 表示该面用于构成当前 CH， flag = 0 表示该面不用于构成当前 CH. 然后 P 在该面的外部（也就是（1））的充要条件就化为如下混合积大于0

(\overrightarrow{ab}\times \overrightarrow{ac})\cdot \overrightarrow{ap} \gt 0 \ \ \ \ \ (2)

于是，我们就很清楚 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 <stdio.h>
#include <iostream>
#include <iomanip>
#include <string>
#include <ctype.h>
#include <string.h>
#include <fstream>
#include <sstream>
#include <math.h>
#include <map>
//#include <unordered_map>
#include <algorithm>
#include <vector>
#include <stack>
#include <queue>
#include <set>
#include <time.h>
#include <stdlib.h>
#include <bitset>
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) < eps)
    {
        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 < num; 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 < num; 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 < num; 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 < num; i++) 
        {
            ok = 1;
            for (re j = 0; j < i; 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 < n; 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 < n; 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 < n; 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 < n; i++) 
        {
            for (re j = 0; j < num; j++) 
            {
                if (surfaces[j].flag && ck(ps[i], surfaces[j])) 
                {
                    dfs(i, j); 
                    break;
                }
            }
        }
        int tmp = num; num = 0;
        for (re i = 0; i < tmp; 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 < num; 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.n; i++) read(ch.ps[i].x), read(ch.ps[i].y), read(ch.ps[i].z);
        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情况