首页 > 代码库 > 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 }
hdu4273Rescue(三维凸包重心)
声明:以上内容来自用户投稿及互联网公开渠道收集整理发布,本网站不拥有所有权,未作人工编辑处理,也不承担相关法律责任,若内容有误或涉及侵权可进行投诉: 投诉/举报 工作人员会在5个工作日内联系你,一经查实,本站将立刻删除涉嫌侵权内容。