首页 > 代码库 > hdu4273Rescue(三维凸包重心)

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             else102             {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

 

hdu4273Rescue(三维凸包重心)