增量法求解:
首先任选4个点形成的一个四面体,然后每次新加一个点,分两种情况:
1> 在凸包内,则可以跳过2> 在凸包外,找到从这个点可以"看见"的面S(看不看得见可以用法向量，看点是否在面外侧),删除这些面S,然后对于S的每条边E进行判断，看该点还能否看到这些边E的另一侧的面，这样深度搜索判断。
程序模板：
#include"stdio.h"
#include"string.h"
#include"iostream"
#include"map"
#include"string"
#include"queue"
#include"stack"
#include"vector"
#include"stdlib.h"
#include"algorithm"
#include"math.h"
#define M 533
#define eps 1e-10
#define inf 0x3f3f3f3f
#define mod 1070000009
#define PI acos(-1.0)
using namespace std;
struct node
{
double x,y,z,dis;
node(){}
node(double xx,double yy,double zz):x(xx),y(yy),z(zz){}
node operator +(const node p)//向量间求和操作
{
return node(x+p.x,y+p.y,z+p.z);
}
node operator -(const node p)//向量间相减操作
{
return node(x-p.x,y-p.y,z-p.z);
}
node operator *(const node p)//向量间叉乘操作
{
return node(y*p.z-z*p.y,z*p.x-x*p.z,x*p.y-y*p.x);
}
node operator *(const double p)//向量乘以一个数
{
return node(x*p,y*p,z*p);
}
node operator /(const double p)//向量除以一个数
{
return node(x/p,y/p,z/p);
}
double operator ^(const node p)//向量间点乘操作
{
return x*p.x+y*p.y+z*p.z;
}
};
struct threeD_convex_hull//三维凸包
{
struct face
{
int a,b,c;
int ok;
};
int n;//初始点数
int cnt;//凸包三角形数
node p[M];//初始点
face f[M*8];//凸包三角形
int to[M][M];//点i到j是属于哪个面
double len(node p)//向量的长度
{
return sqrt(p.x*p.x+p.y*p.y+p.z*p.z);
}
double area(node a,node b,node c)//三个点的面积*2
{
return len((b-a)*(c-a));
}
double volume(node a,node b,node c,node d)//四面体面积*6
{
return (b-a)*(c-a)^(d-a);
}
double ptof(node q,face f)//点与面同向
{
node m=p[f.b]-p[f.a];
node n=p[f.c]-p[f.a];
node t=q-p[f.a];
return m*n^t;
}
void dfs(int q,int cur)//维护凸包，若点q在凸包外则更新凸包
{
f[cur].ok=0;//删除当前面，因为此时它在更大的凸包内部
deal(q,f[cur].b,f[cur].a);
deal(q,f[cur].c,f[cur].b);
deal(q,f[cur].a,f[cur].c);
}
//因为每个三角形的的三边是按照逆时针记录的，所以把边反过来后对应的就是与ab边共线的另一个面
void deal(int q,int a,int b)
{
int fa=to[a][b];//与当前面cnt共边的另一个面
face add;
if(f[fa].ok)//若fa面目前是凸包的表面则继续
{
if(ptof(p[q],f[fa])>eps)//若点q能看到fa面继续深搜fa的三条边，更新新的凸包面
dfs(q,fa);
else//当q点可以看到cnt面的同时看不到a，b共边的fa面，则p和a，b点组成一个新的表面三角形
{
add.a=b;
add.b=a;
add.c=q;
add.ok=1;
to[b][a]=to[a][q]=to[q][b]=cnt;
f[cnt++]=add;
}
}
}
int same(int s,int t)//判断两个三角形是否共面
{
node a=p[f[s].a];
node b=p[f[s].b];
node c=p[f[s].c];
if(fabs(volume(a,b,c,p[f[t].a]))<eps
&&fabs(volume(a,b,c,p[f[t].b]))<eps
&&fabs(volume(a,b,c,p[f[t].c]))<eps)
return 1;
return 0;
}
void make()//构建3D凸包
{
cnt=0;
if(n<4)
return;
int sb=1;
for(int i=1;i<n;i++)//保证前两个点不共点
{
if(len(p[0]-p[i])>eps)
{
swap(p[1],p[i]);
sb=0;
break;
}
}
if(sb)return;
sb=1;
for(int i=2;i<n;i++)//保证前三个点不共线
{
if(len((p[1]-p[0])*(p[i]-p[0]))>eps)
{
swap(p[2],p[i]);
sb=0;
break;
}
}
if(sb)return;
sb=1;
for(int i=3;i<n;i++)//保证前四个点不共面
{
if(fabs(volume(p[0],p[1],p[2],p[i]))>eps)
{
swap(p[3],p[i]);
sb=0;
break;
}
}
if(sb)return;
face add;
for(int i=0;i<4;i++)//构建初始四面体
{
add.a=(i+1)%4;
add.b=(i+2)%4;
add.c=(i+3)%4;
add.ok=1;
if(ptof(p[i],add)>eps)
swap(add.c,add.b);
to[add.a][add.b]=to[add.b][add.c]=to[add.c][add.a]=cnt;
f[cnt++]=add;
}
for(int i=4;i<n;i++)//倍增法更新凸包
{
for(int j=0;j<cnt;j++)//判断每个点是在当前凸包的内部或者外部
{
if(f[j].ok&&ptof(p[i],f[j])>eps)//若在外部且看到j面继续
{
dfs(i,j);
break;
}
}
}
int tmp=cnt;//把不是凸包上的面删除即ok=0；
cnt=0;
for(int i=0;i<tmp;i++)
if(f[i].ok)
f[cnt++]=f[i];
}
double Area()//表面积
{
double S=0;
if(n==3)
{
S=area(p[0],p[1],p[2])/2.0;
return S;
}
for(int i=0;i<cnt;i++)
S+=area(p[f[i].a],p[f[i].b],p[f[i].c]);
return S/2.0;
}
double Volume()//体积
{
double V=0;
node mid(0,0,0);
for(int i=0;i<cnt;i++)
V+=volume(p[f[i].a],p[f[i].b],p[f[i].c],mid);
V=fabs(V)/6.0;
return V;
}
int tringleCnt()//凸包表面三角形数目
{
return cnt;
}
int faceCnt()//凸包表面多边形数目
{
int num=0;
for(int i=0;i<cnt;i++)
{
int flag=1;
for(int j=0;j<i;j++)
{
if(same(i,j))
{
flag=0;
break;
}
}
num+=flag;
}
return num;
}
double pf_dis(face f,node q)//点到面的距离
{
double V=volume(p[f.a],p[f.b],p[f.c],q);
double S=area(p[f.a],p[f.b],p[f.c]);
return fabs(V/S);
}
double min_dis(node q)//暴力搜索内部的点q到面的最短距离即体积/面积
{
double mini=inf;
for(int i=0;i<cnt;i++)
{
double h=pf_dis(f[i],q);
if(mini>h)
mini=h;
}
return mini;
}
node barycenter()//凸包的重心
{
node ret(0,0,0),mid(0,0,0);
double sum=0;
for(int i=0;i<cnt;i++)
{
double V=volume(p[f[i].a],p[f[i].b],p[f[i].c],mid);
ret=ret+(mid+p[f[i].a]+p[f[i].b]+p[f[i].c])/4.0*V;
sum+=V;
}
ret=ret/sum;
return ret;
}
}hull;
int main()
{
while(scanf("%d",&hull.n)!=EOF)
{
for(int i=0;i<hull.n;i++)
scanf("%lf%lf%lf",&hull.p[i].x,&hull.p[i].y,&hull.p[i].z);
hull.make();
}
return 0;
}