# hdu4449Building Design（三维凸包+平面旋转）

  1 #include <iostream>
2 #include<cstdio>
3 #include<cstring>
4 #include<algorithm>
5 #include<stdlib.h>
6 #include<vector>
7 #include<cmath>
8 #include<queue>
9 #include<set>
10 using namespace std;
11 #define N 510
12 #define INF 1e20
13 #define max(a,b) (a>b?a:b)
14 #define min(a,b) (a<b?a:b)
15 #define eps 1e-8
16 #define MAXV 505
17 const double pi = acos(-1.0);
18 const double inf = ~0u>>2;
19 //三维点
20 struct point3
21 {
22     double x, y,z;
23     point3() {}
24     point3(double _x, double _y, double _z): x(_x), y(_y), z(_z) {}
25     point3 operator +(const point3 p1)
26     {
27         return point3(x+p1.x,y+p1.y,z+p1.z);
28     }
29     point3 operator - (const point3 p1)
30     {
31         return point3(x - p1.x, y - p1.y, z - p1.z);
32     }
33     point3 operator * (point3 p)
34     {
35         return point3(y*p.z-z*p.y, z*p.x-x*p.z, x*p.y-y*p.x);    //叉乘
36     }
37     point3 operator *(double d)
38     {
39         return point3(x*d,y*d,z*d);
40     }
41     point3 operator /(double d)
42     {
43         return point3(x/d,y/d,z/d);
44     }
45     double operator ^ (point3 p)
46     {
47         return x*p.x+y*p.y+z*p.z;    //点乘
48     }
49
50 } pp[N],rp[N];
51 struct point
52 {
53     double x,y;
54     point(double x=0,double y=0):x(x),y(y) {}
55     point operator -(point b)
56     {
57         return point(x-b.x,y-b.y);
58     }
59 } p[N],ch[N];
60 struct _3DCH
61 {
62     struct fac
63     {
64         int a, b, c;    //表示凸包一个面上三个点的编号
65         bool ok;        //表示该面是否属于最终凸包中的面
66     };
67
68     int n;    //初始点数
69     point3 P[MAXV];    //初始点
70
71     int cnt;    //凸包表面的三角形数
72     fac F[MAXV*8]; //凸包表面的三角形
73
74     int to[MAXV][MAXV];
75     double vlen(point3 a)
76     {
77         return sqrt(a.x*a.x+a.y*a.y+a.z*a.z);
78     }  //向量长度
79     double area(point3 a, point3 b, point3 c)
80     {
81         return vlen((b-a)*(c-a));
82     }    //三角形面积*2
83     double volume(point3 a, point3 b, point3 c, point3 d)
84     {
85         return (b-a)*(c-a)^(d-a);    //四面体有向体积*6
86     }
87     //正：点在面同向
88     double point3of(point3 &p, fac &f)
89     {
90         point3 m = P[f.b]-P[f.a], n = P[f.c]-P[f.a], t = p-P[f.a];
91         return (m * n) ^ t;
92     }
93     void deal(int p, int a, int b)
94     {
95         int f = to[a][b];
97         if (F[f].ok)
98         {
99             if (point3of(P[p], F[f]) > eps)
100                 dfs(p, f);
101             else
102             {
104                 to[p][b] = to[a][p] = to[b][a] = cnt;
106             }
107         }
108     }
109     void dfs(int p, int cur)
110     {
111         F[cur].ok = 0;
112         deal(p, F[cur].b, F[cur].a);
113         deal(p, F[cur].c, F[cur].b);
114         deal(p, F[cur].a, F[cur].c);
115     }
116     bool same(int s, int t)
117     {
118         point3 &a = P[F[s].a], &b = P[F[s].b], &c = P[F[s].c];
119         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;
120     }
121     //构建三维凸包
122     void construct()
123     {
124         cnt = 0;
125         if (n < 4)
126             return;
127         bool sb = 1;
128         //使前两点不公点
129         for (int i = 1; i < n; i++)
130         {
131             if (vlen(P[0] - P[i]) > eps)
132             {
133                 swap(P[1], P[i]);
134                 sb = 0;
135                 break;
136             }
137         }
138         if (sb)return;
139         sb = 1;
140         //使前三点不公线
141         for (int i = 2; i < n; i++)
142         {
143             if (vlen((P[0] - P[1]) * (P[1] - P[i])) > eps)
144             {
145                 swap(P[2], P[i]);
146                 sb = 0;
147                 break;
148             }
149         }
150         if (sb)return;
151         sb = 1;
152         //使前四点不共面
153         for (int i = 3; i < n; i++)
154         {
155             if (fabs((P[0] - P[1]) * (P[1] - P[2]) ^ (P[0] - P[i])) > eps)
156             {
157                 swap(P[3], P[i]);
158                 sb = 0;
159                 break;
160             }
161         }
162         if (sb)return;
164         for (int i = 0; i < 4; i++)
165         {
167             if (point3of(P[i], add) > 0)
171         }
172         for (int i = 4; i < n; i++)
173         {
174             for (int j = 0; j < cnt; j++)
175             {
176                 if (F[j].ok && point3of(P[i], F[j]) > eps)
177                 {
178                     dfs(i, j);
179                     break;
180                 }
181             }
182         }
183         int tmp = cnt;
184         cnt = 0;
185         for (int i = 0; i < tmp; i++)
186         {
187             if (F[i].ok)
188             {
189                 F[cnt++] = F[i];
190             }
191         }
192     }
193     //表面积
194     double area()
195     {
196         double ret = 0.0;
197         for (int i = 0; i < cnt; i++)
198         {
199             ret += area(P[F[i].a], P[F[i].b], P[F[i].c]);
200         }
201         return ret / 2.0;
202     }
203     double ptoface(point3 p,int i)
204     {
205         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])));
206     }
207     //体积
208     double volume()
209     {
210         point3 O(0, 0, 0);
211         double ret = 0.0;
212         for (int i = 0; i < cnt; i++)
213         {
214             ret += volume(O, P[F[i].a], P[F[i].b], P[F[i].c]);
215         }
216         return fabs(ret / 6.0);
217     }
218     //表面三角形数
219     int facetCnt_tri()
220     {
221         return cnt;
222     }
223
224     //表面多边形数
225     int facetCnt()
226     {
227         int ans = 0;
228         for (int i = 0; i < cnt; i++)
229         {
230             bool nb = 1;
231             for (int j = 0; j < i; j++)
232             {
233                 if(same(i, j))
234                 {
235                     nb = 0;
236                     break;
237                 }
238             }
239             ans += nb;
240         }
241         return ans;
242     }
243
244 } hull;
245 point3 get_point(point3 st,point3 ed,point3 tp)//tp在直线st-en上的垂足
246 {
247     double t1=(tp-st)^(ed-st);
248     double t2=(ed-st)^(ed-st);
249     double t=t1/t2;
250     point3 tt = (ed-st)*t;
251     point3 ans=st + tt;
252     return ans;
253 }
254 point3 rotate(point3 st,point3 ed,point3 tp,double A)//将点tp绕st-ed逆时针旋转A度  从ed往st看为逆时针
255 {
256     point3 root=get_point(st,ed,tp);
257     point3 e=(ed-st)/hull.vlen(ed-st);
258     point3 r=tp-root;
259     point3 vec=e*r;
260     point3 ans=r*cos(A)+vec*sin(A)+root;
261     return ans;
262 }
263 int dcmp(double x)
264 {
265     if(fabs(x)<eps) return 0;
266     return x<0?-1:1;
267 }
268 double cross(point a,point b)
269 {
270     return a.x*b.y-a.y*b.x;
271 }
272 double mul(point p0,point p1,point p2)
273 {
274     return cross(p1-p0,p2-p0);
275 }
276 double dis(point a)
277 {
278     return sqrt(a.x*a.x+a.y*a.y);
279 }
280 bool cmp(point a,point b)
281 {
282     if(dcmp(mul(p[0],a,b))==0)
283         return dis(a-p[0])<dis(b-p[0]);
284     else
285         return dcmp(mul(p[0],a,b))>0;
286 }
287 double Polyarea(int n,point p[])
288 {
289     double area = 0;
290     for(int i = 1 ; i < n-1 ; i++)
291         area+=cross(p[i]-p[0],p[i+1]-p[0]);
292     return fabs(area)/2;
293 }
294 int Graham(int n)
295 {
296     int i,k = 0,top;
297     point tmp;
298     for(i = 0 ; i < n; i++)
299     {
300         if(p[i].y<p[k].y||(p[i].y==p[k].y&&p[i].x<p[k].x))
301             k = i;
302     }
303     if(k!=0)
304     {
305         tmp = p[0];
306         p[0] = p[k];
307         p[k] = tmp;
308     }
309     sort(p+1,p+n,cmp);
310     ch[0] = p[0];
311     ch[1] = p[1];
312     top = 1;
313     for(i = 2; i < n ; i++)
314     {
315         while(top>0&&dcmp(mul(ch[top-1],ch[top],p[i]))<0)
316             top--;
317         top++;
318         ch[top] = p[i];
319     }
321 }
322 void solve()
323 {
324     int i,j;
325     double h = 0,sa = INF;
326     int cnt = hull.cnt,n = hull.n;
327     for(i = 0; i < cnt ; i++)
328     {
329         for(j = 0 ; j < n; j++)
330             rp[j] = hull.P[j];//rp数组为待旋转点
331
332         point3 p1 = (rp[hull.F[i].b]-rp[hull.F[i].a])*(rp[hull.F[i].c]-rp[hull.F[i].a]);//平面的法向量，注意法向量的方向。
333         point3 e = point3(0,0,1);//z轴上的向量
334         point3 vec = p1*e;//垂直于p1和e的向量，即所有点将要绕其旋转的向量
335
336         double A = p1^e/hull.vlen(p1);
337         A = acos(A); //p1与e的夹角
338
339         if(dcmp(A)!=0&&dcmp(A-pi)!=0)
340         {
341             point3 p0 = point3(0,0,0);
342             for(j = 0 ; j < n; j++)
343                 rp[j] = rotate(p0,vec,rp[j],A);//绕直线p0-vec旋转
344         }
345         double tt = rp[hull.F[i].a].z;
346         for(j = 0 ; j < n; j++)
347             rp[j].z-=tt;
348         double th = 0,ts;
349         for(j = 0 ; j < n; j++)
350         {
351             th = max(th,hull.ptoface(hull.P[j],i));
352         }
353         for(j = 0 ; j< n; j++)
354         {
355             p[j].x = rp[j].x;
356             p[j].y = rp[j].y;
357         }
358         int m = Graham(n);
359         ch[++m] = ch[0];
360
361         ts = Polyarea(m,ch);//cout<<ts<<endl;
362         if(dcmp(th-h)>0||(dcmp(th-h)==0&&dcmp(ts-sa)<0))
363         {
364             h = th;
365             sa = ts;
366         }
367     }
368     printf("%.3f %.3f\n",h,sa);
369 }
370 int main()
371 {
372     int n,i;
373     while(scanf("%d",&n)&&n)
374     {
375         hull.n = n;
376         for(i = 0 ; i < n; i++)
377         {
378             scanf("%lf%lf%lf",&pp[i].x,&pp[i].y,&pp[i].z);
379             hull.P[i] = pp[i];
380         }
381         hull.construct();
382         solve();
383     }
384 }
View Code

posted @ 2014-09-04 18:40  _雨  阅读(373)  评论(0编辑  收藏  举报