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