hdu - 4316 - Mission Impossible - 计算几何

题目模型可以抽象为：有三个点光源，一个放在z = 0 平面的凸多面体，求三个点光源照射凸多面体在z = 0平面上的影子的公共部分。

对于每个点光源，它在z = 0平面上产生的影子都是一个凸多边形，三个凸多边形的交集即为所求。

凸多边形求法：每个点光源和凸多面体的每个点进行连线，和z = 0平面求交点，对于交点求凸包，即为所求凸多边形。

三个凸多边形求交集：半平面交http://acm.hdu.edu.cn/showproblem.php?pid=4316

#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;
const int MAXN = 1001;
const double eps = 1e-12;

int dblcmp(double x)
{
if(fabs(x) < eps)
{
return 0;
}
return x > 0 ? 1 : -1;
}

struct Point
{
double x, y, z;
inline void input2d()
{
scanf("%lf%lf",&x,&y);
z = 100.0;
}
inline void input3d()
{
scanf("%lf%lf%lf",&x,&y,&z);
}
bool operator < (const Point &point) const
{
if(y == point.y)
{
return x < point.x;
}
return y < point.y;
}
}machine[MAXN], light[3], shadow[3][MAXN], final[MAXN * 3];

int n, shadowNumber[3], shadowPoly[3][MAXN], shadowPolyTop[3], finalNumber;

double operator * (const Point &x, const Point &y)
{
return x.x * y.y - x.y * y.x;
}

Point operator - (const Point &x, const Point &y)
{
Point r;
r.x = x.x - y.x;
r.y = x.y - y.y;
return r;
}

bool operator == (const Point &a, const Point &b)
{
return fabs(a.x - b.x) < eps && fabs(a.y - b.y) < eps;
}

struct Line
{
Point a, b;
double ang;
}line[MAXN*3], stack[MAXN*3];
int lineNumber, stackTop;

bool operator < (const Line &x, const Line &y)
{
if(fabs(x.ang - y.ang) < eps)
{
return (y.b - x.a) * (x.b - y.a) > eps;
}
return x.ang < y.ang;
}

Point operator * (const Line &x, const Line &y)
{
double a1 = (y.b - x.a) * (y.a - x.a);
double a2 = (y.a - x.b) * (y.b - x.b);
Point r;
r.x = (x.a.x * a2 + x.b.x * a1) / (a1 + a2);
r.y = (x.a.y * a2 + x.b.y * a1) / (a1 + a2);
return r;
}

bool mult(const Point &s, const Point &e, const Point &o)
{
return (s.x - o.x) * (e.y - o.y) >= (e.x - o.x) * (s.y - o.y);
}

void graham(Point p[], int n, int res[], int &top)
{
int len;
top = 1;
sort(p, p + n);
if(n == 0) return;
res[0] = 0;
if(n == 1) return;
res[1] = 1;
if(n == 2) return;
res[2] = 2;
for(int i=2;i<n;++i)
{
while(top && mult(p[i], p[res[top]], p[res[top-1]]))
{
-- top;
}
res[++top] = i;
}
len = top;
res[++top] = n - 2;
for(int i=n-3;i>=0;--i)
{
while(top!=len && mult(p[i], p[res[top]], p[res[top-1]]))
{
-- top;
}
res[++top] = i;
}
}

double cross(Point &a, Point &b, Point &o)
{
return (a.x - o.x) * (b.y - o.y) - (a.y - o.y) * (b.x - o.x);
}

bool judgeOut(const Line &x, const Point &p)
{
return (p - x.a) * (x.b - x.a) > eps;
}

bool isParellel(const Line &x, const Line &y)
{
return fabs((x.b - x.a) * (y.b - y.a)) < eps;
}

double getArea(Point p[], int n)
{
double ans = 0.0;
for(int i=2;i<n;++i)
{
ans += fabs(cross(p[i], p[i-1], p[0]));
}
return ans * 0.5;
}

void halfPlane()
{
finalNumber = 0;
sort(line, line + lineNumber);
stackTop = 1;
int bottom = 0, tmp = 1;
for(int i=1;i<lineNumber;++i)
{
if(line[i].ang - line[i-1].ang > eps)
{
line[tmp++] = line[i];
}
}
lineNumber = tmp;
stack[0] = line[0], stack[1] = line[1];
for(int i=2;i<lineNumber;++i)
{
if(isParellel(stack[stackTop], stack[stackTop-1]) || isParellel(stack[bottom], stack[bottom+1]))
{
return;
}
while(bottom < stackTop && judgeOut(line[i], stack[stackTop] * stack[stackTop - 1]))
{
-- stackTop;
}
while(bottom < stackTop && judgeOut(line[i], stack[bottom] * stack[bottom + 1]))
{
++ bottom;
}
stack[++stackTop] = line[i];
}
while(bottom < stackTop && judgeOut(stack[bottom], stack[stackTop] * stack[stackTop - 1]))
{
-- stackTop;
}
while(bottom < stackTop && judgeOut(stack[stackTop], stack[bottom] * stack[bottom + 1]))
{
++ bottom;
}
if(stackTop <= bottom + 1)
{
return;
}
stack[++stackTop] = stack[bottom];
for(int i=bottom;i<stackTop;++i)
{
final[finalNumber ++] = stack[i] * stack[i+1];
}
}

int main()
{
while(~scanf("%d", &n))
{
for(int i=0;i<n;++i)
{
machine[i].input3d();
}
for(int i=0;i<3;++i)
{
light[i].input2d();
}
lineNumber = 0;
for(int i=0;i<3;++i)
{
shadowNumber[i] = 0;
for(int j=0;j<n;++j)
{
shadow[i][j].x = 100.0 * (light[i].x - machine[j].x) / (machine[j].z - 100.0) + light[i].x;
shadow[i][j].y = 100.0 * (light[i].y - machine[j].y) / (machine[j].z - 100.0) + light[i].y;
}
graham(shadow[i], n, shadowPoly[i], shadowPolyTop[i]);
shadowPoly[i][shadowPolyTop[i]] = shadowPoly[i][0];
for(int j=0;j<shadowPolyTop[i];++j)
{
line[lineNumber].a = shadow[i][shadowPoly[i][j]];
line[lineNumber].b = shadow[i][shadowPoly[i][j+1]];
line[lineNumber].ang = atan2(line[lineNumber].b.y - line[lineNumber].a.y,
line[lineNumber].b.x - line[lineNumber].a.x);
++ lineNumber;
}
}
halfPlane();
printf("%.2lf\n", getArea(final, finalNumber));
}
return 0;
}