计算几何三维入门
Piwry
2020-07-04 20:59:41
本文主要介绍了三维计算几何的基本运算，并取其比较基本的求凸包问题介绍了几种算法。
另外我还整理了我找到的一些资料，可以 [从这](https://pan.baidu.com/s/1wbdMtgas3Ggu4NeJDCj8ww
)（baidu 网盘，6ahh ）下载。
# 目录
0. 目录
1. 三维计算几何基础
2. 三维凸包
3. 三维凸包上的一些运算
4. 后记
---
# 1. 三维计算几何基础
## 1.1. 空间向量
### 向量缩放
同二维
```cpp
vect operator * (double k){ return vect(x*k, y*k, z*k); }
vect operator / (double k){ return vect(x/k, y/k, z/k); }
```
### 向量加减
同二维
意义也一样
```cpp
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); }
```
### 向量点积（数值积）
同二维
意义也和二维一样，是一个向量到另一个向量的投影长乘另一个向量的模长
满足交换律
```cpp
double operator * (vect v){ return x*v.x+y*v.y+z*v.z; }
```
### 向量叉积
```cpp
vect operator / (vect v){ return vect(y*v.z-z*v.y, z*v.x-x*v.z, x*v.y-y*v.x); }
```
返回的是一个空间向量
它的**默认**方向和两个输入向量满足右手定理（形如三维坐标轴：$x$ 逆时针到 $y, z$，$y$ 逆时针到 $z$。返回的向量可以看做 $z$ 轴），
长度是输入向量构成的平行四边形面积，
**正负**由第一个向量到第二个向量的逆时针夹角决定（同二维），
不满足交换律
3b1b 中（[见此](https://www.bilibili.com/video/BV1vs411t7w9)）有对其含义的一种解释：
将输入向量的起点置于原点，若定义运算 $*$ 表示求一个点和输入向量构成的平行四边形构成的平行六面体的有向体积 $^1$，
那么做这个运算 $*$ 就相当于和输入向量的**叉积**（叉积方向为第一个叉第二个）做点积。
---
[1] 有向体积：设输入向量 $u, v$，原点 $O$，则有平面（具体存储规则见下文） $(O, O+u, O+v)$。若点在平面上方则体积为正，否则为负。
### 向量模长
```cpp
double m(){ return sqrt(x*x+y*y+z*z); }
```
投影到坐标面化简下式子就可以。
## 1.2. 直线
### 直线的存储
```cpp
struct line{
vect v[2];
line(){}
line(vect uu, vect vv){ v[0] =uu, v[1] =vv; }
};
```
和二维一样存储直线上的两个点。
该直线的方向定义为从存储的第一个点到第二个点的向量。
### 点与直线的距离
```cpp
inline double dist(vect v, line f){ return ((f.v-f.u)/(v-f.u)).m()/(f.v-f.u).m(); }
```
做叉积得到平行四边形的面积，除掉底长（存储直线也是存点）就是高，这个高即是距离。
## 1.3. 平面
### 面的存储
```cpp
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()`）即为从一个点顺次连向后两个点的向量的叉积（顺次）。这个向量竖直指向面的**上方**。
### 判断点与面关系
```cpp
inline bool isabove(vect v, plane p){ return gtr((v-p.u())*p.normal()/*这是投影长度乘法向量模长，只是因为不需要知道数值只需要符号，所以没做处理*/, 0); }
```
只需判断法向量起点到该点的这个向量在法向量上的投影符号即可。
### 点与面的距离
```cpp
/*带方向距离*/
inline double dist(vect v, plane p){ return (v-p.u())*p.normal()/p.normal().m(); }
```
法向量起点到该点的这个向量在法向量上的投影长。
### 判断面相等
一个方法是判断两个面的法向量是否相等。而判断法向量相等，我们可以将它分别投影到 $x-y, x-z$ 平面求极角，并比较极角（有点像球极坐标系的思路）
```cpp
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)));
}
```
这种方法要求两个相等的面**朝向**也是相等的。
不过考虑到极角为 $-\pi$ 和极角为 $\pi$ 的二维向量的夹角实际上很小，加上求法向量本就有不小的精度误差，这种方法很有可能在实际运用时**误判**。
还可以想到枚举某个面存储的三个代表点，看它们是不是**都**在另一个面上
```cpp
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));
}
```
但这种方法不会判断平面的朝向。如果实在需要，可以做通过平面外一点与两个平面的有向体积来判断
```cpp
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;
}
```
## 1.4. 面积、体积相关
通过叉积的定义，我们可以很轻松地求出两个向量形成的平行四边形的面积：
```cpp
vect u, v;
double S =(u/v).m();
```
体积问题则可以转化成底面积乘高来计算。高可以用点积的意义求出（其实也就是点到面的距离）。
举个最简单的例子，给定顶点，计算三棱锥的体积。
首先给出三棱锥体积的计算公式：$\frac 1 3 ah$（$a$ 为底面积）
那么其代码就为：
```cpp
vect A, B, C, D;/*四个顶点，具体位置见下图*/
vect n =(B-A)/(C-A);
double V =(n.m()/2)*(D*n/n.m())/3;
```
同时，许多复杂多面体体积其实也都可以转化为多个三棱锥体积的和。
## 1.5. 多面体
### 多面体的储存
通常我们通过储存**面**来储存多面体。
对于多面体每个面的面积，我们可以用叉积的含义轻易算出：
```cpp
plane facet;
double S =facet.normal().m()/2;/*其面积*/
```
在某些情况下，我们还需要存储每个面的**邻接面信息**。
为了方便，我们将每个顶点**大于三个**的面划分成若干个三角形的面；视具体情况可能划分后的面的顶点**不属于原来面的顶点**。
于是我们就可以这样保存面：
```cpp
struct facet{
int n[3];/*neighbor，和点对应 (u->v, v->w, w->u)*/
plane p;
};
```
对于一个多面体，我们保存其中一个面作为**索引**，就可以通过邻接信息读取多面体的表面。
### 地平线
（貌似一般都叫 "Horizon"...直译过来就是地平线...）
对于一个点和一个多面体，从该点能 "看到" 的所有面一定是一片连续的面，而**这片面的边缘**就是这个点在这个多面体上的地平线。
（可想象站在高处望向远处的场景...）
若要求得地平线，我们可以从某个该点能看到的面（具体来说就是这个面的**平面**在该点下面）开始 dfs（这个面也可以从索引面开始 dfs 求出），并在图中标记所有**边缘面**；
对于每个边缘面，我们将每条**边缘边**的两个顶点**建图**连边。
可以知道我们最后求出的图一定是一个**环**。
同时为了避免不必要的计算，我们可以在 dfs 顺便返回一个环上的一个点作为索引。
```cpp
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;
}
```
### 多面体面与边的关系
我们可以想到，一条边会且仅会被两个面使用；按我们储存面的方式，这条边一定是被 "逆时针用一次"，再 "顺时针用一次"。
就像下图中那样：
这个性质可以帮助我们**通过统计边来获得多面体面的信息**。
## 1.6 点的随机扰动
### 直接扰动
为了避免可能的四点共面，三点共线等问题增加问题复杂度，我们可以给点的坐标加减一个极小的值，同时也能保证一些问题中求出的答案不会偏差太大。
扰动的具体实现和应用的实际情况有关，
如果要扰动距离较**近**的**少数**点，可以这样实现（其中 `eps` 是判断两个数不同的最小差值，之后（应该）不再解释）：
```cpp
double reps() { return (rand()%2) ? eps : -eps; }
```
如果要扰动距离较**远**的**多数**点，可以这样实现：
```cpp
double reps() {return (rand()/(double)RAND_MAX-0.5)*eps;}
```
（较远的点对会放大精度误差，因此扰动值可以不大于 `eps`；点数目的增加要求扰动的值更加随机）
如果不确定问题具体是哪种情况，其实也可以综合两种实现：
```cpp
double reps2() { return (rand()%2) ? eps : -eps; }
double reps() {return (rand()/(double)RAND_MAX-0.5)*eps+reps2();}
```
### 自适应扰动
我们还是给每个点一个随机的误差，但是我们将每个点的坐标显式地表示成 $p+q\cdot eps$。
在运算时，首先将 $p$ 带入计算，如果相同再考虑 $q\cdot eps$。
这种方法对精度的影响应该会比较小；如果统计贡献时并不要求排除特殊情况（比如四点共面什么），这种方法甚至能做到无精度影响。
### 对结果的影响
点扰动对一些不要求求出**准确值**的问题是影响不大的，例如：
1. 凸包表面积
2. 凸包体积
3. 凸包重心
但有一些问题是存在准确解的，这时候我们可能要改进扰动算法（例如保存原坐标）或选用不需要扰动点的算法，例如：
1. 求凸包表面数
2. 求在凸包表面的点数
---
# 2. 三维凸包
这里先介绍下较为广义的凸包定义：
对于 $N$ 维的凸包，它是包含所有给定点的最小凸集，更确切的说，是所有可能的由 $N+1$ 个由给定点构成的 $N-1$ 维凸包构成的 $N$ 维形的交集。递归定义的终点是单个点。
而对于一般的三维凸包问题，就是所有可能的由给定点集构成的**面**构成的单纯形的交集。
下面将介绍几种代表性的求三维凸包的方法：增量法，卷包裹法（gift wrap），QuickHull，随机增量法。
其中前两个是 $O(n^2)$ 级别的，后两个是 $O(nlogn)$ 级别的（QuickHull 的复杂度仍存疑）；只有增量法支持在线。
## 2.1. 增量法
### 基本思路
我们考虑加入一个点时如何维护凸包
可以发现只有它能 "看到" 的面，或者说**在该点下方**的面（对于点**恰好**在面上，我们**不算**能 "看到"）会被删除。
所以初始先构造一个**任意**单纯形（四个顶点的多面体）。我们每加入一个点，就遍历所有面找出能看到的面，将它们删除；并对 "缺口" 的边缘边和新点连面，加入凸包。
### 求缺口边缘（地平线）
对于边缘的标记，我们可以想到：一条边会且仅会被两个面使用；按我们储存面的方式，这条边一定是被 "逆时针用一次"，再 "顺时针用一次"。（这里 copy 自上面）
于是用一个二维数组标记，只对**只使用一次**的边，按**使用过的方向**和新点连面即可（可见代码）。
### 对特殊情况的考虑
有时可能会出现几个点共面的情况。
我们考虑增量法的实现，可以发现在已有面上的点是无效的，也不会产生贡献。
不过由于我们连的面都是三角形，因此最后求得的面集合可能会把原来一个多边形面划分成**多个**三角形小面（同时也不保证三角形的顶点一定在**边**上）。注意后续运算需要考虑**去重**的问题。
另外，值得一提的是，如果我们通过扰动点来期望最多使三个点共面，我们初始可以**只构造面**。这样可以略微缩减代码量。
不扰动的话这样可能出错。考虑这样的数据：
由于扰动的值不需要太大，因此一般的精度要求用这种方法都可以满足。
### CODE
（这里的实例是 [求给定点集的凸包的表面积](https://www.luogu.com.cn/problem/P4724)）
扰动点（只构造面）：
```cpp
#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);
}
```
不扰动点（仅给出不同部分）：
```cpp
/*------------------------------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` 数组：
```cpp
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）：
随机点集下期望**每一步**的凸包是 “均匀" 的；也就是说，它的面是在各个方向是 ”平均“ 的，
随着凸包的成型，新加入的点可能 “看到” 的面越来越少；而即使 “看到” 的面最多，也只有半个凸包，也就是期望一半的面数。
也就是说，随机点集下，遍历看到的面的次数的**最大值**和遍历看不到的面的次数的**最小值**相等。
实际实现时按心情写就行了）
### 复杂度及其它
可以想到是 $O(n^2)$ 的。由于内部没有太多分支，常数应该不会太大。
这种做法可以很轻松地支持在线
## 2.2. 卷包裹法（gift wrap）
### 基本思路
我们想象用一张纸慢慢地包住凸包。
具体来说，我们想象有一个平面从凸包的**一个面的一条边**开始向**内**折，并取**第一个**碰到的点作为新的面的顶点，构造新面；这是一次打包（wrap）。
我们再将新面新产生的两条边**包含其面的信息**（具体来说是边的方向，参考上文介绍的 面与边的关系）放入队列，重复同样的工作，直到没有边可以再打包即可。
### 如何找到合适的点
比较直观的思路是求出所有点相对当前面的 "夹角"，不过这种方式难以实现。
我们再转化其含义，可以发现这实际上就是求一个新多边形面，使剩下的所有点都在新面**下方**（新面上方朝多面体外侧）。
```cpp
/*单次打包*/
/*大部分变量名实例同上文*/
/*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;
}
```
### 构建初始边
卷包裹法一开始需要一个属于最终凸包的边作为算法起点。
我们可以想到，如果将凸包 "**拍平**" 到一个平面上，其映射在平面上的点集的凸包的边一定也是原凸包的边。
因此我们只需取某个坐标最大的点作为第一个点（它一定是凸包顶点），再将所有点拍平到某个面上（例如坐标面）按二维卷包裹法求一条相邻的边即可。
```cpp
/*将 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）
### CODE
（这里的实例是 [求给定点集的凸包的表面积](https://www.luogu.com.cn/problem/P4724)）
（仅含有扰动实现的代码）
下面的代码进行了多次扰动。
主要原因是增量法**只有在一开始**有共面的危险（只用一个面做初始数据），但卷包裹法在**每一次循环**中都是要求不能超过三个点共面的。
```cpp
#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);
}
```
### 复杂度及其它
是 $O(nh)$ 的，其中 $h$ 是凸包点数；随机数据下减少的循环次数非常可观（自测了几个数据，大概是和 $\sqrt n$ 同级的，仅供参考）
不过不知道为什么我写实现常数非常大，实际效率和我的增量法一样...（可能是因为过多的结构体构造）
这种做法是不支持在线的。
另外，个人实测下来卷包裹法的出错率是比增量法要高不少的（不如说后者几乎不会出错...），如果真的需要在两者间选择还是**建议写增量法**。
## 2.3. QuickHull
这种做法因为其和 快速排序 类似的启发式划分而得名。
（如果有人知道二维的 quickhull 的话，注意其实三维的 quickhull 和二维的区别还是蛮大的...）
### 基本思路
对于一个由**最终凸包顶点**组成的面和另一个点组成的棱锥，所有在**其内部**的点都是可以直接排除的。
具体来说，我们有一个均由最终凸包顶点组成的多面体 $H'$，其一个面（或面的划分，保证是一个三角形） $\pi$，和它被分配到的一个点集 $P_\pi$。
我们取离这个面**最远**的点，和**多面体对于这个点的地平线**连面**更新多面体**，可以获得多个**新面** $F$，可以知道这些面的顶点一定是最终凸包的顶点；
将 $P_\pi$ 剩下的点**加上**被删除的面的点集 $P$，划分给 $F$，并保证每个点在划分到的面的**上方**（注意一个点可能满足划分给多个面，而二维 quickhull 就不会有这种情况）；
对于不满足条件的点，就直接抛弃；
然后我们再将 $F$ 加入队列，继续处理下一个面（注意有时会处理到一个已经删除的面）。
可以参考这个流程图（不太好翻译，直接用原文可能更清晰）：
这里还有一个 quickhull 算法的 demo。
截图：
（可见面）
（地平线）
（新建面）
完整 demo [在这](https://pan.baidu.com/s/194jOL5rhHxeE3eDQxhD7vg
)（baidu网盘，提取码：0eot）
### 初始化
为了找到一个由最终凸包的顶点组成的单纯形作为算法的开始，我们可以执行以下步骤：
1. 我们首先分别在 $x, y, z$ 三个坐标轴方向上找到最大最小点，共 $3$ 对 $6$ 个点（这 $6$ 个点必在凸包上）；
2. 在这 $6$ 个点当中，选取距离最大的两个点为初始点；
3. 遍历点集当中所有点，选取距离第二步所选 $2$ 个点所在直线距离最远的点作为第 $3$ 个点；
4. 遍历点集当中所有点，找到距离这 $3$ 个点所确定的平面最远的点作为第 $4$ 个点；
（[来自](https://zhuanlan.zhihu.com/p/33950536)，不过注意第三步的不同，最坏情况下这 $6$ 个点可能**只包含 $2$ 个不同的点**）
这 $4$ 个点即是最初的单纯形。
然后我们再遍历所有点，给这个单纯形的四个面分配点即可。
### CODE
（这里的实例是 [求给定点集的凸包的表面积](https://www.luogu.com.cn/problem/P4724)）
代码调试非常麻烦。
实例代码的**结构**很难看...还请见谅（暂时没啥时间改）；仅供参考思路。
（另外实现中删除面和流程图没严格对应，因为在查询地平线已经标记了删除的面，维护邻接信息时也抛弃了删除的面）
```cpp
//输入退化为面或点未测试
#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());
}
```
### 复杂度及其它
时间复杂度据说是期望 $O(nlogn)$ 的，最坏可能被卡到 $O(nh)$。常数可能较大。
空间复杂度是 $O(n)$ 的。
按我上面的实现,；我自己测试了下，开 O2，$10^6$ 的随机点大概可以跑到 $5s$ 左右（仅供参考）。
## 2.4. 随机增量法
### 冲突图
我们建立一个图来表示**点与已有面的可见关系**：若一个点能看见一个面（该点在面**上方**），我们在它们代表的节点间连一条无向边。
可以发现它是一个二部图，左边是所有点，右边是已有的面。其中每条边都表示了一个点与一个面间的冲突关系（加入该点一定会删除这个面），因此我们把它叫做 "冲突图"。
考虑如何在加入新面时维护冲突图。
一个简单的思路是枚举所有点并连边。这样做每次维护是 $O(n)$ 的。
在我们新建面时，它一定是由已有凸包上的一边（地平线上）和当前点形成的面。且我们可以证明 $^1$ 所有能看到新面（$f$）的点也是能看到这条边邻接的两个面（$f_1, f_2$）的**其中一个**的。
因此我们只需枚举能看到 $f_1, f_2$ 的点即可。
---
[1] 实际上新面的 "角度" 范围是这样的：
### 基本思路
我们初始先随意构建一个单纯形 $^1$，并根据这个单纯形的面建立初始冲突图。
接着**随机打乱**（这也是算法名称的由来）并遍历剩下的点。
对于每次循环：
根据冲突图找到所有当前点能看到的面，这些面一定是相邻的，因此可以遍历它们求出地平线。
然后再根据地平线新建面。对于每个新面，如果它和地平线当前边邻接的（不可见）面处于同一平面，就在冲突图中 "**合并**" 这两个面；否则就往冲突图中加入新面，并维护冲突图。
最后将当前点标记为已使用。
另外还可以参考这张流程图：（感觉原文更清晰就不翻译了 qwq）
---
[1] 不过这里构建单纯形越大越能在开始排除更多点。
### 复杂度及其它
有证明指出随机点序的期望时间复杂度是 $O(nlogn)$ 的 $^1$。由于点的顺序不由输入决定，因此这个期望复杂度应该是始终成立的。
空间复杂度应该和时间复杂度同阶（冲突图的边信息）。
---
[1] 貌似带了一个 $96$ 的大常数，具体证明见资料中的 Computational Geometry Algorithmsand Applications. Chapter 11.3
---
# 3. 三维凸包上的一些运算
## 3.1. 求凸包表面积
这个比较简单，直接遍历面求出即可。
## 3.2. 求凸包体积
和二维一样，我们同样可以用有向体积割补求出凸包的面积。
具体的来说，定义凸包面和其**下方**的点形成的棱锥体积为正，我们从某个点（例如原点）向每个面做有向体积，其和就是凸包的体积。
（凸包）
（贡献为负的体积）
（贡献为正的体积）
（交互视图的链接 [在这](https://www.geogebra.org/m/rzfnbcmf)）
即每个面对体积的贡献如下：
```cpp
volume +=-(dist(vect(0, 0, 0), facet.p)*(facet.p.normal().m()/2))/3/*V_棱锥 =h*s/3*/;
```
这个方法对任意多面体求体积也适用。
## 3.3. 求凸包表面个数
由于我们存储的面可能仅仅是凸包完整面的一部分，因此不可以直接输出存储面的个数。
首先可以想到直接暴力判断：
```cpp
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;
}
}
```
如果只是单纯的将面存在数组里，这种方式已经是基本最优的了（实在想还可以优化下，不过没什么意义）。
但如果存储了面的邻接信息，实际上我们有一种 $O(n)$ 的做法：
可以想到凸包上相等的面都是连续的，因此可以考虑一次深搜仅遍历同平面上的凸包面，并将它们邻接的不同平面上的面加入队列，再搜索队列下一个面。
要注意队列中的面**也有可能是已被访问过的**。（可以考虑构造一种搜索顺序）
```cpp
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;
}
```
## 3.4. 求凸包顶点个数
我们可以考虑记录每个点是否访问，直接遍历存储的凸包面统计每个小面的三个顶点：
```cpp
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;
}
```
## 3.5. 求凸包重心
首先说一下三棱锥的重心：为其四个顶点的向量和除 $4$
```cpp
/*四个顶点分别为 A, B, C, D*/
vect g =(A+B+C+D)/4;
```
而求凸包重心具体方法是：先定一个点（例如坐标原点），我们从该点向存储的每个小面连一个三棱锥，然后再按有向体积比例带权求和所有重心的坐标。
具体的实现可以看代码：
```cpp
/*为了使代码思路清晰，这里假定已经将存储的凸包所有小面读入一个数组*/
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;
}
```
（[来自](https://blog.csdn.net/qq_42555009/article/details/105279645)）
（其证明我还没找到 qaq，暂时先咕着）
---
# 后记
一开始真的没想到这东西这么复杂...
一时起意就硬着头皮啃下来了，然后觉得既然都学了就顺便写下来当做思路整理了（正好网上也没多少关于三维计算几何的博客...）
这次写这篇文章没什么感想...可能之后还会回来补充后记）qwq。
如果有什么问题可以私信告诉我。