UVA 12303 Composite Transformations [三维仿射矩阵][更正:旋转轴归一化和平面三点式问题]

三维仿射变换什么的,纠结了很多天想把绕任意轴旋转的齐次矩阵手推出来,最后,结果连sample都过不了,只能等回头再调了,目前先用着模板。

所谓仿射变换,咱们通俗的讲就是通过补维把坐标转化为齐次矩阵的形式进行矩阵乘法,这样就能找到平移、拉伸、旋转的通式解。
此处要用到的有三个变换矩阵:
平移:(偏移量(dx,dy,dz))
1 0 0 dx
0 1 0 dy
0 0 1 dz
0 0 0 1
放缩:(比例系数(sx,sy,sz))
sx 0 0 0
0 sy 0 0
0 0 sz 0
0 0 0 1
旋转:(绕过原点的向量(x,y,z)转d弧度角)(注意旋转轴要归一化。。)
(1-cos(d))*x*x+cos(d)     (1-cos(d))*x*y-sin(d)*z   (1-cos(d))*x*z+sin(d)*y   0
(1-cos(d))*y*x+sin(d)*z   (1-cos(d))*y*y+cos(d)     (1-cos(d))*y*z-sin(d)*x   0
(1-cos(d))*z*x-sin(d)*y   (1-cos(d))*z*y+sin(d)*x   (1-cos(d))*z*z+cos(d)     0
          0                             0                                          0                                   1

而变换面的话,就从三维平面的数学表示式入手。
面方程ax+by+cz=d->法向量dir(a,b,c),必在其上的点on(1,1,(z==0)?1:z) 得到第二个点out=on+dir
通过旋转on和out这两个点来重构平面->得到新的dir和新的on,通过dir和on就可得到新的面方程。

上述点法式在缩放矩阵乘的时候会出现问题。因为平面外的点不能受统一的缩放变换。(调了一晚上才想明白,我擦嘞。。。)
所以果断要用三点式确定平面。根据平面方程if else构造出平面上不共线的三点,然后旋转这三个点,然后通过三维叉积算出新的法向量,再利用点法式构造出新平面。
唉,偷懒总会死的。。竟然帮乳鸽查出了他spj的错误。。之前的代码在平面部分没有处理好,而且绕轴旋转矩阵的旋转轴没有归一化,这样竟然都ac了。
随手出个样例都能挂掉之前的代码,而且平面部分算是完全写扯了。。。


而本题的特殊之处在于,50000个点50000个平面,1000次变换。如果裸做的话复杂度5*10^7*大常数。所以TLE是必然的。
那么由于矩阵连乘的特殊性,我们就可以按顺序把变换矩阵连乘起来,然后让点和面这100000个列向量右乘最终变换矩阵,就可得到答案。最终复杂度为50000*4^2左右。

话说用了cos和sin精度必然丢的很厉害。亏着这题比较厚道,才需要两位精度。

计算几何神马的,太有爱了~

下附代码:
#include <iostream>
#include 
<cassert>
#include 
<cmath>
#include 
<cstring>
#include 
<cstdio>
using namespace std;
const double eps = 1e-8;
const double pi = acos(-1.0);
int sgn(const double x){return fabs(x) < eps?0:x < -eps?-1:1;}
#define sqr(x) ((x) * (x))
#define REP(i,x) for(int i = 0;i < x;i ++)
#define rep(i,j,k) for(int (i) = (j);(i) < (k);(i) ++ )
#define maxn 50005
#define maxm maxn
struct matrix
{
    
int r,c;
    
double v[10][10];
    matrix(){}
    matrix(
int a,int b):r(a),c(b)
    {
        REP(i,r) 
            fill(v[i],v[i]
+c,0.0);
        
if(r==c)
        {
            REP(i,r)
//unit matrix as default matrix
                v[i][i]=1;
        }
    }
    
void init()
    {
        REP(i,r) 
            fill(v[i],v[i]
+c,0.0);
        
if(r==c)
        {
            REP(i,r)
//unit matrix as default matrix
                v[i][i]=1;
        }
    }
    matrix 
operator * (const matrix p)
    {
        matrix ans(r,p.c);
        REP(i,r)
            REP(j,ans.c)
            {
                
double tmp = 0.0;
                REP(k,c)
                    tmp 
+= v[i][k] * p.v[k][j];
                ans.v[i][j] 
= tmp;
            }
        
return ans;
    }
    
void print()
    {
        rep(i,
0,r)
            rep(j,
0,c)
            printf(
"%.2lf%c",v[i][j],j==c-1?'\n':' ');
    }
};
struct Tpoint
{
    
double x,y,z;
    Tpoint(){}
    Tpoint(
double a,double b,double c):x(a),y(b),z(c){}
    Tpoint 
operator - (const Tpoint p){return Tpoint(x-p.x,y-p.y,z-p.z);}
    Tpoint 
operator + (const Tpoint p){return Tpoint(x+p.x,y+p.y,z+p.z);}
    
double operator * (const Tpoint p){return x*p.x+y*p.y+z*p.z;}
    Tpoint 
operator ^ (const Tpoint p)
    {
return Tpoint(y * p.z - z * p.y,z * p.x - x * p.z,x * p.y - y * p.x);}
    
bool operator == (const Tpoint p){return sgn(x-p.x)==0&&sgn(y-p.y)==0&&sgn(z-p.z)==0;}
    
double norm2(){return (*this)*(*this);}
    
double norm(){return sqrt(norm2());}
    Tpoint unit(){
return Tpoint(x/norm(),y/norm(),z/norm());}
    Tpoint gao(matrix h)
    {
        matrix o(
4,1);
        o.v[
0][0= x;o.v[1][0= y;o.v[2][0= z;o.v[3][0= 1;
        o 
= h * o;
        
return Tpoint(o.v[0][0],o.v[1][0],o.v[2][0]);
    }
    
void scan(){scanf("%lf%lf%lf",&x,&y,&z);}
    
void print(){printf("%.2lf %.2lf %.2lf\n",x,y,z);}
}pt[maxn];
Tpoint tZero(
0,0,0);
struct Tplane
{
    
double a,b,c,d;
    
//a + b + c = d
    Tplane(){}
    Tplane(
double e,double f,double g,double h):a(e),b(f),c(g),d(-h){}
    
//return a perpendicular vector of the given plane
    Tpoint ndir(){return Tpoint(a,b,c);}
    
//shifting transform
    double norm2(){return sqr(a)+sqr(b)+sqr(c);}
    
double norm(){return sqrt(norm2());}
    Tplane unit(){
return Tplane(a/norm(),b/norm(),c/norm(),-d/norm());}
    
void print()
    {printf(
"%.2lf %.2lf %.2lf %.2lf\n",a,b,c,-d);}
}tp[maxm];
int n,m,t;
char opt[15];
int main()
{
    
double a,b,c,d;
    scanf(
"%d %d %d",&n,&m,&t);
    REP(i,n)
    {
        scanf(
"%lf %lf %lf",&a,&b,&c);
        pt[i] 
= Tpoint(a,b,c);
    }
    REP(i,m)
    {
        scanf(
"%lf %lf %lf %lf",&a,&b,&c,&d);
        tp[i] 
= Tplane(a,b,c,d);
    }
    matrix handle(
4,4),tran(4,4);
    REP(i,t)
    {
        scanf(
"%s",opt);
        
if(opt[0== 'T')
        {
            scanf(
"%lf %lf %lf",&a,&b,&c);
            tran.init();
            tran.v[
0][3= a;
            tran.v[
1][3= b;
            tran.v[
2][3= c;
            handle 
= tran * handle;
        }
        
else if(opt[0== 'R')
        {
            scanf(
"%lf %lf %lf %lf",&a,&b,&c,&d);
            tran.init();
            
double t = d / 180.0 * pi;
            
double ct = cos(t),st = sin(t);
            
double fuck = sqrt(a * a + b * b + c * c);
            a 
/= fuck;
            b 
/= fuck;
            c 
/= fuck;
            tran.v[
0][0= (1.0 - ct) * a * a + ct;
            tran.v[
0][1= (1.0 - ct) * a * b - st * c;
            tran.v[
0][2= (1.0 - ct) * a * c + st * b;
            tran.v[
1][0= (1.0 - ct) * b * a + st * c;
            tran.v[
1][1= (1.0 - ct) * b * b + ct;
            tran.v[
1][2= (1.0 - ct) * b * c - st * a;
            tran.v[
2][0= (1.0 - ct) * c * a - st * b;
            tran.v[
2][1= (1.0 - ct) * c * b + st * a;
            tran.v[
2][2= (1.0 - ct) * c * c + ct;
            handle 
= tran * handle;
        }
        
else if(opt[0== 'S')
        {
            scanf(
"%lf %lf %lf",&a,&b,&c);
            tran.init();
            tran.v[
0][0= a;
            tran.v[
1][1= b;
            tran.v[
2][2= c;
            handle 
= tran * handle;
        }
    }
    REP(i,n)
    {
        pt[i] 
= pt[i].gao(handle);
        pt[i].print();
    }
    REP(i,m)
    {
        a 
= tp[i].a,b = tp[i].b,c = tp[i].c,d = tp[i].d;
        
double limit = max(fabs(a),max(fabs(b),fabs(c)));
        Tpoint on[
3];
        
if(fabs(a) == limit)
        {
            on[
0= Tpoint(d/a,0,0);
            on[
1= Tpoint((d-b)/a,1,0);
            on[
2= Tpoint((d-c)/a,0,1);
        }
        
else if(fabs(b) == limit)
        {
            on[
0= Tpoint(0,d/b,0);
            on[
1= Tpoint(1,(d-a)/b,0);
            on[
2= Tpoint(0,(d-c)/b,1);
        }
        
else
        {
            on[
0= Tpoint(0,0,d/c);
            on[
1= Tpoint(1,0,(d-a)/c);
            on[
2= Tpoint(0,1,(d-b)/c);
        }
        rep(j,
0,3)
            on[j] 
= on[j].gao(handle);
        Tpoint dir 
= (on[2- on[0]) ^ (on[1- on[0]);
        tp[i] 
= Tplane(dir.x,dir.y,dir.z,-(dir * on[0]));
        tp[i].unit().print();
    }
}

posted on 2011-11-04 01:40 BUPT-[aswmtjdsj] @ Penalty 阅读(707) 评论(0)  编辑 收藏 引用 所属分类: 计算几何UVA Solution Report


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


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

导航

统计

常用链接

留言簿(1)

随笔分类(150)

随笔档案(71)

搜索

积分与排名

最新评论

阅读排行榜

评论排行榜