好恶心的一题。
给一个点，一个多面体，一个平面，假设点是光源，那么就问，多面体在平面上的投影面积为多少。
做法：
hint还是比较人道的，提示了旋转平面把三维凸包转成二维凸包的方法。（其实不能算是三维凸包，只不过是三维空间内的二维凸包；hint转化之后就变成了二维空间的二维凸包。）
平面给的是平面方程Ax+By+Cz=D，那么平面法向量为(A,B,C)。目标是把平面旋转到和Z平面平行。根据hint的提示，先把(A,B,C)->(0,sqrt(A^2+B^2),C);再把(0,sqrt(A^2+B^2),C)->(0,0,sqrt(A^2+B^2+C^2))。
这样就把三维仿射变换拆解成了两个二维仿射变换。
根据二维变换推导，两个三维仿射变换矩阵分别为
|cos(a),sin(a),0|
ROTxy=|-sin(a),cos(a),0|
|0,0,1|
|1,0,0|
ROTyz=|0,cos(b),sin(b)|
|0,-sin(b),cos(b)|
然后根据对应关系求出a、b分别是多少。
Acosa-Bsina=0 ,a = atan2(A,B)
sqrt(A^2+B^2)cosb-Csinb=0,b = atan2(sqrt(A^2+B^2),C)
从而解出两个仿射矩阵的旋转角，回代，然后把两个矩阵相乘得到ROTxyz。
再在原平面上选一个基准点，用于构造新平面，我选的是(1,1,?)或者(1,?,1)或者(?,1,1)，此处？根据平面方程算出，并由A、B、C是否为0确定？在哪个维度；接下来将基准点带入新平面方程即可得出新平面的D'，这样新平面就构造出来了。
之后就是用基准点和原图中的多面体和光源点构造向量进行旋转，再根据基准点和新向量构造出这些点的新坐标。
再之后就是由光源点和多面体上的点连线向平面上做投影了，此处用直线、平面相交函数求出（基本空见解析知识，一推即得），同时判断是否有fly现象存在，fly现象即，某条投影线和平面平行，或者光源点和某点的连线向与平面相反的方向发出（多面体的点不被光源点和投影点所夹）。如果所有点皆fly，那么投影面积为0；如果既有fly又有投影点，那么Infi；如果没有fly现象，那么就可以下一步了。
下一步就是平面凸包+叉积算有向面积。。我这做的稍微麻烦了点，还重新构造了二维点什么的。。。略去不表。反正怎么搞都行了。
写了半个晚上，真够难写的。三维几何什么的伤不起。现在三维仿射变换理解的还是有些问题，绕轴旋转什么的还是不太会写。
附巨矬代码：
#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;
#define maxn 105
const double eps = 1e-8;
int comp(double x)
{
return (fabs(x) < eps?0:x < -eps ? -1 : 1);
}
struct point
{
double x,y,z;
point(){}
point(double a,double b,double c):x(a),y(b),z(c){}
point operator -(const point a)
{
return point(x - a.x,y - a.y,z - a.z);
}
point operator +(const point p)
{
return point(x + p.x,y + p.y,z + p.z);
}
double operator *(const point p)
{
return x * p.x + y * p.y + z * p.z;
}
point operator ^(const point p)
{
return point(y * p.z - p.y * z,p.x * z - x * p.z,x * p.y - p.x * y);
}
double len()
{
return sqrt(x * x + y * y + z * z);
}
}poly[maxn],newpoly[maxn],newcenter,center;
point pvec(point p1,point p2,point p3)
{
return (p1 - p2) ^ (p1 - p3);
}
point rotate(point vec,double mat[][4])
{
point temp;
temp.x = vec.x * mat[0][0] + vec.y * mat[1][0] + vec.z * mat[2][0];
temp.y = vec.x * mat[0][1] + vec.y * mat[1][1] + vec.z * mat[2][1];
temp.z = vec.x * mat[0][2] + vec.y * mat[1][2] + vec.z * mat[2][2];
return temp;
}
double A,B,C,D;
int n;
point inter(point p1,point p2,point p3,point u,point v)
{
point ret = pvec(p1,p2,p3);
double t = ret * (p1 - u) / (ret * (v - u));
point vec = v - u;
ret.x = u.x + vec.x * t;
ret.y = u.y + vec.y * t;
ret.z = u.z + vec.z * t;
return ret;
}
bool between(point a,point b,point c)
{
return ((comp(a.x - b.x) <= 0 && comp(b.x - c.x) <= 0) || (comp(c.x - b.x) <= 0 && comp(b.x - a.x) <= 0)) && ((comp(a.y - b.y) <= 0 && comp(b.y - c.y) <= 0) || (comp(c.y - b.y) <= 0 && comp(b.y - a.y) <= 0)) && ((comp(a.z - b.z) <= 0 && comp(b.z - c.z) <= 0) || (comp(c.z - b.z) <= 0 && comp(b.z - a.z) <= 0)) ;
}
struct Dpoint
{
double x,y;
Dpoint(){}
Dpoint(double a,double b):x(a),y(b){}
Dpoint operator -(const Dpoint p)
{
return Dpoint(x - p.x,y - p.y);
}
double operator *(const Dpoint p)
{
return x * p.x + y * p.y;
}
double norm()
{
return sqrt(x * x + y * y);
}
double operator ^(const Dpoint p)
{
return x * p.y - y * p.x;
}
}p[maxn],stack[maxn],c;
bool cmp1(const Dpoint &a,const Dpoint &b)
{
return comp(a.x - b.x) < 0 || (comp(a.x - b.x) == 0 && comp(a.y - b.y) < 0);
}
bool cmp2(Dpoint a,Dpoint b)
{
double phi1 = atan2(a.y - c.y,a.x - c.x);
double phi2 = atan2(b.y - c.y,b.x - c.x);
return comp(phi1 - phi2) < 0 || (comp(phi1 - phi2) == 0 && comp(a.norm() - b.norm()) < 0);
}
double cha(Dpoint a,Dpoint b,Dpoint o)
{
return (b - o) ^ (a - o);
}
int now;
void hull()
{
stack[0] = p[0];
now = 0;
for(int i = 1;i < n;i++)
{
while(now >= 2 && comp(cha(p[i],stack[now],stack[now-1])) <= 0)
now--;
now++;
stack[now] = p[i];
}
}
int main()
{
while(scanf("%lf %lf %lf %lf",&A,&B,&C,&D) == 4)
{
if(!comp(A) && !comp(B) && !comp(C) && !comp(D))
break;
scanf("%d",&n);
for(int i = 0;i < n;i++)
{
double x,y,z;
scanf("%lf %lf %lf",&x,&y,&z);
poly[i] = point(x,y,z);
}
scanf("%lf %lf %lf",¢er.x,¢er.y,¢er.z);
//a * x + b * y + c * z = d
//vertical vector (a,b,c)
//rotate it to (0,0,sqrt(a^2+b^2+c^2))
point vertical(A,B,C);
point mid(0,sqrt(A*A+B*B),C);
point final(0,0,sqrt(A*A+B*B+C*C));
double theta,phi;
phi = atan2(A,B);
theta = atan2(sqrt(A*A+B*B),C);
double ROT1[][4] = {cos(phi),sin(phi),0,0,
-sin(phi),cos(phi),0,0,
0,0,1,0};
double ROT2[][4] = {1,0,0,0,
0,cos(theta),sin(theta),0,
0,-sin(theta),cos(theta),0};
double ROT3[4][4];
for(int i = 0;i < 3;i++)
for(int j = 0;j < 3;j++)
{
double aaa = 0.0;
for(int k = 0;k < 3;k++)
aaa += ROT1[i][k] * ROT2[k][j];
ROT3[i][j] = aaa;
}
point whatelse = rotate(vertical,ROT3);
point on;
//find a point on the original plane
if(comp(A) != 0)
on = point((D - C - B) / A,1,1);
else if(comp(B) != 0)
on = point(1,(D - C - A) / B,1);
else if(comp(C) != 0)
on = point(1,1,(D - B - A) / C);
double newd = final.x * on.x + final.y * on.y + final.z * on.z;
point tempto,to;
for(int i = 0;i < n;i++)
{
tempto = poly[i] - on;
to = rotate(tempto,ROT3);
newpoly[i] = on + to;
}
tempto = center - on;
to = rotate(tempto,ROT3);
newcenter = on + to;
bool onplane = false,fly = false;
int con = 0;
point fuck1(0,0,on.z),fuck2(1,1,on.z),fuck3(1,0,on.z);
for(int i = 0;i < n;i++)
{
point she = newpoly[i] - newcenter;
if(comp(she * final) == 0)
{
fly = true;
continue;
}
point cross = inter(fuck1,fuck2,fuck3,newcenter,newpoly[i]);
if(between(cross,newpoly[i],newcenter))
{
p[con++] = Dpoint(cross.x,cross.y);
onplane = true;
}
else
fly = true;
}
if(fly && onplane)
puts("Infi");
else if((fly && !onplane) || (n == 0 || n == 1 || n == 2))
puts("0.00");
else
{
int mark = 0;
now = 0;
for(int i = 0;i < con;i++)
{
if(i > 0)
{
if(cmp1(p[i],c))
{
c = p[i];
mark = i;
}
}
else
c = p[0];
}
p[mark] = p[0];
p[0] = c;
sort(p + 1,p + con,cmp2);
hull();
double ans = 0;
for(int i = 1;i < now;i++)
ans += fabs(cha(stack[i],stack[i+1],stack[0])/2.0);
printf("%.2lf\n",ans);
}
}
}
Down