Piwry
2020-07-04 20:59:41
本文主要介绍了三维计算几何的基本运算，并取其比较基本的求凸包问题介绍了几种算法。
另外我还整理了我找到的一些资料，可以 从这（baidu 网盘，6ahh ）下载。
同二维
vect operator * (double k){ return vect(x*k, y*k, z*k); }
vect operator / (double k){ return vect(x/k, y/k, z/k); }同二维
意义也一样
vect operator + (vect v){ return vect(x+v.x, y+v.y, z+v.z); }
vect operator - (vect v){ return vect(x-v.x, y-v.y, z-v.z); }同二维
意义也和二维一样，是一个向量到另一个向量的投影长乘另一个向量的模长
满足交换律
double operator * (vect v){ return x*v.x+y*v.y+z*v.z; }vect operator / (vect v){ return vect(y*v.z-z*v.y, z*v.x-x*v.z, x*v.y-y*v.x); }返回的是一个空间向量
它的默认方向和两个输入向量满足右手定理（形如三维坐标轴：
长度是输入向量构成的平行四边形面积，
正负由第一个向量到第二个向量的逆时针夹角决定（同二维），
不满足交换律
3b1b 中（见此）有对其含义的一种解释：
将输入向量的起点置于原点，若定义运算
那么做这个运算
[1] 有向体积：设输入向量
double m(){ return sqrt(x*x+y*y+z*z); }投影到坐标面化简下式子就可以。
struct line{
vect v[2];
line(){}
line(vect uu, vect vv){ v[0] =uu, v[1] =vv; }
};和二维一样存储直线上的两个点。
该直线的方向定义为从存储的第一个点到第二个点的向量。
inline double dist(vect v, line f){ return ((f.v-f.u)/(v-f.u)).m()/(f.v-f.u).m(); }做叉积得到平行四边形的面积，除掉底长（存储直线也是存点）就是高，这个高即是距离。
struct plane{
vect v[3];
plane(){}
plane(vect uu, vect vv, vect ww){ v[0] =uu, v[1] =vv, v[2] =ww; }
vect normal(){ return (v[1]-v[0])/(v[2]-v[0]);}
};存储面上的三个点。
若定义面的 "上方" 和 "下方"，则这三个点是按从上方看到的逆时针顺序排列的
其法向量（vect normal()）即为从一个点顺次连向后两个点的向量的叉积（顺次）。这个向量竖直指向面的上方。
inline bool isabove(vect v, plane p){ return gtr((v-p.u())*p.normal()/*这是投影长度乘法向量模长，只是因为不需要知道数值只需要符号，所以没做处理*/, 0); }只需判断法向量起点到该点的这个向量在法向量上的投影符号即可。
/*带方向距离*/
inline double dist(vect v, plane p){ return (v-p.u())*p.normal()/p.normal().m(); }法向量起点到该点的这个向量在法向量上的投影长。
一个方法是判断两个面的法向量是否相等。而判断法向量相等，我们可以将它分别投影到
inline bool eq(plane p, plane q){
vect n1 =p.normal(), n2 =q.normal();
return (eq(atan2(n1.y, n1.x), atan2(n2.y, n2.x)) && eq(atan2(n1.z, n1.x), atan2(n2.z, n2.x)));
}这种方法要求两个相等的面朝向也是相等的。
不过考虑到极角为
还可以想到枚举某个面存储的三个代表点，看它们是不是都在另一个面上
bool eq(plane p, plane q){
vect a =p.v[0], b =p.v[1],c =p.v[2];
return (eq(dist(a, q), 0) && eq(dist(b, q), 0) && eq(dist(c, q), 0));
}但这种方法不会判断平面的朝向。如果实在需要，可以做通过平面外一点与两个平面的有向体积来判断
bool eq(plane p, plane q){
vect a =p.v[0], b =p.v[1],c =p.v[2];
bool flg =(eq(dist(a, q), 0) && eq(dist(b, q), 0) && eq(dist(c, q), 0));
if(flg){
vect n =p.normal()+p.u();/*可以知道这个点一定在平面外*/
flg &=!(gtr(dist(n, p), 0) && !gtr(dist(n, q), 0)/*体积符号不同则面朝向不同*/);
}
return flg;
}通过叉积的定义，我们可以很轻松地求出两个向量形成的平行四边形的面积：
vect u, v;
double S =(u/v).m();
体积问题则可以转化成底面积乘高来计算。高可以用点积的意义求出（其实也就是点到面的距离）。
举个最简单的例子，给定顶点，计算三棱锥的体积。
首先给出三棱锥体积的计算公式：
那么其代码就为：
vect A, B, C, D;/*四个顶点，具体位置见下图*/
vect n =(B-A)/(C-A);
double V =(n.m()/2)*(D*n/n.m())/3;同时，许多复杂多面体体积其实也都可以转化为多个三棱锥体积的和。
通常我们通过储存面来储存多面体。
对于多面体每个面的面积，我们可以用叉积的含义轻易算出：
plane facet;
double S =facet.normal().m()/2;/*其面积*/
在某些情况下，我们还需要存储每个面的邻接面信息。
为了方便，我们将每个顶点大于三个的面划分成若干个三角形的面；视具体情况可能划分后的面的顶点不属于原来面的顶点。
于是我们就可以这样保存面：
struct facet{
int n[3];/*neighbor，和点对应 (u->v, v->w, w->u)*/
plane p;
};对于一个多面体，我们保存其中一个面作为索引，就可以通过邻接信息读取多面体的表面。
（貌似一般都叫 "Horizon"...直译过来就是地平线...）
对于一个点和一个多面体，从该点能 "看到" 的所有面一定是一片连续的面，而这片面的边缘就是这个点在这个多面体上的地平线。
（可想象站在高处望向远处的场景...）
若要求得地平线，我们可以从某个该点能看到的面（具体来说就是这个面的平面在该点下面）开始 dfs（这个面也可以从索引面开始 dfs 求出），并在图中标记所有边缘面；
对于每个边缘面，我们将每条边缘边的两个顶点建图连边。
可以知道我们最后求出的图一定是一个环。
同时为了避免不必要的计算，我们可以在 dfs 顺便返回一个环上的一个点作为索引。
struct vect{
double x, y, z;
int id;
};
struct plane{
vect vec[3];
};
struct facet{
int n[3];
int id, vistime /*相当于一个时间戳*/;
plane p;
};
struct edge{
int netid /*连向点的 id*/, facetid /*当前这条边缘边邻接的不可看见的面的 id*/;
};
bool isabove(vect v, plane p){ /*return (点是在面上方 );*/ }
vector<facet> FAC;/*存储所有面*/
int TIME =0;/*全局时间*/
/*e1 e2 存储一个点连的两条边，resfdel 存储所有看到的面，vistime 是每个点访问的时间戳*/
/*调用前需要更新时间 (++TIME)，清空 resfdel(至于 e1 e2 不清也没关系)*/
/*最后会返回边界线上的一个点，并保存地平线信息和所有能看到的面到指定数组*/
int getHorizon(int f, vect &p, int vistime[], edge e1[], edge e2[], vector<int> &resfdel){
if(!isabove(p, FAC[f].p)/*恰好在平面上应当视为不能看见*/) return 0;/*如果返回 0 则代表这个面无法看见*/
if(FAC[f].vistime == TIME) return -1;/*如果返回 -1 则代表这个面已经见过*/
FAC[f].vistime =TIME;/*记录当前时间戳*/
FAC[f].isdel =1;
resfdel.push_back(FAC[f].id);
int ret =-2;/*如果返回 -2 则代表这个面的 dfs 树无法接触到地平线*/
/*--(思考下 dfs 一个被夹在中心的面，周围一圈的面都已搜索过 )*/
for(int i =0; i < 3; ++i){
int res =getHorizon(FAC[f].n[i], p, vistime, e1, e2, resfdel);
if(res == 0){
int pt[2];
pt[0] =FAC[f].p.vec[i].id, pt[1] =FAC[f].p.vec[(i+1)%3].id;
for(int j =0; j < 2; ++j){
/*我们无法得知这条边是沿缺口逆时针还是顺时针的*/
/*根据被访问的次数判断这是该点连的第一条还是第二条边*/
if(vistime[pt[j]] != TIME){
vistime[pt[j]] =TIME;
e1[pt[j]].netid =pt[(j+1)%2];
e1[pt[j]].facetid =FAC[f].n[i];
}
else{
e2[pt[j]].netid =pt[(j+1)%2];
e2[pt[j]].facetid =FAC[f].n[i];
}
}
ret =pt[0];
}
else if(res != -1 && res != -2) ret =res;/*保存地平线上的点*/
}
return ret;
}我们可以想到，一条边会且仅会被两个面使用；按我们储存面的方式，这条边一定是被 "逆时针用一次"，再 "顺时针用一次"。
就像下图中那样：
这个性质可以帮助我们通过统计边来获得多面体面的信息。
为了避免可能的四点共面，三点共线等问题增加问题复杂度，我们可以给点的坐标加减一个极小的值，同时也能保证一些问题中求出的答案不会偏差太大。
扰动的具体实现和应用的实际情况有关，
如果要扰动距离较近的少数点，可以这样实现（其中 eps 是判断两个数不同的最小差值，之后（应该）不再解释）：
double reps() { return (rand()%2) ? eps : -eps; }如果要扰动距离较远的多数点，可以这样实现：
double reps() {return (rand()/(double)RAND_MAX-0.5)*eps;}（较远的点对会放大精度误差，因此扰动值可以不大于 eps；点数目的增加要求扰动的值更加随机）
如果不确定问题具体是哪种情况，其实也可以综合两种实现：
double reps2() { return (rand()%2) ? eps : -eps; }
double reps() {return (rand()/(double)RAND_MAX-0.5)*eps+reps2();}我们还是给每个点一个随机的误差，但是我们将每个点的坐标显式地表示成
在运算时，首先将
这种方法对精度的影响应该会比较小；如果统计贡献时并不要求排除特殊情况（比如四点共面什么），这种方法甚至能做到无精度影响。
点扰动对一些不要求求出准确值的问题是影响不大的，例如：
但有一些问题是存在准确解的，这时候我们可能要改进扰动算法（例如保存原坐标）或选用不需要扰动点的算法，例如：
这里先介绍下较为广义的凸包定义：
对于
而对于一般的三维凸包问题，就是所有可能的由给定点集构成的面构成的单纯形的交集。
下面将介绍几种代表性的求三维凸包的方法：增量法，卷包裹法（gift wrap），QuickHull，随机增量法。
其中前两个是
我们考虑加入一个点时如何维护凸包
可以发现只有它能 "看到" 的面，或者说在该点下方的面（对于点恰好在面上，我们不算能 "看到"）会被删除。
所以初始先构造一个任意单纯形（四个顶点的多面体）。我们每加入一个点，就遍历所有面找出能看到的面，将它们删除；并对 "缺口" 的边缘边和新点连面，加入凸包。
对于边缘的标记，我们可以想到：一条边会且仅会被两个面使用；按我们储存面的方式，这条边一定是被 "逆时针用一次"，再 "顺时针用一次"。（这里 copy 自上面）
于是用一个二维数组标记，只对只使用一次的边，按使用过的方向和新点连面即可（可见代码）。
有时可能会出现几个点共面的情况。
我们考虑增量法的实现，可以发现在已有面上的点是无效的，也不会产生贡献。
不过由于我们连的面都是三角形，因此最后求得的面集合可能会把原来一个多边形面划分成多个三角形小面（同时也不保证三角形的顶点一定在边上）。注意后续运算需要考虑去重的问题。
另外，值得一提的是，如果我们通过扰动点来期望最多使三个点共面，我们初始可以只构造面。这样可以略微缩减代码量。
不扰动的话这样可能出错。考虑这样的数据：
由于扰动的值不需要太大，因此一般的精度要求用这种方法都可以满足。
（这里的实例是 求给定点集的凸包的表面积）
扰动点（只构造面）：
#include <cstdio>
#include <cstdlib>
#include <cmath>
const int MAXN =2500;
/*------------------------------Computational geometry------------------------------*/
const double eps =1e-8;
struct vect{
double x, y, z;
vect(){}
vect(double xx, double yy, double zz):x(xx), y(yy), z(zz){}
vect operator - (vect v){ return vect(x-v.x, y-v.y, z-v.z); }
vect operator / (vect v){ return vect(y*v.z-z*v.y, z*v.x-x*v.z, x*v.y-y*v.x); }
double operator * (vect v){ return x*v.x+y*v.y+z*v.z; }
double m(){ return sqrt(x*x+y*y+z*z); }
double rd(){ return (rand()%2) ? eps : -eps; }
void shake(){ x +=rd(), y +=rd(), z +=rd(); }
}pts[MAXN];
struct plane{
int v[3];/*节省空间*/
plane(){}
plane(int uu, int vv, int ww){ v[0] =uu, v[1] =vv, v[2] =ww; }
vect normal(){ return (pts[v[1]]-pts[v[0]])/(pts[v[2]]-pts[v[0]]); }
};
inline bool gtr(double a, double b){ return (a-b > eps); }
inline bool eq(double a, double b){ return (a-b > -eps && a-b < eps); }
inline double dist(vect x, vect y){ return (x-y).m(); }
inline double dist(vect v, vect f1, vect f2){ return ((f2-f1)/(v-f1)).m()/(f2-f1).m(); }
inline bool isabove(vect v, plane p){ return gtr((v-pts[p.v[0]])*p.normal(), 0); }
/*------------------------------Convex Hulls------------------------------*/
int vise[MAXN][MAXN];
plane res[MAXN<<1], del[MAXN<<1];
inline int getConvexHulls(int totp, plane facets[]){
int s[3];
s[0] =0, s[1] =1, s[2] =2;
int totf =0;
facets[totf++] =plane(s[0], s[1], s[2]);
facets[totf++] =plane(s[0], s[2], s[1]);
for(int i =3; i < totp; ++i){
/*重复的点按逻辑不会产生贡献，但考虑精度误差建议跳过这些点*/
int totr =0, totd =0;
for(int j =0; j < totf; ++j){
if(!isabove(pts[i], facets[j])) res[totr++] =facets[j];/*这里保存不需要删除的面*/
else{
del[totd++] =facets[j];
/*由于每个点只会循环一次，所以这个 i+1 就相当于一个时间戳*/
for(int k =0; k < 3; ++k)
vise[facets[j].v[k]][facets[j].v[(k+1)%3]] =i+1;
}
}
for(int j =0; j < totd; ++j){
plane f =del[j];
for(int k =0; k < 3; ++k)
/*可以证明反向边可以不判（即 vise[f.v[k]][f.v[(k+1)%3]] == i+1）*/
if(vise[f.v[(k+1)%3]][f.v[k]] != i+1)
res[totr++] =plane(f.v[k], f.v[(k+1)%3], i);
}
totf =totr;
for(int j =0; j < totr; ++j) facets[j] =res[j];
}
return totf;
}
/*------------------------------Main------------------------------*/
plane facets[MAXN<<1];
int main(){
int n; scanf("%d", &n);
for(int i =0; i < n; ++i) scanf("%lf%lf%lf", &pts[i].x, &pts[i].y, &pts[i].z), pts[i].shake();/*扰动点*/
int h =getConvexHulls(n, facets);
double area =0;
for(int i =0; i < h; ++i) area +=facets[i].normal().m()/2;
printf("%.3lf", area);
}
不扰动点（仅给出不同部分）：
/*------------------------------Convex Hulls------------------------------*/
inline int getConvexHulls(int totp, plane facets[]){
int s[4];
s[0] =0, s[1] =1, s[2] =2, s[3] =0;
while(s[2] < totp-1 && eq(dist(pts[s[2]], pts[s[0]], pts[s[1]]), 0)) ++s[2];/*确保不共线*/
while(s[3] < totp-1 && eq(dist(pts[s[3]], plane(s[0], s[1], s[2])), 0)) ++s[3];/*确保不共面*/
if(gtr(0, dist(pts[s[3]], plane(s[0], s[1], s[2])))) s[1] ^=s[2] ^=s[1] ^=s[2];/*为了保证下面构造的面向外*/
int totf =0;
facets[totf++] =plane(s[0], s[2], s[1]);
facets[totf++] =plane(s[0], s[1], s[3]);
facets[totf++] =plane(s[1], s[2], s[3]);
facets[totf++] =plane(s[2], s[0], s[3]);
for(int i =0; i < totp; ++i){
/*虽然按理来说重复点不会有贡献，但考虑有时的精度误差问题还是特判比较好*/
if(i == s[0] || i == s[1] || i == s[2] || i == s[3]) continue;
/*...*/
}
return totf;
}
/*------------------------------Main------------------------------*/
int main(){
int n; scanf("%d", &n);
for(int i =0; i < n; ++i) scanf("%lf%lf%lf", &pts[i].x, &pts[i].y, &pts[i].z);/*没有扰动*/
/*...*/
printf("%.3lf", area);
}
另外，其实我们也可以不使用 del 数组：
inline int getConvexHulls(int totp, plane facets[]){
int s[4];
/*...*/
/*不扰动点的实例*/
for(int i =0; i < totp; ++i){
if(i == s[0] || i == s[1] || i == s[2] || i == s[3]) continue;
int totr =0;
for(int j =0; j < totf; ++j){
if(isabove(pts[i], facets[j]))
for(int k =0; k < 3; ++k)
vise[facets[j].v[k]][facets[j].v[(k+1)%3]] =i+1;
else res[totr++] =facets[j];
}
totf =0;
/*直接从保留的面找地平线*/
for(int j =0; j < totr; ++j){
facets[totf++] =res[j];
for(int k =0; k < 3; ++k)
if(vise[res[j].v[(k+1)%3]][res[j].v[k]] == i+1)
facets[totf++] =plane(res[j].v[(k+1)%3], res[j].v[k], i);
}
}
return totf;
}但个人觉得这种写法在随机点集下可能会慢些，
口胡的证明（仅供参考qaq）：
随机点集下期望每一步的凸包是 “均匀" 的；也就是说，它的面是在各个方向是 ”平均“ 的，
随着凸包的成型，新加入的点可能 “看到” 的面越来越少；而即使 “看到” 的面最多，也只有半个凸包，也就是期望一半的面数。
也就是说，随机点集下，遍历看到的面的次数的最大值和遍历看不到的面的次数的最小值相等。
实际实现时按心情写就行了）
可以想到是
这种做法可以很轻松地支持在线
我们想象用一张纸慢慢地包住凸包。
具体来说，我们想象有一个平面从凸包的一个面的一条边开始向内折，并取第一个碰到的点作为新的面的顶点，构造新面；这是一次打包（wrap）。
我们再将新面新产生的两条边包含其面的信息（具体来说是边的方向，参考上文介绍的 面与边的关系）放入队列，重复同样的工作，直到没有边可以再打包即可。
比较直观的思路是求出所有点相对当前面的 "夹角"，不过这种方式难以实现。
我们再转化其含义，可以发现这实际上就是求一个新多边形面，使剩下的所有点都在新面下方（新面上方朝多面体外侧）。
/*单次打包*/
/*大部分变量名实例同上文*/
/*ne 是指当前边，边的方向保存了原来面的信息*/
pair<int, int> ne;/*存储两个顶点的 id*/
int p =-1;
for(int i =0; i < totp; ++i){
if(i == ne.first || i == ne.second) continue;/*避免作为 p 的初始值 (p == -1)*/
if(p == -1 || isabove(pts[i], plane(ne.first, ne.second, p)))
p =i;
}卷包裹法一开始需要一个属于最终凸包的边作为算法起点。
我们可以想到，如果将凸包 "拍平" 到一个平面上，其映射在平面上的点集的凸包的边一定也是原凸包的边。
因此我们只需取某个坐标最大的点作为第一个点（它一定是凸包顶点），再将所有点拍平到某个面上（例如坐标面）按二维卷包裹法求一条相邻的边即可。
/*将 z 拍平的向量叉积*/
double mul2(vect u, vect v){ return u.x*v.y-u.y*v.x; }
/*将 z 拍平*/
inline bool onright(vect v, vect f1, vect f2){ return gtr(mul2(v-f1, f2-f1), 0); }
{/*求初始边*/
int s[2];
s[0] =0, s[1] =-1;
for(int i =1; i < totp; ++i)
if(gtr(pts[i].x, pts[s[0]].x))/*取 x 坐标最大的点*/
s[0] =i;
for(int i =0; i < totp; ++i){
if(i == s[0]) continue;
if(s[1] == -1 || onright(pts[i], pts[s[0]], pts[s[1]])/*点是否在直线 (拍平 )右侧*/)
s[1] =i;
}
}我们同样也考虑多个点共面的特殊情况。
首先可以想到一种比较 "友好" 的共面情况：
可以想到我们只需取和当前边所做三角形面积最大的的点就可以了。
但接下来再考虑这种情况：
虽然我们对求得面要求并不严格，但这种重叠的面一定是不合法的。
或许可以想到将这些点映射到它们的平面上并求一次二维凸包，再一次性做它们面。由于每个点至多只会被带入做一次二维凸包，所以这种方法的复杂度是可以接受的。
但其代码复杂度难以想象...
相对而言更有效的方法是直接对点扰动，这样就可以直接不考虑共面的问题。
（如果有比映射求二维凸包更简洁的方法请务必联系我 QAQ）
（这里的实例是 求给定点集的凸包的表面积）
（仅含有扰动实现的代码）
下面的代码进行了多次扰动。
主要原因是增量法只有在一开始有共面的危险（只用一个面做初始数据），但卷包裹法在每一次循环中都是要求不能超过三个点共面的。
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <algorithm>
using std::pair;
typedef pair<int, int> pad;
const int MAXN =2500;
/*------------------------------Computational geometry------------------------------*/
const double eps =1e-9;
struct vect{
double x, y, z;
vect(){}
vect(double xx, double yy, double zz):x(xx), y(yy), z(zz){}
vect operator - (vect v){ return vect(x-v.x, y-v.y, z-v.z); }
vect operator / (vect v){ return vect(y*v.z-z*v.y, z*v.x-x*v.z, x*v.y-y*v.x); }
double operator * (vect v){ return x*v.x+y*v.y+z*v.z; }
double m(){ return sqrt(x*x+y*y+z*z); }
}pts[MAXN];
struct plane{
int v[3];
plane(){}
plane(int uu, int vv, int ww){ v[0] =uu, v[1] =vv, v[2] =ww; }
vect normal(){ return (pts[v[1]]-pts[v[0]])/(pts[v[2]]-pts[v[0]]); }
};
inline bool gtr(double a, double b){ return (a-b > eps); }
inline double dist(vect x, vect y){ return (x-y).m(); }
inline bool isabove(vect v, plane p){ return gtr((v-pts[p.v[0]])*p.normal(), 0); }
double mul2(vect u, vect v){ return u.x*v.y-u.y*v.x; }
/*将 z 拍平*/
inline bool onright(vect v, vect f1, vect f2){ return gtr(mul2(v-f1, f2-f1), 0); }
/*------------------------------Convex Hulls------------------------------*/
bool vise[MAXN][MAXN];
pad q[MAXN<<1]; int head, tail;
inline int getConvexHulls(int totp, plane facets[]){
int s[2];
s[0] =0, s[1] =-1;
for(int i =1; i < totp; ++i)
if(gtr(pts[i].x, pts[s[0]].x))
s[0] =i;
for(int i =0; i < totp; ++i){
if(i == s[0]) continue;
if(s[1] == -1 || onright(pts[i], pts[s[0]], pts[s[1]]))
s[1] =i;
}
int totf =0;
q[tail++] =pad(s[0], s[1]);
q[tail++] =pad(s[1], s[0]);
while(tail-head){
pad ne =q[head++];
if(vise[ne.first][ne.second]) continue;
vise[ne.first][ne.second] =1;
int p =-1;
plane nf;
for(int i =0; i < totp; ++i){
/*避免当前边顶点作为 p 的初始值 (p == -1)，不然会导致面的三点共线*/
if(i == ne.first || i == ne.second) continue;
if(p == -1 || isabove(pts[i], nf))
p =i, nf =plane(ne.first, ne.second, p);
}
facets[totf++] =nf;
vise[ne.second][p] =vise[p][ne.first] =1;
q[tail++] =pad(p, ne.second); q[tail++] =pad(ne.first, p);
}
return totf;
}
/*------------------------------Main------------------------------*/
plane facets[MAXN<<1];
double reps2() { return (rand()%2) ? eps : -eps; }
double reps() {return (rand()/(double)RAND_MAX-0.5)*eps+reps2();}
int main(){
int n; scanf("%d", &n);
for(int i =0; i < n; ++i) scanf("%lf%lf%lf", &pts[i].x, &pts[i].y, &pts[i].z);
srand(time(0));
/*这个程序可以拍过约一百个随机数据*/
/*具体扰动次数和大小可以因情况而定*/
for(int k =0; k < 5; ++k) for(int i =0; i < n; ++i) pts[i].x +=reps(), pts[i].y +=reps(), pts[i].z +=reps();
int h =getConvexHulls(n, facets);
double area =0;
for(int i =0; i < h; ++i) area +=facets[i].normal().m()/2;
printf("%.3lf", area);
}是
不过不知道为什么我写实现常数非常大，实际效率和我的增量法一样...（可能是因为过多的结构体构造）
这种做法是不支持在线的。
另外，个人实测下来卷包裹法的出错率是比增量法要高不少的（不如说后者几乎不会出错...），如果真的需要在两者间选择还是建议写增量法。
这种做法因为其和 快速排序 类似的启发式划分而得名。
（如果有人知道二维的 quickhull 的话，注意其实三维的 quickhull 和二维的区别还是蛮大的...）
对于一个由最终凸包顶点组成的面和另一个点组成的棱锥，所有在其内部的点都是可以直接排除的。
具体来说，我们有一个均由最终凸包顶点组成的多面体
我们取离这个面最远的点，和多面体对于这个点的地平线连面更新多面体，可以获得多个新面
将
对于不满足条件的点，就直接抛弃；
然后我们再将
可以参考这个流程图（不太好翻译，直接用原文可能更清晰）：
这里还有一个 quickhull 算法的 demo。
截图：
（可见面）
（地平线）
（新建面）
完整 demo 在这（baidu网盘，提取码：0eot）
为了找到一个由最终凸包的顶点组成的单纯形作为算法的开始，我们可以执行以下步骤：
（来自，不过注意第三步的不同，最坏情况下这
这
然后我们再遍历所有点，给这个单纯形的四个面分配点即可。
（这里的实例是 求给定点集的凸包的表面积）
代码调试非常麻烦。
实例代码的结构很难看...还请见谅（暂时没啥时间改）；仅供参考思路。
（另外实现中删除面和流程图没严格对应，因为在查询地平线已经标记了删除的面，维护邻接信息时也抛弃了删除的面）
//输入退化为面或点未测试
#include <cstdio>
#include <cmath>
#include <vector>
#include <list>
#include <queue>
using std::vector;
using std::list;
using std::queue;
using std::pair;
using std::max;
typedef pair<int, int> pad;
const int MAXN =2500;
/*------------------------------Computational geometry------------------------------*/
const double eps =1e-8;
struct vect{
double x, y, z;
int id;
vect(){}
vect(double xx, double yy, double zz):x(xx), y(yy), z(zz){}
vect operator - (vect v){ return vect(x-v.x, y-v.y, z-v.z); }
vect operator / (vect v){ return vect(y*v.z-z*v.y, z*v.x-x*v.z, x*v.y-y*v.x); }
double operator * (vect v){ return x*v.x+y*v.y+z*v.z; }
double m(){ return sqrt(x*x+y*y+z*z); }
};
struct line{
vect u, v;
line(){}
line(vect uu, vect vv):u(uu), v(vv){}
};
struct plane{
vect vec[3];
plane(){}
plane(vect uu, vect vv, vect ww){ vec[0] =uu, vec[1] =vv, vec[2] =ww; }
vect normal(){ return (vec[1]-vec[0])/(vec[2]-vec[0]);}
vect u(){ return vec[0]; }
};
inline bool gtr(double a, double b){ return a-b > eps; }
inline bool eq(double a, double b){ return (a-b > -eps && a-b < eps); }
inline bool eq(vect u, vect v){ return (eq(u.x, v.x) && eq(u.y, v.y) && eq(u.z, v.z)); }
inline double Abs(double x){ return gtr(0, x) ? -x : x; }
/*带符号距离*/
inline double dist(vect v, plane p){ return (v-p.u())*p.normal()/p.normal().m(); }
/*不带符号*/
inline double dist(vect v, line f){ return ((f.v-f.u)/(v-f.u)).m()/(f.v-f.u).m(); }
inline double dist(vect u, vect v){ return (u-v).m(); }
inline bool isabove(vect v, plane p){ return gtr((v-p.u())*p.normal(), 0); }
/*------------------------------Convex Hulls------------------------------*/
int TIME =0;/*全局时间戳*/
struct facet{
int n[3];/*neighbor，和点对应 (u->v, v->w, w->u)*/
int id, vistime/*访问的时间戳*/;
bool isdel;
plane p;
facet(){ vistime =isdel =0; }
facet(plane pp):p(pp){ vistime =isdel =0; }
facet(int idd, plane pp):id(idd), p(pp){ vistime =isdel =0; }
void in(int n1, int n2, int n3){ n[0] =n1, n[1] =n2, n[2] =n3; }
};
/*地平线的边*/
struct edge{
int netid, facetid;
};
/*存储所有面*/
vector<facet> FAC;
struct ConvexHulls3d{
int index/*索引面*/;
double surfacearea;
ConvexHulls3d(int indd):index(indd){ surfacearea =0; }
void dfsArea(int nf){
if(FAC[nf].vistime == TIME) return;
FAC[nf].vistime =TIME;
surfacearea +=FAC[nf].p.normal().m()/2;
for(int i =0; i < 3; ++i)
dfsArea(FAC[nf].n[i]);
}
double getSurfaceArea(){
if(gtr(surfacearea, 0)) return surfacearea;
++TIME;
dfsArea(index);
return surfacearea;
}
int getHorizon(int f, vect &p, int vistime[], edge e1[], edge e2[], vector<int> &resfdel){
if(!isabove(p, FAC[f].p)) return 0;
if(FAC[f].vistime == TIME) return -1;
FAC[f].vistime =TIME;
FAC[f].isdel =1;/*顺便标记删除的面*/
resfdel.push_back(FAC[f].id);
int ret =-2;
for(int i =0; i < 3; ++i){
int res =getHorizon(FAC[f].n[i], p, vistime, e1, e2, resfdel);
if(res == 0){
int pt[2];
pt[0] =FAC[f].p.vec[i].id, pt[1] =FAC[f].p.vec[(i+1)%3].id;
for(int j =0; j < 2; ++j){
if(vistime[pt[j]] != TIME){
vistime[pt[j]] =TIME;
e1[pt[j]].netid =pt[(j+1)%2];
e1[pt[j]].facetid =FAC[f].n[i];
}
else{
e2[pt[j]].netid =pt[(j+1)%2];
e2[pt[j]].facetid =FAC[f].n[i];
}
}
ret =pt[0];
}
else if(res != -1 && res != -2/*被围在中间的面*/) ret =res;
}
return ret;
}
};
/*----------------------------------------------------------------*/
/*全局点*/
vector<vector<vect> > pts;
/*构造初始单纯形*/
inline ConvexHulls3d getStart(vect point[], int totp){
vect pt[6], s[4];
for(int i =0; i < 6; ++i) pt[i] =point[1];
/*取坐标轴最大点*/
for(int i =2; i <= totp; ++i){
if(gtr(point[i].x, pt[0].x)) pt[0] =point[i];
if(gtr(pt[1].x, point[i].x)) pt[1] =point[i];
if(gtr(point[i].y, pt[2].y)) pt[2] =point[i];
if(gtr(pt[3].y, point[i].y)) pt[3] =point[i];
if(gtr(point[i].z, pt[4].z)) pt[4] =point[i];
if(gtr(pt[5].z, point[i].z)) pt[5] =point[i];
}
s[0] =pt[0], s[1] =pt[0], s[2] =pt[0], s[3] =pt[0];
/*取距离最大的两个点*/
for(int i =0; i < 6; ++i) for(int j =i+1; j < 6; ++j)
if(gtr(dist(pt[i], pt[j]), dist(s[0], s[1])))
s[0] =pt[i], s[1] =pt[j];
/*取距离上两个点所连直线距离最远的点*/
for(int i =0; i < 6; ++i)
if(gtr(dist(pt[i], line(s[0], s[1])), dist(s[2], line(s[0], s[1]))))
s[2] =pt[i];
/*取所有点集中距离该面最远的点*/
for(int i =1; i <= totp; ++i)/*!!*/
if(gtr(Abs(dist(point[i], plane(s[0], s[1], s[2]))), Abs(dist(s[3], plane(s[0], s[1], s[2])))))
s[3] =point[i];
/*确保接下来构造的面是朝单纯形外的*/
if(gtr(0, dist(s[3], plane(s[0], s[1], s[2])))){
vect tmp =s[1]; s[1] =s[2]; s[2] =tmp;
}
/*构造单纯形*/
int f[4];
for(int i =0; i < 4; ++i) FAC.push_back(facet()), f[i] =FAC.size()-1, FAC[f[i]].id =f[i];
FAC[f[0]].p =plane(s[0], s[2], s[1]),/*底面*/
FAC[f[1]].p =plane(s[0], s[1], s[3]),
FAC[f[2]].p =plane(s[1], s[2], s[3]),
FAC[f[3]].p =plane(s[2], s[0], s[3]);
FAC[f[0]].in(f[3], f[2], f[1]);
FAC[f[1]].in(f[0], f[2], f[3]);
FAC[f[2]].in(f[0], f[3], f[1]);
FAC[f[3]].in(f[0], f[1], f[2]);
/*给四个面分配点集空间*/
for(int i =0; i < 4; ++i)
pts.push_back(vector<vect>());
/*给四个面分配点*/
for(int i =1; i <= totp; ++i){
if(eq(point[i], s[0]) || eq(point[i], s[1]) || eq(point[i], s[2]) || eq(point[i], s[3])) continue;
for(int j =0; j < 4; ++j)
if(isabove(point[i], FAC[f[j]].p)){
pts[f[j]].push_back(point[i]);
break;
}
}
/*返回初始单纯形，以一个面作为索引*/
return ConvexHulls3d(f[0]);
}
edge e[2][MAXN] /*边界线的图信息*/;
int vistime[MAXN] /*每个点访问的时间戳*/;
queue<int> que;
vector<int> resfnew /*保存新构造的面*/, resfdel /*保存删除的面*/;
vector<vect> respt /*保存需要分配的点*/;
inline ConvexHulls3d quickHull3d(vect point[], int totp){
ConvexHulls3d hull =getStart(point, totp);
/*将初始单纯形的面加入队列*/
que.push(hull.index);
for(int i =0; i < 3; ++i)
que.push(FAC[hull.index].n[i]);
/*snew 保存最后返回的凸包的索引面*/
int snew =0;
while(que.size()){
int nf =que.front(); que.pop();
/*当前面已被删除，跳过*/
if(FAC[nf].isdel) continue;
/*当前面没有分配到顶点，跳过*/
if(pts[nf].size() == 0){
snew =nf;/*确保面最后存在*/
continue;
}
/*取距离该面最远的点*/
vect p =pts[nf][0];
for(int i =1; i < (int)pts[nf].size(); ++i)
if(gtr(dist(pts[nf][i], FAC[nf].p), dist(p, FAC[nf].p)))
p =pts[nf][i];
/*求地平线，可以知道得到的地平线至少有三个点*/
++TIME;
resfdel.clear();
/*当前面一定会被删除，因此直接从当前面 dfs*/
int s =hull.getHorizon(nf, p, vistime, e[0], e[1], resfdel);
/*遍历地平线（绕一圈），构造新面*/
/*在求地平线时我们无法得知某条边是逆时针还是顺时针的，因此需要这里判断*/
resfnew.clear();
++TIME;
int from =0 /*上一个访问的点*/, lastf =0 /*上一个新建的面*/, fstf =0 /*第一个新建的面*/;
while(vistime[s] != TIME){
/*用时间戳记录当前点是否访问*/
vistime[s] =TIME;
int net/*下一个点*/;
int f /*地平线上当前边所接的无法看见的面*/, fnew /*新面*/;
/*确保遍历方向正确*/
if(e[0][s].netid == from) net =e[1][s].netid, f =e[1][s].facetid;
else net =e[0][s].netid, f =e[0][s].facetid;
/*求出这两个点在邻接面上的逆顺时针信息*/
int pt1 =-1, pt2 =-1;
for(int i =0; i < 3; ++i){
if(eq(point[s], FAC[f].p.vec[i])) pt1 =i;
if(eq(point[net], FAC[f].p.vec[i])) pt2 =i;
}
/*确保 pt1->pt2 是按邻接面的点的逆时针排列*/
if((pt1+1)%3 != pt2) pt1 ^=pt2 ^=pt1 ^=pt2;/*交换*/
/*这样构造的面是朝凸包外的*/
FAC.push_back(facet(plane(FAC[f].p.vec[pt2], FAC[f].p.vec[pt1], p)));
fnew =FAC.size()-1, FAC[fnew].id =fnew;
pts.push_back(vector<vect>());
resfnew.push_back(fnew);
/*维护邻接信息*/
FAC[fnew].n[0] =f, FAC[f].n[pt1] =fnew;
if(lastf){
/*不能预先确定是顺时针遍历还是逆时针遍历*/
/*维护新建的面之间的邻接信息*/
if(eq(FAC[fnew].p.vec[1], FAC[lastf].p.vec[0])) FAC[fnew].n[1] =lastf, FAC[lastf].n[2] =fnew;
else FAC[fnew].n[2] =lastf, FAC[lastf].n[1] =fnew;
}
else fstf =fnew;/*还未新建面*/
lastf =fnew;
from =s;
s =net;
}
/*给新建的面头尾维护临界信息*/
if(eq(FAC[fstf].p.vec[1], FAC[lastf].p.vec[0])) FAC[fstf].n[1] =lastf, FAC[lastf].n[2] =fstf;
else FAC[fstf].n[2] =lastf, FAC[lastf].n[1] =fstf;
/*取得所有需要分配的点*/
respt.clear();
for(int i =0; i < (int)resfdel.size(); ++i){
for(int j =0; j < (int)pts[resfdel[i]].size(); ++j)
respt.push_back(pts[resfdel[i]][j]);
pts[resfdel[i]].clear();
}
/*分配点*/
for(int i =0; i < (int)respt.size(); ++i){
if(eq(respt[i], p)) continue;/*跳过用于新建面的点*/
for(int j =0; j < (int)resfnew.size(); ++j)
if(isabove(respt[i], FAC[resfnew[j]].p)){
pts[resfnew[j]].push_back(respt[i]);
break;/*确保点不会被重复分配*/
}
}
/*将新的面加入队列*/
for(int i =0; i < (int)resfnew.size(); ++i)
que.push(resfnew[i]);
}
hull.index =snew;
return hull;
}
inline void preConvexHulls(){
/*0 位做保留*/
pts.push_back(vector<vect>());
FAC.push_back(facet());
}
/*------------------------------Main------------------------------*/
vect point[MAXN];
int main(){
preConvexHulls();
int n; scanf("%d", &n);
for(int i =1; i <= n; ++i){
double x, y, z; scanf("%lf%lf%lf", &x, &y, &z);
point[i] =vect(x, y, z);
point[i].id =i;
}
ConvexHulls3d hull =quickHull3d(point, n);
printf("%.3lf", hull.getSurfaceArea());
}时间复杂度据说是期望
空间复杂度是
按我上面的实现,；我自己测试了下，开 O2，
我们建立一个图来表示点与已有面的可见关系：若一个点能看见一个面（该点在面上方），我们在它们代表的节点间连一条无向边。
可以发现它是一个二部图，左边是所有点，右边是已有的面。其中每条边都表示了一个点与一个面间的冲突关系（加入该点一定会删除这个面），因此我们把它叫做 "冲突图"。
考虑如何在加入新面时维护冲突图。
一个简单的思路是枚举所有点并连边。这样做每次维护是
在我们新建面时，它一定是由已有凸包上的一边（地平线上）和当前点形成的面。且我们可以证明
因此我们只需枚举能看到
[1] 实际上新面的 "角度" 范围是这样的：
我们初始先随意构建一个单纯形
接着随机打乱（这也是算法名称的由来）并遍历剩下的点。
对于每次循环：
根据冲突图找到所有当前点能看到的面，这些面一定是相邻的，因此可以遍历它们求出地平线。
然后再根据地平线新建面。对于每个新面，如果它和地平线当前边邻接的（不可见）面处于同一平面，就在冲突图中 "合并" 这两个面；否则就往冲突图中加入新面，并维护冲突图。
最后将当前点标记为已使用。
另外还可以参考这张流程图：（感觉原文更清晰就不翻译了 qwq）
[1] 不过这里构建单纯形越大越能在开始排除更多点。
有证明指出随机点序的期望时间复杂度是
空间复杂度应该和时间复杂度同阶（冲突图的边信息）。
[1] 貌似带了一个
这个比较简单，直接遍历面求出即可。
和二维一样，我们同样可以用有向体积割补求出凸包的面积。
具体的来说，定义凸包面和其下方的点形成的棱锥体积为正，我们从某个点（例如原点）向每个面做有向体积，其和就是凸包的体积。
（凸包）
（贡献为负的体积）
（贡献为正的体积）
（交互视图的链接 在这）
即每个面对体积的贡献如下：
volume +=-(dist(vect(0, 0, 0), facet.p)*(facet.p.normal().m()/2))/3/*V_棱锥 =h*s/3*/;
这个方法对任意多面体求体积也适用。
由于我们存储的面可能仅仅是凸包完整面的一部分，因此不可以直接输出存储面的个数。
首先可以想到直接暴力判断：
int cnt =0;
for(int i =0; i < totf; ++i){
++cnt;
for(int j =i+1; j < totf; ++j)
if(eq(facets[i].p, facets[j].p)){
/*数组右侧有相同面就去掉贡献*/
/*这样同样的一组面只有最右边的有贡献*/
--cnt;
break;
}
}如果只是单纯的将面存在数组里，这种方式已经是基本最优的了（实在想还可以优化下，不过没什么意义）。
但如果存储了面的邻接信息，实际上我们有一种
可以想到凸包上相等的面都是连续的，因此可以考虑一次深搜仅遍历同平面上的凸包面，并将它们邻接的不同平面上的面加入队列，再搜索队列下一个面。
要注意队列中的面也有可能是已被访问过的。（可以考虑构造一种搜索顺序）
bool eq(plane p, plane q){
vect a =p.v[0], b =p.v[1],c =p.v[2];
return (eq(dist(a, q), 0) && eq(dist(b, q), 0) && eq(dist(c, q), 0));
}
bool dfsCnt(int nf, int pref, queue<int> &Q){
if(FAC[nf].vistime == TIME) return 0;/*已访问，扩张失败*/
if(pref != 0 && !eq(FAC[nf].p, FAC[pref].p)){
Q.push(nf);
return 0;/*不在同一平面，扩张失败*/
}
FAC[nf].vistime =TIME;
for(int i =0; i < 3; ++i)
dfsCnt(FAC[nf].n[i], nf, Q);
return 1;
}
int getFacetCnt(){
if(cntface) return cntface;
++TIME;
resQ.push(index);
while(resQ.size()){
int nf =resQ.front(); resQ.pop();
if(dfsCnt(nf, 0, resQ)) ++cntface;
}
return cntface;
}我们可以考虑记录每个点是否访问，直接遍历存储的凸包面统计每个小面的三个顶点：
bool vistime[MAXNVECT];
void dfsCnt2(int nf){
if(FAC[nf].vistime == TIME) return;
for(int i =0; i < 3; ++i)
if(vistime[FAC[nf].p.v[i].id] != TIME){
vistime[FAC[nf].p.v[i].id] =TIME;
++cntVect;
}
for(int i =0; i < 3; ++i)
dfsCnt(FAC[nf].n[i], nf, Q);
}
int getVectCnt(){
if(cntVect) return cntVect;
++TIME;
TST =0;
dfsCnt2(index);
return cntVect;
}首先说一下三棱锥的重心：为其四个顶点的向量和除
/*四个顶点分别为 A, B, C, D*/
vect g =(A+B+C+D)/4;而求凸包重心具体方法是：先定一个点（例如坐标原点），我们从该点向存储的每个小面连一个三棱锥，然后再按有向体积比例带权求和所有重心的坐标。
具体的实现可以看代码：
/*为了使代码思路清晰，这里假定已经将存储的凸包所有小面读入一个数组*/
facets[MAXN];/*存储所有面*/
int totf;/*面总数*/
vect getG(){
vect G(0, 0, 0);
double sumV =0;
for(int i =0; i < totf; ++i){
vect ng =(vect(0, 0, 0)+facets[i].v[0]+facets[i].v[1]+facets[i].v[2])/4;
double nv =-(dist(vect(0, 0, 0), facets[i])*(facets[i].normal().m()/2))/3;
sumV +=nv;
G =G+ng*nv;
}
return G/sumV;
}（来自）
（其证明我还没找到 qaq，暂时先咕着）
一开始真的没想到这东西这么复杂...
一时起意就硬着头皮啃下来了，然后觉得既然都学了就顺便写下来当做思路整理了（正好网上也没多少关于三维计算几何的博客...）
这次写这篇文章没什么感想...可能之后还会回来补充后记）qwq。
如果有什么问题可以私信告诉我。