最佳包裹

求三维凸包用增量法求解

假设我们已经维护好了前面的点的凸包，对于新加入的点：

如果这个点在凸包内部显然那就不用管了

如果这个点在凸包外部，那么考虑如下情况：

将新点\(P_r\)当做光源，照的到的面全部删掉，照不到的面保留下来即可

如何判断一个面是否能够被照到：取多面体外侧为正方向，则对于某一个面，如果光源在其正侧那么能被照到，否则不行；具体判断的时候，在平面上取一个点\(a\)，则判断由\(a,P_r\)组成的向量与平面的法向量的叉积即可

维护面的时候，我们都是维护的三角形（如果一个面不是三角形我们对其进行三角剖分），存储三角形的三个点就好了（按照逆时针方向存储）；一共有\(n\)个点，每个面对应三条边，但是每条边会被重复数两次，所以设面有\(m\)个，边有\(k\)个，则\(3m=2k\)，代入多面体欧拉定理得，\(m=2n-4,k=3n-6\)，可知时间复杂度为\(O(n^2)\)

如何找出分界线：设\(g[i][j]\)表示由\(i\)指向\(j\)的这条向量是否是分界线，那么由于我们是按照逆时针方向存储的三个点，如下图

我们只需要判断\(\vec{ij}\)和\(\vec{ji}\)表示的两个面中，一个可以被照到，另一个不能被照到就好了

代码见下：

#include<bits/stdc++.h>
#define ll long long
#define Vector Point
using namespace std;
const int N=2010;
const double eps=1e-10;
double Rand()
{return rand()/(double)RAND_MAX;//RAND_MAX是宏定义，表示rand()能返回的最大整数 
}
double reps()
{return (Rand()-0.5)*eps;
}
struct Point
{double x,y,z;Point(double a=0,double b=0,double c=0):x(a),y(b),z(c){}void shake(){x+=reps(),y+=reps(),z+=reps();}
}p[N];
Vector Cross(Vector a,Vector b)
{return Vector{a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x};
}
bool vis[N][N];
int cnt,n;
Point operator - (Point a,Point b){return Point(a.x-b.x,a.y-b.y,a.z-b.z);}
double get_length(Vector a)
{return sqrt(a.x*a.x+a.y*a.y+a.z*a.z);
}
struct Face
{int v[3];//存储这一个面的三个点的编号，逆时针方向 double area(){return get_length(Cross(p[v[1]]-p[v[0]],p[v[2]]-p[v[0]]))/2;}
}f[N<<2],C[N<<2];
double operator * (Vector a,Vector b){return a.x*b.x+a.y*b.y+a.z*b.z;}
bool see(Face a,Point b)
{return (b-p[a.v[0]])*Cross(p[a.v[1]]-p[a.v[0]],p[a.v[2]]-p[a.v[0]])>0;//进行了微小扰动，不会出现等于0的情况
}
void Convex_3D()
{f[++cnt]={1,2,3},f[++cnt]={3,2,1};//对于前三个点的凸包，我们认为是两张一模一样的平面，但是方向不同//如果{1,2,3}从某一个方向看是顺时针，那么我们就认为我们从另一个方向看//总之可以找到一个方向使得{1,2,3}为逆时针，{3,2,1}是顺时针//于是此时认为这个方向为正方向，{3,2,1}的正方向就是这个方向的反方向 for(int i=4,cc=0;i<=n;i++){bool v;for(int j=1;j<=cnt;j++){v=see(f[j],p[i]);//判断是否看得见 if(!v) C[++cc]=f[j];for(int k=0;k<3;k++) vis[f[j].v[k]][f[j].v[(k+1)%3]]=v;//vis[i][j]用来记录向量ij所表示的面（左手边）是否被看得见 } 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]={x,y,i};//添加新的面，注意三个点的顺序，要保证逆时针 } for(int j=1;j<=cc;j++) f[j]=C[j];cnt=cc;cc=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();//这一步在随机微小扰动，为什么要下文说明 }Convex_3D();double ans=0;for(int i=1;i<=cnt;i++) ans+=f[i].area();printf("%.3f\n",ans);return 0;
}

说明一下为什么要微小扰动：这是为了防止出现四点共面的情况。如果出现了四点共面的情况，可能会导致各种错误。举一个例子，比如对于前四个点共面，那么在加入第四个点的时候，由上面的代码判断的是之前三个点形成的凸包的两个面都看不见，就都会算进去，而且现在不知道新加入的三角形如何剖分（剖分方式很多），所以不行