三维凸包模板

分析：给出三维空间中的n个顶点，构建凸包

增量法求解:

首先任选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;
}