三维凸包 学习笔记
一些计算几何基础
三维空间内给一堆点，求其三维凸包的表面积。
其实体积也能求，最后会写。
首先要一个三维向量来表示点，方便计算。
const double eps=1e-9;
double randeps() {return (rand()/(double)RAND_MAX-0.5)*eps;}
struct vec {
double x,y,z;
vec(){x=y=z=0;}
vec(double x_,double y_,double z_){x=x_,y=y_,z=z_;}
void shake() {x+=randeps(),y+=randeps(),z+=randeps();}
double len() {return sqrt(x*x+y*y+z*z);}
vec operator + (const vec &t) {return vec(x+t.x,y+t.y,z+t.z);}
vec operator - (const vec &t) {return vec(x-t.x,y-t.y,z-t.z);}
}p[N];
double dot(const vec &a,const vec &b) {return a.x*b.x+a.y*b.y+a.z*b.z;}
vec crs(const vec &a,const vec &b) {return vec(a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x);}
几个注意事项：
- 三维向量的叉积不是数，还是向量，方向垂直于叉积的两个向量所在的平面，模长就是有向面积。
- 大家会发现有个shake函数。那个是为了防止四点共面，否则在后面求凸包的过程会出错。
关于如何记录凸包：考虑到凸包是由一堆平面，直接按照逆时针方向记录每个面的三个点（之前防止四点共面的原因）
struct face {
int v[3];
face(){v[0]=v[1]=v[2]=0;}
face(int a,int b,int c) {v[0]=a,v[1]=b,v[2]=c;}
vec normal() {return crs(p[v[1]]-p[v[0]],p[v[2]]-p[v[0]]);}//叉积，模长是面积的2倍
double area() {return normal().len()/2.0;}//面积
}cv[N],h[N];
构造凸包
(引自博客https://www.cnblogs.com/-sunshine/archive/2012/08/25/2656794.html)
大概就是这么张图。
增量法构造，每次加入一个点，把它看不见的点保留在凸包内，然后把上图加粗的棱与新增的点组成一个面，加入凸包。
- 判定能否看到：把那个face的normal（法向量）与新增的点点乘，看正负。看得到则为正。（其实就是看与法向量的夹角，判断点是否在平面上方）
- 判断粗边。代码中体现为 (vis[x][y]&&!vis[y][x]) 。那么xy就是分割线。
模拟上述过程即可，复杂度 (O(n^2))：
int cansee(face A,vec b) {return dot(b-p[A.v[0]],A.normal())>0;}
int n,th,cnt,vis[N][N];
double ans;
void convex() {
cnt=th=0;
cv[++cnt]=face(1,2,3);
cv[++cnt]=face(3,2,1);
for(rint i=4;i<=n;++i) {
for(rint j=1,v;j<=cnt;++j) {
if(!(v=cansee(cv[j],p[i])))h[++th]=cv[j];
for(rint k=0;k<3;++k)
vis[cv[j].v[k]][cv[j].v[k>1?0:k+1]]=v;
}
for(rint j=1;j<=cnt;++j)
for(rint k=0;k<3;++k) {
int x=cv[j].v[k],y=cv[j].v[k>1?0:k+1];
if(vis[x][y]&&!vis[y][x])h[++th]=face(x,y,i);
}
for(rint j=1;j<=th;++j)cv[j]=h[j];
cnt=th,th=0;
}
}
凸包构造完毕！
计算表面积的话，把每个面的normal的模长加起来除以2.
模板题代码：
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef unsigned int uint;
#define rint register int
#define pb push_back
//#define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2) EOF:*p1++)
//char buf[1<<21],*p1=buf,*p2=buf;
inline int rd() {
int x=0,f=1;char ch=getchar();
while(!isdigit(ch)) {if(ch=='-')f=-1;ch=getchar();}
while(isdigit(ch))x=x*10+(ch^48),ch=getchar();
return x*f;
}
const int N=2010;
const double eps=1e-9;
double randeps() {return (rand()/(double)RAND_MAX-0.5)*eps;}
struct vec {
double x,y,z;
vec(){x=y=z=0;}
vec(double x_,double y_,double z_){x=x_,y=y_,z=z_;}
void shake() {x+=randeps(),y+=randeps(),z+=randeps();}
double len() {return sqrt(x*x+y*y+z*z);}
vec operator + (const vec &t) {return vec(x+t.x,y+t.y,z+t.z);}
vec operator - (const vec &t) {return vec(x-t.x,y-t.y,z-t.z);}
}p[N];
double dot(const vec &a,const vec &b) {return a.x*b.x+a.y*b.y+a.z*b.z;}
vec crs(const vec &a,const vec &b) {return vec(a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x);}
struct face {
int v[3];
face(){v[0]=v[1]=v[2]=0;}
face(int a,int b,int c) {v[0]=a,v[1]=b,v[2]=c;}
vec normal() {return crs(p[v[1]]-p[v[0]],p[v[2]]-p[v[0]]);}
double area() {return normal().len()/2.0;}
}cv[N],h[N];
int cansee(face A,vec b) {return dot(b-p[A.v[0]],A.normal())>0;}
int n,th,cnt,vis[N][N];
double ans;
void convex() {
cnt=th=0;
cv[++cnt]=face(1,2,3);
cv[++cnt]=face(3,2,1);
for(rint i=4;i<=n;++i) {
for(rint j=1,v;j<=cnt;++j) {
if(!(v=cansee(cv[j],p[i])))h[++th]=cv[j];
for(rint k=0;k<3;++k)
vis[cv[j].v[k]][cv[j].v[k>1?0:k+1]]=v;
}
for(rint j=1;j<=cnt;++j)
for(rint k=0;k<3;++k) {
int x=cv[j].v[k],y=cv[j].v[k>1?0:k+1];
if(vis[x][y]&&!vis[y][x])h[++th]=face(x,y,i);
}
for(rint j=1;j<=th;++j)cv[j]=h[j];
cnt=th,th=0;
}
}
void calc() {for(rint i=1;i<=cnt;++i)ans+=cv[i].area();}
signed main() {
srand(time(0));
n=rd();
for(rint i=1;i<=n;++i)
scanf("%lf%lf%lf",&p[i].x,&p[i].y,&p[i].z),p[i].shake();
convex();calc();
printf("%.3lf
",ans);
return 0;
}
计算体积：darkbzoj #1964. hull 三维凸包
只需要在构造出凸包后更改calc函数，每次计算 p[1] 与每个面组成三棱锥的体积，相加即可。
这个用向量很方便：
double vol6(vec a,vec b,vec c,vec d) {
return dot(crs(b-a,c-a),d-a);
} //计算体积的6倍。理解：先算出底的有向面积的2倍，向量的方向变为高，然后点乘就是另一条边在高方向的投影，由体积公式(V=1/3*Sh)，这个值就是体积的6倍
void calc() {
for(rint i=1;i<=cnt;++i)
ans+=fabs(vol6(p[1],p[cv[i].v[0]],p[cv[i].v[1]],p[cv[i].v[2]]));
ans/=6;
}