题意:在天花板上有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;
}