【学习笔记】立体计算几何

【基本知识】

向量运算：
模长：$len=\sqrt{x^2+y^2+z^2}$
加减：对应坐标加减，结果为一个向量；
点积：依旧定义为$\vec{a}$到$\vec{b}$的投影，计算方式为将对应坐标相乘后相加，结果为一个值。如
$(x_1,y_1,z_1)\cdot (x_2,y_2,z_2)=x_1{x}_2+y_1{y}_2+z_1{z}_2$。
叉积：**叉积的结果是一个向量。**方向可用“右手定则”~~（瞎编的名字）~~确定：食指指向一个向量，中指指向另一个，大拇指与中指、食指所在平面垂直。则大拇指的方向即为叉积结果的方向。

坐标计算：
$(x_1,y_1,z_1)\times (x_2,y_2,z_2)=(y_1{z}_2-y_2{z}_1,z_1{x}_2-x_1{z}_2,x_1{y}_2-x_2{y}_1)$
平面的法向量：见向量叉积的定义。

---

基本操作：
·判断点D是否在平面ABC内： 求平面ABC的法向量$\vec{n}$，若$\vec{n}\cdot \overrightarrow{AD}=0$，则$D\in 平面ABC$。
点D到平面ABC距离： 由向量点积的定义：$$dis=\frac{\vec{n}\cdot \overrightarrow{AD}}{|\vec{n}|}$$

---

算法（持续更新）：
求凸包。
step1：扰动。
将每一个点的坐标加/减一个足够小的随机值，使其不会引起误差，并能保证不会出现四点共面的情况。
step2：记录平面
由于最多三点在一个平面上，故通过记录平面上的三个点记录一个平面。
step3：构造
假如我们构造出了一个凸包。现在尝试扩展。
我们从待选点中拿一个出来，记为$P$。从$P$向凸包作射线。
不难发现，这若干条射线一定与凸包交于若干个棱上。

将可见面删去，并将$P$加入凸包。这样，我们就完成了一次扩展。

为保证正确性，需要对每一个平面判断是否可见，故时间复杂度为$O(n^2)$。
判断某个点是否可见：用平面的法向量与P到平面上的任意一点做点积。若结果$>0$（即夹角为锐角），则该平面可见；反之不可见。

#include<bits/stdc++.h>
using namespace std;
const int mn = 2005;
const double eps = 1e-10;
double Rand() {return rand() / (double) RAND_MAX;}
double reps() {return (Rand() - 0.5) * eps;}
struct point{
    double x, y, z;
    void shake() {x += reps(), y += reps(), z += reps();}
    double len() {return sqrt(x * x + y * y + z * z);}
    point operator -(const point b) const {return (point) {x - b.x, y - b.y, z - b.z};}
    point operator *(const point b) const {return (point) {y * b.z - b.y * z, z * b.x - b.z * x, x * b.y - b.x * y};}
    double operator&(const point b) const {return x * b.x + y * b.y + z * b.z;}
}p[mn];
struct surf{
    int v[3];
    point normal() {return (p[v[1]] - p[v[0]]) * (p[v[2]] - p[v[0]]);}
    double S() {return normal().len() / 2.0;}
}f[mn], c[mn];
int n, cnt, vis[mn][mn];
int check(surf f, point pp) {return ((pp - p[f.v[0]]) & f.normal()) > 0;}
inline void convey()
{
    f[++cnt] = (surf) {1, 2, 3}, f[++cnt] = (surf) {3, 2, 1};
    for(int i = 4, tmp = 0; i <= n; i++)
    {
        for(int j = 1, v; j <= cnt; j++)
        {
            if(!(v = check(f[j], p[i]))) c[++tmp] = f[j];
            for(int k = 0; k < 3; k++) vis[f[j].v[k]][f[j].v[(k + 1) % 3]] = v;
        }
        for(int j = 1; j <= cnt; j++)
            for(int k = 0; k < 3; k++)
            {
                int x = f[j].v[k], y = f[j].v[(k + 1) % 3];
                if(vis[x][y] && !vis[y][x]) c[++tmp] = (surf) {x, y, i};
            }
        for(int j = 1; j <= tmp; j++) f[j] = c[j];
        cnt = tmp, tmp = 0;
    }
}
int main()
{
    scanf("%d", &n);
    for(int i = 1; i <= n; i++)
        scanf("%lf%lf%lf", &p[i].x, &p[i].y, &p[i].z), p[i].shake();
    convey(); double ans = 0;
    for(int i = 1; i <= cnt; i++) ans += f[i].S();
    printf("%.6f\n", ans);
}