HDU 4087 三维坐标旋转(仿射变换) 矩阵加速 && 2011 Asia Beijing Regional Contest

这是第一次做仿射变换的题目,搞了一下午。

题意：给你一系列对点的平移，缩放，绕任意轴旋转的操作，然后给你点要求对所有点都进行这些操作，输出操作后所有点的坐标。

如果把每个点对每个操作依次进行变换会超时，因为有重复操作最大是2^32次操作。可以把对坐标的操作转换成矩阵乘以矩阵，这个操作是仿射变换。

例如：

1 0 0 tx    x      x+tx

0 1 0 ty * y  =  y+ty

0 0 1 tz    z      z+tz

0 0 0 1    1     1

缩放与旋转也是同理 不过旋转是抄袭网上某大神导出的公式

以下是仿射变换矩阵：

平移

translate tx ty tz

1 0 0 tx

0 1 0 ty

0 0 1 tz

0 0 0 1

缩放

scale a b c

a  0  0  0

0  b  0  0

0  0  c  0

0  0  0  1

绕任意轴(过原点)旋转(注意要把轴向量归一化，否则点在旋转轴上时有问题)

rotate 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

(ps:如果是不过原点的情况，可以把所轴线平移到过原点的情况，并把所有的点平移过来，旋转后再平移回去，可能会有精度误差。)

附上代码:

#include<stdio.h>
#include<algorithm>
#include<string.h>
#include<math.h>
#include<stdlib.h>
#include<iostream>
#include<stack>
#include<vector>
using namespace std;
const int maxn=500000;
const double eps=1e-6;
const double pi = acos(-1.0);

struct node	//4*4矩阵
{
double a[4][4];

void init()	//初始化
{
memset(a,0,sizeof a);
a[0][0] = 1;
a[1][1] = 1;
a[2][2] = 1;
a[3][3] = 1;
}

node()
{
memset(a,0,sizeof a);
a[0][0] = 1;
a[1][1] = 1;
a[2][2] = 1;
a[3][3] = 1;
}

node(double b[][4])
{
for(int i = 0 ; i < 4; ++ i)
{
for(int j = 0 ; j < 4 ; ++ j)
{
a[i][j] = b[i][j];
}
}
}

node operator * (const node& n) //矩阵乘法
{
node note(a);
for(int i = 0 ; i < 4 ; ++ i)
{
for(int j = 0 ; j < 4 ; ++ j)
{
double tmp = 0;
for(int k = 0 ; k < 4 ; ++ k)
{
tmp += a[i][k]*n.a[k][j];
}
note.a[i][j] = tmp;
}
}
return note;
}
}n;

void translate(node& n, double x, double y, double z)	//平移
{
node tmp;
tmp.a[0][3] = x;
tmp.a[1][3] = y;
tmp.a[2][3] = z;
n = tmp * n;
}

void scale(node& n, double a, double b, double c)	//缩放
{
node tmp;
tmp.a[0][0] = a;
tmp.a[1][1] = b;
tmp.a[2][2] = c;
n = tmp*n;
}

void rotate(node& n, double x, double y, double z, double d)	//旋转
{
d = d/180*pi;
double di = sqrt(x*x+y*y+z*z);	//单位化
x/=di, y/=di, z/=di;
node tmp;
tmp.a[0][0] = (1-cos(d))*x*x + cos(d);
tmp.a[0][1] = (1-cos(d))*x*y - z*sin(d);
tmp.a[0][2] = (1-cos(d))*x*z + y*sin(d);
tmp.a[1][0] = (1-cos(d))*x*y + z*sin(d);
tmp.a[1][1] = (1-cos(d))*y*y + cos(d);
tmp.a[1][2] = (1-cos(d))*y*z - x*sin(d);
tmp.a[2][0] = (1-cos(d))*x*z - y*sin(d);
tmp.a[2][1] = (1-cos(d))*y*z + x*sin(d);
tmp.a[2][2] = (1-cos(d))*z*z + cos(d);
n = tmp*n;
}

void repeat(node& n, int b)	//矩阵加速(快速幂思维)
{
node tmp;
while(b)
{
if(b&1)
{
tmp = n*tmp;
}
n = n*n;
b/=2;
}
n = tmp;
}

node rep(int n)	//递归思维
{
node tmp; tmp.init();
char s[100];
while(~scanf("%s", s))
{
if(s[0]=='t')
{
double x, y, z;
scanf("%lf%lf%lf",&x,&y,&z);
translate(tmp,x,y,z);
continue;
}
if(s[0]=='s')
{
double x, y, z;
scanf("%lf%lf%lf",&x,&y,&z);
scale(tmp,x,y,z);
continue;
}
if(s[0]=='r'&&s[1]=='o' )
{
double x, y, z, d;
scanf("%lf%lf%lf%lf",&x,&y,&z,&d);
rotate(tmp,x,y,z,d);
continue;
}
if(s[0]=='r'&&s[1]=='e' )
{
int n;
scanf("%d", &n);
tmp = rep(n)*tmp;
}
if(s[0]=='e')
{
repeat(tmp,n);
return tmp;
}
}
}

struct Point //点
{
double a[4];
void read() {scanf("%lf%lf%lf", a, a+1, a+2);a[3]=1;}

void operator*= (const node& n)  //矩阵*点
{
double tx = a[0], ty = a[1], tz =a[2], td = a[3];
for(int i = 0 ; i < 4 ; ++ i)
{
a[i] = n.a[i][0]*tx + n.a[i][1]*ty + n.a[i][2]*tz + n.a[i][3]*td;
}
}
};
int main()
{
//freopen("in.txt","r",stdin);
int n;
while(scanf("%d", &n)&&n)
{
char s[100];
node tmp;
tmp = rep(1)*tmp;
Point p;
for(int i = 0 ; i < n  ; ++ i)
{
p.read();
p*=tmp;
printf("%.2lf %.2lf %.2lf\n",p.a[0]+eps,p.a[1]+eps,p.a[2]+eps);//小心 -0.00 !
}
printf("\n");

}
return 0;
}