hdu4273Rescue(三维凸包重心)

链接

模板题已不叫题。。

三维凸包+凸包重心+点到平面距离（体积/点积）  体积-->混合积（先点乘再叉乘）

  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];
 96         fac add;
 97         if (F[f].ok)
 98         {
 99             if (point3of(P[p], F[f]) > eps)
100                 dfs(p, f);
101             else
102             {
103                 add.a = b, add.b = a, add.c = p, add.ok = 1;
104                 to[p][b] = to[a][p] = to[b][a] = cnt;
105                 F[cnt++] = add;
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;
163         fac add;
164         for (int i = 0; i < 4; i++)
165         {
166             add.a = (i+1)%4, add.b = (i+2)%4, add.c = (i+3)%4, add.ok = 1;
167             if (point3of(P[i], add) > 0)
168                 swap(add.b, add.c);
169             to[add.a][add.b] = to[add.b][add.c] = to[add.c][add.a] = cnt;
170             F[cnt++] = add;
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     point3 barycenter()
245     {
246         point3 ans(0,0,0),o(0,0,0);
247         double all=0;
248         for(int i=0;i<cnt;i++)
249         {
250             double vol=volume(o,P[F[i].a],P[F[i].b],P[F[i].c]);
251             ans=ans+(o+P[F[i].a]+P[F[i].b]+P[F[i].c])/4.0*vol;
252             all+=vol;
253         }
254         ans=ans/all;
255         return ans;
256     }
257 
258 }hull;
259 
260 void solve()
261 {
262     double ans = INF;
263     int i;
264     int cnt = hull.cnt;
265     point3 pp = hull.barycenter();
266     for(i = 0 ; i < cnt ; i++)
267     {
268         ans = min(ans,hull.ptoface(pp,i));
269     }
270     printf("%.3f\n",ans);
271 }
272 int main()
273 {
274     int n,i;
275     while(scanf("%d",&n)!=EOF)
276     {
277         hull.n = n;
278         for(i = 0 ; i < n; i++)
279         {
280             scanf("%lf%lf%lf",&pp[i].x,&pp[i].y,&pp[i].z);
281             hull.P[i] = pp[i];
282         }
283         hull.construct();
284         solve();
285     }
286 }
View Code

 

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