POJ 3862 POJ 3528 HDU 3662 HDU 4273 三維凸包問題解決模板

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