UVA 5100 && FZU 2002 Shade of Hallelujah Mountain 【三维旋转、投影+凸包】

好恶心的一题。
给一个点,一个多面体,一个平面,假设点是光源,那么就问,多面体在平面上的投影面积为多少。

做法:
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",&center.x,&center.y,&center.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);
        }
    }
}

posted on 2011-09-07 10:03 BUPT-[aswmtjdsj] @ Penalty 阅读(339) 评论(1)  编辑 收藏 引用 所属分类: 计算几何UVAlive Solution ReportFZU Solution Report

评论

# re: UVA 5100 && FZU 2002 Shade of Hallelujah Mountain 【三维旋转、投影+凸包】 2011-10-01 13:55 BUPT-[aswmtjdsj] @ Penalty

唉。。话说这题的精度损失有点严重。。同一份代码在HDU上面过不去啊唉。。  回复  更多评论   


只有注册用户登录后才能发表评论。
【推荐】超50万行VC++源码: 大型组态工控、电力仿真CAD与GIS源码库
网站导航: 博客园   IT新闻   BlogJava   知识库   博问   管理


<2011年9月>
28293031123
45678910
11121314151617
18192021222324
2526272829301
2345678

导航

统计

常用链接

留言簿(1)

随笔分类(150)

随笔档案(71)

搜索

积分与排名

最新评论

阅读排行榜

评论排行榜