題目：Mission Impossible

?

題意：在天花板上有3個攝像頭，下面有一個凸多面體，問攝像頭觀測不到的在地面上的面積是多少。

?

思路：先求出每個攝像頭對于凸多面體在xoy平面的投影，然后求凸包，然后利用半平面交來求面積交即可，注意求投影時用到一個結論：

如果在空間有3點共線則滿足：(z3-z2)/(z2-z1)=(y3-y2)/(y2-y1)=(x3-x2)/(x2-x1)

代碼太亂，有時間再整理。

#include <math.h> #include <stdio.h> #include <iostream> #include <algorithm>using namespace std;const int N=105; const double EPS=1e-8;typedef double DIY;DIY Area;struct Point {DIY x,y;Point() {}Point(DIY _x,DIY _y):x(_x),y(_y){} } p[N];Point MakeVector(Point &P,Point &Q) {return Point(Q.x-P.x,Q.y-P.y); }DIY CrossProduct(Point P,Point Q) {return P.x*Q.y-P.y*Q.x; }Point MinA;Point stack1[N],stack2[N],stack3[N]; int top1,top2,top3;DIY dist(Point A,Point B) {return sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)); }DIY cross(Point A,Point B,Point C) {return (B.x-A.x)*(C.y-A.y)-(B.y-A.y)*(C.x-A.x); }bool cmp1(Point a,Point b) {DIY k=cross(MinA,a,b);if(k>0) return 1;if(k<0) return 0;return dist(MinA,a)>dist(MinA,b); }void Graham1(Point *p,int n,int &top) {int i;for(i=0; i<n; i++)if(p[i].y<p[0].y||(p[i].y==p[0].y&&p[i].x<p[0].x))swap(p[i],p[0]);MinA=p[0];sort(p+1,p+n,cmp1);p[n]=p[0];stack1[0]=p[0];stack1[1]=p[1];stack1[2]=p[2];top=2;for(i=3; i<=n; i++){while(cross(stack1[top-1],stack1[top],p[i])<=0&&top>=2) --top;stack1[++top]=p[i];} }void Graham2(Point *p,int n,int &top) {int i;for(i=0; i<n; i++)if(p[i].y<p[0].y||(p[i].y==p[0].y&&p[i].x<p[0].x))swap(p[i],p[0]);MinA=p[0];sort(p+1,p+n,cmp1);p[n]=p[0];stack2[0]=p[0];stack2[1]=p[1];stack2[2]=p[2];top=2;for(i=3; i<=n; i++){while(cross(stack2[top-1],stack2[top],p[i])<=0&&top>=2) --top;stack2[++top]=p[i];} }void Graham3(Point *p,int n,int &top) {int i;for(i=0; i<n; i++)if(p[i].y<p[0].y||(p[i].y==p[0].y&&p[i].x<p[0].x))swap(p[i],p[0]);MinA=p[0];sort(p+1,p+n,cmp1);p[n]=p[0];stack3[0]=p[0];stack3[1]=p[1];stack3[2]=p[2];top=2;for(i=3; i<=n; i++){while(cross(stack3[top-1],stack3[top],p[i])<=0&&top>=2) --top;stack3[++top]=p[i];} }DIY MultiCross(Point P,Point Q,Point R) {return CrossProduct(MakeVector(Q,P),MakeVector(Q,R)); }struct halfPlane {Point s,t;DIY angle;halfPlane(){}halfPlane(Point _s,Point _t):s(_s),t(_t){}halfPlane(DIY sx,DIY sy,DIY tx,DIY ty):s(sx,sy),t(tx,ty){}void GetAngle(){angle=atan2(t.y-s.y,t.x-s.x);} } hp[N],q[N];Point IntersectPoint(halfPlane P,halfPlane Q) {DIY a1=CrossProduct(MakeVector(P.s,Q.t),MakeVector(P.s,Q.s));DIY a2=CrossProduct(MakeVector(P.t,Q.s),MakeVector(P.t,Q.t));return Point((P.s.x*a2+P.t.x*a1)/(a2+a1),(P.s.y*a2+P.t.y*a1)/(a2+a1)); }bool cmp2(halfPlane P,halfPlane Q) {if(fabs(P.angle-Q.angle)<EPS)return MultiCross(P.s,P.t,Q.s)>0;return P.angle<Q.angle; }bool IsParallel(halfPlane P,halfPlane Q) {return fabs(CrossProduct(MakeVector(P.s,P.t),MakeVector(Q.s,Q.t)))<EPS; }void HalfPlaneIntersect(int n,int &m) {sort(hp,hp+n,cmp2);int i,l=0,r=1;for(m=i=1; i<n; ++i)if(hp[i].angle-hp[i-1].angle>EPS) hp[m++]=hp[i];n=m; m=0;q[0]=hp[0];q[1]=hp[1];for(i=2; i<n; i++){if(IsParallel(q[r],q[r-1])||IsParallel(q[l],q[l+1])) return;while(l<r&&MultiCross(hp[i].s,hp[i].t,IntersectPoint(q[r],q[r-1]))>0) --r;while(l<r&&MultiCross(hp[i].s,hp[i].t,IntersectPoint(q[l],q[l+1]))>0) ++l;q[++r]=hp[i];}while(l<r&&MultiCross(q[l].s,q[l].t,IntersectPoint(q[r],q[r-1]))>0) --r;while(l<r&&MultiCross(q[r].s,q[r].t,IntersectPoint(q[l],q[l+1]))>0) ++l;q[++r]=q[l];for(i=l; i<r; ++i)p[m++]=IntersectPoint(q[i],q[i+1]); }void Solve(Point *p1,Point *p2,Point *p3,int &top1,int &top2,int &top3,int &m) {int i,j;Point a,b;Point O;O.x=O.y=0;int num=0;for(i=0;i<top1;i++){hp[num]=halfPlane(p1[i],p1[(i+1)%top1]);hp[num].GetAngle();num++;}for(i=0;i<top2;i++){hp[num]=halfPlane(p2[i],p2[(i+1)%top2]);hp[num].GetAngle();num++;}for(i=0;i<top3;i++){hp[num]=halfPlane(p3[i],p3[(i+1)%top3]);hp[num].GetAngle();num++;}HalfPlaneIntersect(num,m);Area=0;p[m]=p[0];for(i=0;i<m;++i)Area+=cross(O,p[i],p[i+1]);if(Area<0) Area=-Area;Area/=2.0; }int main() {int n,m,t,i,j,k;Point p1[N],p2[N],p3[N];top1=top2=top3=0;Point tmp[N];DIY x[N],y[N],z[N];while(cin>>n){for(i=0; i<n; i++){cin>>x[i]>>y[i]>>z[i];}for(i=0;i<3;i++)cin>>tmp[i].x>>tmp[i].y;for(i=0;i<n;i++){p1[i].y=y[i]+(0-z[i])*(y[i]-tmp[0].y)/(z[i]-100);p1[i].x=x[i]+(0-z[i])*(x[i]-tmp[0].x)/(z[i]-100);}for(i=0;i<n;i++){p2[i].y=y[i]+(0-z[i])*(y[i]-tmp[1].y)/(z[i]-100);p2[i].x=x[i]+(0-z[i])*(x[i]-tmp[1].x)/(z[i]-100);}for(i=0;i<n;i++){p3[i].y=y[i]+(0-z[i])*(y[i]-tmp[2].y)/(z[i]-100);p3[i].x=x[i]+(0-z[i])*(x[i]-tmp[2].x)/(z[i]-100);}int k=0;Graham1(p1,n,top1);Graham2(p2,n,top2);Graham3(p3,n,top3);Solve(stack1,stack2,stack3,top1,top2,top3,m);printf("%.2lf\n",Area);}return 0; }

?

總結

以上是生活随笔為你收集整理的NJUST4316(立体几何投影的面积交)的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。