三维凸包模板

  1 #include<stdio.h>
  2 #include<algorithm>
  3 #include<string.h>
  4 #include<math.h>
  5 #include<stdlib.h>
  6 using namespace std;
  7 const int MAXN=550;
  8 const double eps=1e-8;
  9 
 10 struct Point
 11 {
 12     double x,y,z;
 13     Point(){}
 14 
 15     Point(double xx,double yy,double zz):x(xx),y(yy),z(zz){}
 16 
 17     //两向量之差
 18     Point operator -(const Point p1)
 19     {
 20         return Point(x-p1.x,y-p1.y,z-p1.z);
 21     }
 22 
 23     //两向量之和
 24     Point operator +(const Point p1)
 25     {
 26         return Point(x+p1.x,y+p1.y,z+p1.z);
 27     }
 28 
 29     //叉乘
 30     Point operator *(const Point p)
 31     {
 32         return Point(y*p.z-z*p.y,z*p.x-x*p.z,x*p.y-y*p.x);
 33     }
 34 
 35     Point operator *(double d)
 36     {
 37         return Point(x*d,y*d,z*d);
 38     }
 39 
 40     Point operator / (double d)
 41     {
 42         return Point(x/d,y/d,z/d);
 43     }
 44 
 45     //点乘
 46     double  operator ^(Point p)
 47     {
 48         return (x*p.x+y*p.y+z*p.z);
 49     }
 50 };
 51 
 52 struct CH3D
 53 {
 54     struct face
 55     {
 56         //表示凸包一个面上的三个点的编号
 57         int a,b,c;
 58         //表示该面是否属于最终凸包上的面
 59         bool ok;
 60     };
 61     //初始顶点数
 62     int n;
 63     //初始顶点
 64     Point P[MAXN];
 65     //凸包表面的三角形数
 66     int num;
 67     //凸包表面的三角形
 68     face F[8*MAXN];
 69     //凸包表面的三角形
 70     int g[MAXN][MAXN];
 71     //向量长度
 72     double vlen(Point a)
 73     {
 74         return sqrt(a.x*a.x+a.y*a.y+a.z*a.z);
 75     }
 76     //叉乘
 77     Point cross(const Point &a,const Point &b,const Point &c)
 78     {
 79         return Point((b.y-a.y)*(c.z-a.z)-(b.z-a.z)*(c.y-a.y),
 80                      (b.z-a.z)*(c.x-a.x)-(b.x-a.x)*(c.z-a.z),
 81                      (b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x)
 82                      );
 83     }
 84     //三角形面积*2
 85     double area(Point a,Point b,Point c)
 86     {
 87         return vlen((b-a)*(c-a));
 88     }
 89     //四面体有向体积*6
 90     double volume(Point a,Point b,Point c,Point d)
 91     {
 92         return (b-a)*(c-a)^(d-a);
 93     }
 94     //正：点在面同向
 95     double dblcmp(Point &p,face &f)
 96     {
 97         Point m=P[f.b]-P[f.a];
 98         Point n=P[f.c]-P[f.a];
 99         Point t=p-P[f.a];
100         return (m*n)^t;
101     }
102     void deal(int p,int a,int b)
103     {
104         int f=g[a][b];//搜索与该边相邻的另一个平面
105         face add;
106         if(F[f].ok)
107         {
108             if(dblcmp(P[p],F[f])>eps)
109               dfs(p,f);
110             else
111             {
112                 add.a=b;
113                 add.b=a;
114                 add.c=p;//这里注意顺序，要成右手系
115                 add.ok=true;
116                 g[p][b]=g[a][p]=g[b][a]=num;
117                 F[num++]=add;
118             }
119         }
120     }
121     void dfs(int p,int now)//递归搜索所有应该从凸包内删除的面
122     {
123          F[now].ok=0;
124          deal(p,F[now].b,F[now].a);
125          deal(p,F[now].c,F[now].b);
126          deal(p,F[now].a,F[now].c);
127     }
128     bool same(int s,int t)
129     {
130         Point &a=P[F[s].a];
131         Point &b=P[F[s].b];
132         Point &c=P[F[s].c];
133         return fabs(volume(a,b,c,P[F[t].a]))<eps &&
134                fabs(volume(a,b,c,P[F[t].b]))<eps &&
135                fabs(volume(a,b,c,P[F[t].c]))<eps;
136     }
137     //构建三维凸包
138     void create()
139     {
140         int i,j,tmp;
141         face add;
142 
143         num=0;
144         if(n<4)return;
145     //**********************************************
146         //此段是为了保证前四个点不共面
147         bool flag=true;
148         for(i=1;i<n;i++)
149         {
150             if(vlen(P[0]-P[i])>eps)
151             {
152                 swap(P[1],P[i]);
153                 flag=false;
154                 break;
155             }
156         }
157         if(flag)return;
158         flag=true;
159         //使前三个点不共线
160         for(i=2;i<n;i++)
161         {
162             if(vlen((P[0]-P[1])*(P[1]-P[i]))>eps)
163             {
164                 swap(P[2],P[i]);
165                 flag=false;
166                 break;
167             }
168         }
169         if(flag)return;
170         flag=true;
171         //使前四个点不共面
172         for(int i=3;i<n;i++)
173         {
174             if(fabs((P[0]-P[1])*(P[1]-P[2])^(P[0]-P[i]))>eps)
175             {
176                 swap(P[3],P[i]);
177                 flag=false;
178                 break;
179             }
180         }
181         if(flag)return;
182     //*****************************************
183         for(i=0;i<4;i++)
184         {
185             add.a=(i+1)%4;
186             add.b=(i+2)%4;
187             add.c=(i+3)%4;
188             add.ok=true;
189             if(dblcmp(P[i],add)>0)swap(add.b,add.c);
190             g[add.a][add.b]=g[add.b][add.c]=g[add.c][add.a]=num;
191             F[num++]=add;
192         }
193         for(i=4;i<n;i++)
194         {
195             for(j=0;j<num;j++)
196             {
197                 if(F[j].ok&&dblcmp(P[i],F[j])>eps)
198                 {
199                     dfs(i,j);
200                     break;
201                 }
202             }
203         }
204         tmp=num;
205         for(i=num=0;i<tmp;i++)
206           if(F[i].ok)
207             F[num++]=F[i];
208 
209     }
210     //表面积
211     double area()
212     {
213         double res=0;
214         if(n==3)
215         {
216             Point p=cross(P[0],P[1],P[2]);
217             res=vlen(p)/2.0;
218             return res;
219         }
220         for(int i=0;i<num;i++)
221           res+=area(P[F[i].a],P[F[i].b],P[F[i].c]);
222         return res/2.0;
223     }
224     double volume()
225     {
226         double res=0;
227         Point tmp(0,0,0);
228         for(int i=0;i<num;i++)
229            res+=volume(tmp,P[F[i].a],P[F[i].b],P[F[i].c]);
230         return fabs(res/6.0);
231     }
232     //表面三角形个数
233     int triangle()
234     {
235         return num;
236     }
237     //表面多边形个数
238     int polygon()
239     {
240         int i,j,res,flag;
241         for(i=res=0;i<num;i++)
242         {
243             flag=1;
244             for(j=0;j<i;j++)
245               if(same(i,j))
246               {
247                   flag=0;
248                   break;
249               }
250             res+=flag;
251         }
252         return res;
253     }
254     //三维凸包重心
255     Point barycenter()
256     {
257         Point ans(0,0,0),o(0,0,0);
258         double all=0;
259         for(int i=0;i<num;i++)
260         {
261             double vol=volume(o,P[F[i].a],P[F[i].b],P[F[i].c]);
262             ans=ans+(o+P[F[i].a]+P[F[i].b]+P[F[i].c])/4.0*vol;
263             all+=vol;
264         }
265         ans=ans/all;
266         return ans;
267     }
268     //点到面的距离
269     double ptoface(Point p,int i)
270     {
271         return fabs(volume(P[F[i].a],P[F[i].b],P[F[i].c],p)/vlen((P[F[i].b]-P[F[i].a])*(P[F[i].c]-P[F[i].a])));
272     }
273 };
274 CH3D hull;
275 int main()
276 {
277     while(scanf("%d",&hull.n) != EOF)
278     {
279         for(int i=0;i<hull.n;i++)
280         {
281             scanf("%lf%lf%lf",&hull.P[i].x,&hull.P[i].y,&hull.P[i].z);
282         }
283         hull.create();
284         printf("%d\n", hull.polygon());
285     }
286     return 0;
287 }