三维凸包

来源:互联网 发布:百度的人工智能 编辑:程序博客网 时间:2024/05/12 06:20

三维凸包模板 hdu4266

#include <iostream>#include <cstdio>#include <cstdlib>#include <cstring>#include <cmath>#include <algorithm>using namespace std;const int MAXN = 1100;const double eps = 1e-9;struct Point{     double x,y,z;     Point(){}     Point(double _x,double _y,double _z):x(_x),y(_y),z(_z){}     Point operator -(const Point &P)     {           return Point(x-P.x, y-P.y, z-P.z);     }       Point operator ^ (const Point &P)//叉积      {           return Point(y*P.z - z * P.y, z * P.x - x * P.z, x*P.y - y * P.x);     }     double operator * (const Point &P)//点积      {           return (x * P.x + y * P.y + z * P.z);      }      };struct CH3D{       struct face       {            int a,b,c;            bool ok ;// 判断是否是凸包的面             face(){}            face(int _a,int _b,int _c, bool _ok):a(_a),b(_b),c(_c),ok(_ok){}            void init(int _a,int _b,int _c, bool _ok){                 a = _a; b = _b; c =_c; ok = _ok;            }       };       int n,num;       Point p[MAXN];        int g[MAXN][MAXN];//凸包表面的三角形       face F[MAXN * 8];       double vlen(const Point &a) //向量长度       {              return sqrt(a.x*a.x+a.y*a.y+a.z*a.z);       }       double area(Point &a,Point &b,Point &c)// 三角形面积*2       {              return vlen((b-a) ^ (c-a));       }        double volume(Point &a, Point &b , Point &c ,Point &d)//四面体体积       {              return ((b-a)^(c-a)) * (d-a);       }          double dblcmp(Point &P,face &f)                                           {              Point m=p[f.b]-p[f.a];              Point n=p[f.c]-p[f.a];              Point t=P-p[f.a];              return (m ^ n) * t;       }         void deal(int P,int a,int b)       {            int f=g[a][b];               if(F[f].ok)            {                 if(dblcmp(p[P],F[f])>eps) dfs(P,f);                 else                 {                     face add(b,a,P,1);                      g[P][b]=g[a][P]=g[b][a]=num;                     F[num++]=add;                 }            }       }       void dfs(int p,int now)       {            F[now].ok=0;            deal(p,F[now].b,F[now].a);            deal(p,F[now].c,F[now].b);            deal(p,F[now].a,F[now].c);       }       bool same(int s,int t)       {            Point &a=p[F[s].a];            Point &b=p[F[s].b];            Point &c=p[F[s].c];            return 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;       }       void solve()//构建三维凸包       {            face add; bool flag=true;            num=0;            if(n<4) return;            for(int i=1;i<n;i++)//此段是为了保证前四个点不共面,若以保证,则可去掉            {                if(vlen(p[0]-p[i])>eps){                       swap(p[1],p[i]);                       flag=false; break;                }            }            if(flag)   return;            flag=true;            for(int i=2;i<n;i++)//使前三点不共线            {                 if(vlen((p[0]-p[1])^(p[1]-p[i]))>eps)                 {                       swap(p[2],p[i]);                       flag=false; break;                 }            }            if(flag)   return;            flag=true;            for(int i=3;i<n;i++)//使前四点不共面            {                  if(fabs(((p[0]-p[1])^(p[1]-p[2]))*(p[0]-p[i]))>eps){                        swap(p[3],p[i]);                        flag=false; break;                  }            }            if(flag) return;            for(int i=0;i<4;i++)            {                   add.init((i+1)%4,(i+2)%4,(i+3)%4,true);                   if(dblcmp(p[i],add)>0) swap(add.b,add.c);                   g[add.a][add.b]=g[add.b][add.c]=g[add.c][add.a]=num;                   F[num++]=add;            }            for(int i=4;i<n;i++){                for(int j=0;j<num;j++) {                     if(F[j].ok && dblcmp(p[i],F[j])>eps){                          dfs(i,j); break;                     }                }            }            int tmp=num; num = 0;            for(int i=0;i<tmp;i++)                    if(F[i].ok) F[num++]=F[i];       }       double area()//表面积       {              double res=0.0;              if(n==3){                   //Point p=cross(P[0],P[1],P[2]);                   res=vlen((p[1]-p[0])^(p[2]-p[0]))/2.0;                   return res;              }                      for(int i=0;i<num;i++)                 res+=vlen((p[F[i].b]-p[F[i].a]) ^ (p[F[i].c]-p[F[i].a]));              return res/2.0;       }       double volume() //体积       {              double res=0.0;              Point temp(0,0,0);              for(int i=0;i<num;i++)                 res+=((p[F[i].b]-p[F[i].a])^(p[F[i].c]-p[F[i].a])) * (temp-p[F[i].a]);               return fabs(res/6.0);       }       int triangle()//表面三角形个数        {              return num;       }       int polygon()//表面多边形个数       {           int res = 0;           for(int i=0,flag = 1; i<num;i++)           {                for(int j=0;j<i;j++)                if(same(i,j)){                      flag=0;break;                 }                res+=flag;           }           return res;       }       double getMinDis(Point Q)       {              double ans = 0x7fffffff;              for (int i = 0; i < num;++i){              ans = min(ans,              volume(Q,p[F[i].a],p[F[i].b],p[F[i].c]) /area(p[F[i].a],p[F[i].b],p[F[i].c]));              }              return ans;       }       }hull; int main(){    while (scanf("%d",&hull.n)!=EOF)    {          if (hull.n ==0) break;          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.solve();          int q; cin>>q;          for (int i = 0; i < q;++i)          {              Point Q;              scanf("%lf%lf%lf",&Q.x,&Q.y,&Q.z);              printf("%.4lf\n", hull.getMinDis(Q));          }    }    return 0;}


2013 ACM/ICPC 长沙邀请赛- I 题   Throw the stones

#include <iostream>#include <cstdio>#include <cstdlib>#include <cstring>#include <cmath>#include <algorithm>#include <map>using namespace std;const int MAXN = 10002;const double eps = 1e-9;int child;double last,D;struct Point{     double x,y,z;     Point(){}     Point(double _x,double _y,double _z):x(_x),y(_y),z(_z){}     Point operator -(const Point &P)     {           return Point(x-P.x, y-P.y, z-P.z);     }       Point operator ^ (const Point &P)//叉积      {           return Point(y*P.z - z * P.y, z * P.x - x * P.z, x*P.y - y * P.x);     }     double operator * (const Point &P)//点积      {           return (x * P.x + y * P.y + z * P.z);      }      };struct CH3D{       struct face       {            int a,b,c;            bool ok ;// 判断是否是凸包的面             face(){}            face(int _a,int _b,int _c, bool _ok):a(_a),b(_b),c(_c),ok(_ok){}            void init(int _a,int _b,int _c, bool _ok){                 a = _a; b = _b; c =_c; ok = _ok;            }       };       int n,num;       Point p[MAXN];        map<int,int> g[MAXN];//凸包表面的三角形       face F[MAXN * 10];       double vlen(const Point &a) //向量长度       {              return sqrt(a.x*a.x+a.y*a.y+a.z*a.z);       }       double area(Point &a,Point &b,Point &c)// 三角形面积*2       {              return vlen((b-a) ^ (c-a));       }        double volume(Point &a, Point &b , Point &c ,Point &d)//四面体体积       {              return ((b-a)^(c-a)) * (d-a);       }          double dblcmp(Point &P,face &f)                                           {              Point m=p[f.b]-p[f.a];              Point n=p[f.c]-p[f.a];              Point t=P-p[f.a];              return (m ^ n) * t;       }         void deal(int P,int a,int b)       {            int f=g[a][b];               if(F[f].ok)            {                 if(dblcmp(p[P],F[f])>eps) dfs(P,f);                 else                 {                     D += (p[a]^p[P]) * p[b];                       face add(b,a,P,1);                      g[P][b]=g[a][P]=g[b][a]=num;                     F[num++]=add;                 }            }       }      void dfs(int P,int now)       {            D-=(p[F[now].b] ^ p[F[now].c] )* p[F[now].a];            F[now].ok=0;            deal(P,F[now].b,F[now].a);            deal(P,F[now].c,F[now].b);            deal(P,F[now].a,F[now].c);       }       bool same(int s,int t)       {            Point &a=p[F[s].a];            Point &b=p[F[s].b];            Point &c=p[F[s].c];            return 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;       }       void init()       {            face add;            for(int i=0;i<4;i++)            {                   add.init((i+1)%4,(i+2)%4,(i+3)%4,true);                   if(dblcmp(p[i],add)>0) swap(add.b,add.c);                   g[add.a][add.b]=g[add.b][add.c]=g[add.c][add.a]=num;                   F[num++]=add;            }       }       void update(int i)       {            for(int j=0;j<num;j++)             if(F[j].ok && dblcmp(p[i],F[j])>eps){                dfs(i,j); break;            }       }       void solve()//构建三维凸包       {            face add; bool flag=true;            num=0;            if(n<4) return;            for(int i=1;i<n;i++)//此段是为了保证前四个点不共面,若以保证,则可去掉            {                if(vlen(p[0]-p[i])>eps){                       swap(p[1],p[i]);                       flag=false; break;                }            }            if(flag)   return;            flag=true;            for(int i=2;i<n;i++)//使前三点不共线            {                 if(vlen((p[0]-p[1])^(p[1]-p[i]))>eps)                 {                       swap(p[2],p[i]);                       flag=false; break;                 }            }            if(flag)   return;            flag=true;            for(int i=3;i<n;i++)//使前四点不共面            {                  if(fabs(((p[0]-p[1])^(p[1]-p[2]))*(p[0]-p[i]))>eps){                        swap(p[3],p[i]);                        child = i; init();                        for (int j = 4; j<=child;++j) update(j);                        last = volume();                        flag=false; break;                  }            }            if(flag) return;            for(int i=child+1;i<n;i++) {                    D = 0;update(i);D/=6.0;                    if (D>last+eps){                       child = i;                       last = D;                    }            }            int tmp=num; num = 0;            for(int i=0;i<tmp;i++)                    if(F[i].ok) F[num++]=F[i];       }       double area()//表面积       {              double res=0.0;              if(n==3){                   //Point p=cross(P[0],P[1],P[2]);                   res=vlen((p[1]-p[0])^(p[2]-p[0]))/2.0;                   return res;              }                      for(int i=0;i<num;i++)                 res+=vlen((p[F[i].b]-p[F[i].a]) ^ (p[F[i].c]-p[F[i].a]));              return res/2.0;       }       double volume() //体积       {              double res=0.0;              Point temp(0,0,0);              for(int i=0;i<num;i++)                if (F[i].ok)                 res+=((p[F[i].b]-p[F[i].a])^(p[F[i].c]-p[F[i].a])) * (temp-p[F[i].a]);                  //res +=(p[F[i].b] ^ p[F[i].c]) *p[F[i].a];              return fabs(res/6.0);       }       int triangle()//表面三角形个数        {              return num;       }       int polygon()//表面多边形个数       {           int res = 0;           for(int i=0,flag = 1; i<num;i++)           {                for(int j=0;j<i;j++)                if(same(i,j)){                      flag=0;break;                 }                res+=flag;           }           return res;       }       double getMinDis(Point Q)       {              double ans = 0x7fffffff;              for (int i = 0; i < num;++i){              ans = min(ans,              volume(Q,p[F[i].a],p[F[i].b],p[F[i].c]) /area(p[F[i].a],p[F[i].b],p[F[i].c]));              }              return ans;       }       }hull; int main(){    int test =0;    while (scanf("%d",&hull.n)!=EOF)    {          test ++;          child = 0;last = 0;          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.solve();          printf("Case #%d:\n%d %.2lf\n",test, child+1,last);    }    return 0;}