POJ4449(三维凸包+空间坐标旋转+二维凸包)
生活随笔
收集整理的這篇文章主要介紹了
POJ4449(三维凸包+空间坐标旋转+二维凸包)
小編覺得挺不錯(cuò)的,現(xiàn)在分享給大家,幫大家做個(gè)參考.
題目:Building Design
?
題意:
題目就是給了一個(gè)凸多面題。 要把這個(gè)凸多面體放到地上,一個(gè)面要接觸到地面。 求一種放法,使得最高點(diǎn)最高,最高點(diǎn)相同時(shí)使投影面積最小。輸出最高點(diǎn)高度H和投影面積S。 先三維凸包一下求得所有的三角形表面。 然后枚舉每一個(gè)表面,進(jìn)行坐標(biāo)變換,旋轉(zhuǎn)使得這個(gè)表面變水平,然后投影求投影面積,算最高點(diǎn)高度。 ? 題目的數(shù)據(jù)好像沒有n==1 ,2,3的情況,可以忽略。 有個(gè)問題就是如果按照旋轉(zhuǎn)之后的z坐標(biāo)值來算高度的話,會(huì)產(chǎn)生很大的精度差。最好就是直接算點(diǎn)到面的距離。 //============================================================================ // 三維凸包+二維凸包+三維坐標(biāo)旋轉(zhuǎn)變換 //============================================================================#include <iostream> #include <string.h> #include <math.h> #include <algorithm> #include <stdio.h> #include <vector> #include <map> #include <set> #include <queue> using namespace std; const int MAXN = 110; const double eps=1e-8; const double PI=acos(-1.0); inline int sign(double p) {if(p>eps)return 1;else if(p<-eps)return -1;else return 0; } struct Point {double x,y,z;Point(double xx=0,double yy=0,double zz=0):x(xx),y(yy),z(zz){}Point operator +(const Point p1){return Point(x+p1.x,y+p1.y,z+p1.z);}Point operator -(const Point p1){return Point(x-p1.x,y-p1.y,z-p1.z);}Point operator *(const Point p1){return Point(y*p1.z-z*p1.y,z*p1.x-x*p1.z,x*p1.y-y*p1.x);}Point operator *(double d){return Point(d*x,d*y,d*z);}Point operator /(double d){return Point(x/d,y/d,z/d);}double operator ^(const Point p1){return x*p1.x+y*p1.y+z*p1.z;}double getLen(){return sqrt(x*x+y*y+z*z);}void read(){scanf("%lf%lf%lf",&x,&y,&z);} }; struct face {int a,b,c;bool ok; };struct CH3D {int n;//初始頂點(diǎn)數(shù)Point P[MAXN];//初始頂點(diǎn)int num;//凸包表面的三角形個(gè)數(shù)face F[8*MAXN];//凸包表面的三角形int g[MAXN][MAXN];double vlen(Point a){return sqrt(a.x*a.x+a.y*a.y+a.z*a.z);}Point cross(Point a,Point b,Point c){return Point((b.y-a.y)*(c.z-a.z)-(b.z-a.z)*(c.y-a.y),(b.z-a.z)*(c.x-a.x)-(b.x-a.x)*(c.z-a.z),(b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x));}//三角形面積*2double area(Point a,Point b,Point c){return vlen((b-a)*(c-a));}//四面體面積*6double volume(Point a,Point b,Point c,Point d){return ((b-a)*(c-a))^(d-a);}double dblcmp(Point &p,face &f){Point m=P[f.b]-P[f.a];Point n=P[f.c]-P[f.a];Point t=p-P[f.a];return (m*n)^t;}void deal(int p,int a,int b){int f=g[a][b];face add;if(F[f].ok){if(dblcmp(P[p],F[f])>eps)dfs(p,f);else{add.a=b;add.b=a;add.c=p;add.ok=true;g[p][b]=g[a][p]=g[b][a]=num;F[num++]=add;}}}void dfs(int p,int now){F[now].ok=false;deal(p,F[now].b,F[now].a);deal(p,F[now].c,F[now].b);deal(p,F[now].a,F[now].c);}bool same(int s,int t){Point &a=P[F[s].a];Point &b=P[F[s].b];Point &c=P[F[s].c];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;}void create(){int i,j,tmp;face add;num=0;if(n<4)return;bool flag=true;for(i=1;i<n;i++){if(vlen(P[0]-P[i])>eps){swap(P[1],P[i]);flag=false;break;}}if(flag)return;flag=true;for(i=2;i<n;i++){if(vlen((P[0]-P[1])*(P[1]-P[i]))>eps){swap(P[2],P[i]);flag=false;break;}}if(flag)return;flag=true;for(i=3;i<n;i++){if(fabs((P[0]-P[1])*(P[1]-P[2])^(P[0]-P[i]))>eps){swap(P[3],P[i]);flag=false;break;}}if(flag)return;for(i=0;i<4;i++){add.a=(i+1)%4;add.b=(i+2)%4;add.c=(i+3)%4;add.ok=true;if(dblcmp(P[i],add)>0)swap(add.b,add.c);g[add.a][add.b]=g[add.b][add.c]=g[add.c][add.a]=num;F[num++]=add;}for(i=4;i<n;i++){for(j=0;j<num;j++){if(F[j].ok && dblcmp(P[i],F[j])>eps){dfs(i,j);break;}}}tmp=num;for(i=num=0;i<tmp;i++)if(F[i].ok)F[num++]=F[i];}double ptoface(Point p,int i){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])));}}; CH3D hull; Point ps[MAXN]; int n; inline Point get_point(Point st,Point ed,Point tp) {double t1=(tp-st)^(ed-st);double t2=(ed-st)^(ed-st);double t=t1/t2;Point ans=st + ((ed-st)*t);return ans; } inline double dist(Point st,Point ed) {return sqrt((ed-st)^(ed-st)); } Point rotate(Point st,Point ed,Point tp,double A) {Point root=get_point(st,ed,tp);Point e=(ed-st)/dist(ed,st);Point r=tp-root;Point vec=e*r;Point ans=r*cos(A)+vec*sin(A)+root;return ans; }struct point {double x,y; }; point list[MAXN]; int stack[MAXN],top; double cross(point p0,point p1,point p2) {return (p1.x-p0.x)*(p2.y-p0.y)-(p1.y-p0.y)*(p2.x-p0.x); } double dis(point p1,point p2) {return sqrt((p2.x-p1.x)*(p2.x-p1.x)+(p2.y-p1.y)*(p2.y-p1.y)); } bool cmp(point p1,point p2) {double tmp=cross(list[0],p1,p2);if(sign(tmp)>0)return true;else if(sign(tmp)==0 && dis(list[0],p1)<dis(list[0],p2))return true;else return false; } void init(int n) {point p0;p0=list[0];int k=0;for(int i=1;i<n;i++){if( (p0.y+eps > list[i].y) || ( sign(p0.y-list[i].y)==0 && (p0.x+eps > list[i].x))){p0=list[i];k=i;}}list[k]=list[0];list[0]=p0;sort(list+1,list+n,cmp); } void graham(int n) {init(n);if(n==1){top=0;stack[0]=0;return;}if(n==2){top=1;stack[0]=0;stack[1]=1;return;}for(int i=0;i<=1;i++)stack[i]=i;top=1;for(int i=2;i<n;i++){while(top>0 && cross(list[stack[top-1]],list[stack[top]],list[i])<=eps)top--;top++;stack[top]=i;} }bool input() {scanf("%d",&n);if(n==0)return false;hull.n=n;for(int i=0;i<n;i++)hull.P[i].read();return true; }void solve() {double ansH=0,ansS=1e200;double H,S;if(n<=2){printf("0.000 0.000\n");return;}if(n==3){ansH=0;ansS=(hull.P[1]-hull.P[0])^(hull.P[2]-hull.P[0]);ansS/=2.0;printf("%.3lf %.3lf\n",ansH,ansS);return;}hull.create();for(int i=0;i< hull.num;i++){for(int j=0;j<n;j++)ps[j]=hull.P[j];Point p1=(ps[hull.F[i].b]-ps[hull.F[i].a])*(ps[hull.F[i].c]-ps[hull.F[i].a]);Point e(0,0,1);Point vec=p1*e;double A=p1^e/p1.getLen();A=acos(A);if(sign(A)!=0 && sign(A-PI)!=0){Point p0(0,0,0);for(int j=0;j<n;j++)ps[j]=rotate(p0,vec,ps[j],A);}double tt=ps[hull.F[i].a].z;for(int j=0;j<n;j++)ps[j].z -= tt;H=0;for(int j=0;j<n;j++){H=max(H,hull.ptoface(hull.P[j],i));}for(int j=0;j<n;j++){list[j].x=ps[j].x;list[j].y=ps[j].y;}graham(n);S=0;for(int j=0;j<top;j++)S+=(list[stack[j]].x*list[stack[j+1]].y - list[stack[j]].y*list[stack[j+1]].x);S+=(list[stack[top]].x*list[stack[0]].y - list[stack[top]].y*list[stack[0]].x);S=fabs(S);S/=2;if(H > ansH || ( sign(H-ansH)==0 && ansS > S )){ansH=H;ansS=S;}}printf("%.3lf %.3lf\n",ansH,ansS); }int main() {while(input())solve();return 0; }
?
總結(jié)
以上是生活随笔為你收集整理的POJ4449(三维凸包+空间坐标旋转+二维凸包)的全部內(nèi)容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: HDU4082(相似三角形的个数)
- 下一篇: UVA10173(求凸包的面积最小外接矩