三維凸包學習小記

三維凸包

Tags：高級算法

---

Part 1 平面幾何基礎

出門右拐：http://www.javashuo.com/article/p-ruqcwyeu-ee.html （附計算幾何題單）

Part 2 立體幾何基礎

向量運算

加減運算

同平面向量，對應座標相加減

模長

\(|a|=\sqrt{x^2+y^2+z^2}\)

點積

兩個向量的點積仍然表示 a到b的投影×b的模長
仍然知足\(a·b=|a||b|cos<a,b>\)
座標下有\((x_1,y_1,z_1)·(x_2,y_2,z_2)=x_1x_2+y_1y_2+z_1z_2\)，對應座標相乘

叉積

兩個三維向量叉積仍然是一個三維向量（不一樣於平面向量，乘積是實數）
其模長仍然表示以這兩個三維向量做爲鄰邊的平行四邊形面積
方向符合：對於\(a*b\)，伸出右手，食指指向\(a\)，中指指向\(b\)，大拇指所對的方向爲叉積後的向量方向

如上圖，\(AC*AB=AD\)

座標表示就是：\((y_1z_2-z_1y_2，z_1x_2-x_1z_2，x_1y_2-y_1x_2)\)

平面的法向量

在平面上任選兩個向量作叉積便可

判斷點是否在平面上

平面ABC，判斷D是否在平面上。
法向量n，則若AD&n=0，點積爲零，說明D在平面上。

點到平面的距離

該點到平面上任意一點的向量 點積 平面的法向量
而後除以法向量的模長

double Dis(Node a) {Node w=Normal();return fabs((w&(a-A[v[0]]))/w.len());}

求凸包

擾動

首先對其微小擾動，避免出現四點共面的狀況

平面的記錄

擾動以後各個平面必定是一個三角形，逆時針方向記錄三個頂點表示一個面

增量構造

借用網上這篇博客的圖片方便理解
對於一個已知凸包，新增一個點P
將P視做一個點光源，向凸包作射線
能夠知道，光線的可見面和不可見面必定是由若干條棱隔開的

將光的可見面刪去，並新增由其分割棱與P構成的平面

重複此過程便可，複雜度\(O(n^2)\)，分析見Pick定理、歐拉公式和圓的反演。

代碼

洛谷模板：求三維凸包面積
先放上樣例的兩張靚照：

強烈推薦畫圖軟件Geogebra！

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cmath>
using namespace std;
const int N=2010;
const double eps=1e-9;
int n,cnt,vis[N][N];
double ans;
double Rand() {return rand()/(double)RAND_MAX;}
double reps() {return (Rand()-0.5)*eps;}
struct Node
{
    double x,y,z;
    void shake() {x+=reps();y+=reps();z+=reps();}
    double len() {return sqrt(x*x+y*y+z*z);}
    Node operator - (Node A) {return (Node){x-A.x,y-A.y,z-A.z};}
    Node operator * (Node A) {return (Node){y*A.z-z*A.y,z*A.x-x*A.z,x*A.y-y*A.x};}
    double operator & (Node A) {return x*A.x+y*A.y+z*A.z;}
}A[N];
struct Face
{
    int v[3];
    Node Normal() {return (A[v[1]]-A[v[0]])*(A[v[2]]-A[v[0]]);}
    double area() {return Normal().len()/2.0;}
}f[N],C[N];
int see(Face a,Node b) {return ((b-A[a.v[0]])&a.Normal())>0;}
void Convex_3D()
{
    f[++cnt]=(Face){1,2,3};
    f[++cnt]=(Face){3,2,1};
    for(int i=4,cc=0;i<=n;i++)
    {
        for(int j=1,v;j<=cnt;j++)
        {
            if(!(v=see(f[j],A[i]))) C[++cc]=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[++cc]=(Face){x,y,i};
            }
        for(int j=1;j<=cc;j++) f[j]=C[j];
        cnt=cc;cc=0;
    }
}
int main()
{
    cin>>n;
    for(int i=1;i<=n;i++) cin>>A[i].x>>A[i].y>>A[i].z,A[i].shake();
    Convex_3D();
    for(int i=1;i<=cnt;i++) ans+=f[i].area();
    printf("%.3f\n",ans);
}