题意:一个房间,在房间顶部有三个摄像头,房间有放有一个凸多面体的机器,问这个机器对摄像头在地面造成盲点的面积。房间顶部z=100,地面z=0。
首先根据空间直线方程(x-x1)/(x1-x2)=(y-y1)/(y1-y2)=(z-z1)/(z1-z2) 求出每个摄像头与凸多面体连线与地面(z=0)的交点。在对这个摄像头在地面的交点求凸包,这个区域就个该摄像头盲区的面积,但是需要求三个盲区面积的交,所以把3个凸包的边变成半面,求最后的半面交,得出相交的多边形再求面积即为答案。上来就是甩模板啊。。
#include<cmath> #include<cstdio> #include<iostream> #include<algorithm> using namespace std; const int mm=111; typedef double DIY; struct point { DIY x,y; point() {} point(DIY _x,DIY _y):x(_x),y(_y) {} } g[mm]; 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; } DIY MultiCross(point P,point Q,point R) { return CrossProduct(MakeVector(Q,P),MakeVector(Q,R)); } struct halfPlane { point s,t; double 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[mm<<2],q[mm<<2]; 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 cmp1(halfPlane P,halfPlane Q) { if(fabs(P.angle-Q.angle)<1e-8) 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)))<1e-8; } void HalfPlaneIntersect(int n,int &m) { sort(hp,hp+n,cmp1); int i,l=0,r=1; for(m=i=1; i<n; ++i) if(hp[i].angle-hp[i-1].angle>1e-8)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) g[m++]=IntersectPoint(q[i],q[i+1]); } point data[3][mm],stack[mm],MinA; int top; DIY Direction(point pi,point pj,point pk) //判断向量PiPj在向量PiPk的顺逆时针方向 +顺-逆0共线 { return (pj.x-pi.x)*(pk.y-pi.y)-(pk.x-pi.x)*(pj.y-pi.y); } DIY Dis(point a,point b) { return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)); } bool cmp(point a,point b) { DIY k=Direction(MinA,a,b); if(k>0) return 1; if(k<0) return 0; return Dis(MinA,a)>Dis(MinA,b); } void Graham_Scan(point *a,int numa) { for(int i=0; i<numa; i++) if(a[i].y<a[0].y||(a[i].y==a[0].y&&a[i].x<a[0].x)) swap(a[i],a[0]); MinA=a[0],top=0; sort(a+1,a+numa,cmp); stack[top++]=a[0],stack[top++]=a[1],stack[top++]=a[2]; for(int i=3; i<numa; i++) { while(Direction(stack[top-2],stack[top-1],a[i])<0&&top>=2) top--; stack[top++]=a[i]; } } int main() { int n; double x[mm],y[mm],z[mm],lix[3],liy[3]; while(~scanf("%d",&n)) { int numd=0; for(int i=0; i<n; i++) scanf("%lf%lf%lf",&x[i],&y[i],&z[i]); for(int i=0; i<3; i++) scanf("%lf%lf",&lix[i],&liy[i]); for(int j=0; j<3; j++) for(int i=0; i<n; i++) data[j][i].x=(-100)/(z[i]-100)*(x[i]-lix[j])+lix[j], data[j][i].y=(-100)/(z[i]-100)*(y[i]-liy[j])+liy[j]; int numm=0; for(int j=0; j<3; j++) { Graham_Scan(data[j],n); for(int i=0; i<top; i++) { hp[numm]=halfPlane(stack[i],stack[(i+1)%top]); hp[numm].GetAngle(); numm++; } } int m; double s1=0; point o(0,0); HalfPlaneIntersect(numm,m); for(int i=0; i<m; i++) s1+=Direction(o,g[i],g[(i+1)%m]); s1=fabs(s1)/2; printf("%.2f\n",s1); } return 0; }