hdu4449Building Design（三维凸包+平面旋转）
看了几小时也没看懂代码表示的何意。。无奈下来问问考研舍友。
还是考研舍友比较靠谱，分分钟解决了我的疑问。
可能三维的东西在纸面上真的不好表示，网上没有形象的题解，只有简单"明了"的讲解。
这题说起来很简单，求下三维凸包，枚举每一个面，进行坐标旋转，使得当前面作为xoy面时的其他坐标，然后求下投影面的凸包的面积。
为什么要旋转面而不直接算点到面的距离，是因为投影的面积没有办法算。
面旋转时是怎么旋转的，首先求得当前面的法向量p1，再求得它与向量e(0,0,1)的法向量pp,所有的点都是绕pp这个向量旋转的，并且旋转的角度是p1与e的夹角。
如果还有跟我一样不懂空间叉积代表的是什么的同学向后看---》法向量 = 平面上两条不平行的向量的叉积。
这样求出每个点的都是相对于当前平面的点，只取点的x,y坐标即每个点的投影点-->求二维凸包-->求面积。
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 } 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 } 320 return top; 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 }
