http://blog.csdn.net/popy007/article/details/376952
<6>3-D空间中的基变换与坐标变换
-Twinsen编写
-本人水平有限,疏忽错误在所难免,还请各位数学高手、编程高手不吝赐教
-我的Email-address: popyy@netease.com
一、空间坐标系的基和基矩阵
在3-D空间中,我们用空间坐标系来规范物体的位置,空间坐标系由3个相互垂直的坐标轴组成,我们就把它们作为我们观察3-D空间的基础,空间中物体的位置可以通过它们来衡量。当我们把这3个坐标轴上单位长度的向量记为3个相互正交的单位向量i,j,k,空间中每一个点的位置都可以被这3个向量线性表出,如P<1,-2,3>这个点可以表为i-2j+3k。
我们把这3个正交的单位向量称为空间坐标系的基,它们单位长度为1且正交,所以可以成为标准正交基。三个向量叫做基向量。现在我们用矩阵形式写出基向量和基。
i = | 1 0 0 |
j = | 0 1 0 |
k = | 0 0 1 |
| i | | 1 0 0 |
B = | j | = | 0 1 0 |
| k | | 0 0 1 |
这样的矩阵我们叫它基矩阵。有了基矩阵,我们就可以把空间坐标系中的一个向量写成坐标乘上基矩阵的形式,比如上面的向量P可以写成:
P = C x B
=>
| 1 0 0 |
| 1 -2 3 | = | 1 -2 3 | x | 0 1 0 |
| 0 0 1 |
这样的话,空间坐标系下的同一个向量在不同的基下的坐标是不同的。
二、局部坐标系和局部坐标
和空间坐标系(也可以叫做全局坐标系或者世界坐标系)并存的称为局部坐标系(也叫坐标架——coordinate frame),它有自己的基,这些基向量把空间坐标系作为参考系。比如
| x'| | -1 0 0 |
B' = | y'| = | 0 1 0 |
| z'| | 0 0 -1 |
| x''| | 2^½ /2 0 2^½ /2 |
B'' = | y''| = | 0 -1 0 |
| z''| | -(2^½) /2 0 2^½ /2 |
就是两个局部坐标系的基,如图:
现在我们可以把上面那个空间坐标中的向量P|1 -2 3|(以后都用矩阵表示)表示在不同的基下,我把它写成一个大长串的式子:
| x' | | x''|
P = | Px' Py' Pz' | x | y' | = | Px'' Py'' Pz'' | x | y''|
| z' | | z''|
这里| Px' Py' Pz'|是P在B'下的坐标,| Px'' Py'' Pz''|是P在B''下的坐标,我把它写的具体点吧:
| -1 0 0 | | 2^½ /2 0 2^½ /2|
| 1 -2 3 | = | -1 -2 -3 | x | 0 1 0 | = | 2*2^½ -2 2^½ | x | 0 -1 0 |
| 0 0 -1 | | -(2^½) /2 0 2^½ /2|
这就是说,在空间坐标系下面的向量| 1 -2 3 |在基B'下的坐标为|-1 -2 -3|,在B''下的坐标为| 2*2^½ -2 2^½ |。当然空间坐标系也有自己的基B|i j k|^T(因为是列向量,所以写成行向量的转置),但我们现在是拿它当作一个参考系。
在研究了局部坐标系之后,我现在要分析两个应用它们的例子,先来看
三、空间坐标系中一个点围绕任一轴的旋转
上一篇讨论3-D空间旋转的时候说到有一个高档的方法做3-D空间任意轴旋转,现在我们的知识储备已经足够理解这个方法了(Quake引擎使用的就是这个方法)。
如上所示,空间坐标系中的一个局部坐标系xyz中有一个向量a(2,5,3)和一个点p(8,4,2)现在我要让p点围绕a向量旋转60度,得到p’点,该如何做呢?从目前掌握的旋转知识来看,我们有两个理论基础:
1)在一个坐标系中的一个点,如果要它围绕该坐标系中一个坐标轴旋转,就给它的坐标值乘相应的旋转矩阵,如
[cosA -sinA 0 ]
[sinA cosA 0 ]
[0 0 1 ]
等等。
2)我们已经学习了局部坐标系的理论了,知道空间中一个点在不同的坐标系中的坐标不同。利用这一点,我们可以很方便的让一个点或者向量在不同的坐标系之间转换。
我们联系这两个理论根据,得出我们的思路:
1构造另一个局部坐标系abc,使得a成为该坐标系的一个坐标轴。
2 把p的坐标变换到abc中,得到p’,用旋转公式让p’围绕已经成为坐标轴的a旋转,得到p’’。
3把p’’再变换回坐标系xyz,得到p’’’,则p’’’就是p围绕a旋转后的点。
下面我们逐步说明。
首先我们构造abc,我们有无数种方法构造,因为只要保证b、c之间以及他们和a之间都正交就可以了,但我们只要一个。根据上图,我们首先产生一个和a正交的b。这可以通过向量的叉乘来完成:我们取另一个向量v(显然,这个向量是不能和a共线的任何非零向量),让它和a决定一个平面x,然后让v叉乘a得到一个垂直于x的向量b,因为b垂直于x,而a在平面x上,因此b一定垂直于a,然后用a叉乘b得到c,最后单位化a、b、c,这样就得到了局部坐标系abc。
然后我们把p点变换到abc坐标系中,得到p’,即p’就是p在abc中的坐标:
|a b c| * p’= |x y z| * p
p’ = |a b c|^-1 * |x y z| * p
|ax bx cx| |1 0 0| |px|
p’ = |ay by cy| ^-1 * |0 1 0| * |py|
|az bz cz| |0 0 1| |pz|
注意这里|a b c|^-1即矩阵|a b c|的逆矩阵,因为a、b、c是三个正交向量,并且是单位向量,因此|a b c|是一个正交矩阵,正交矩阵的转置和逆相等,这是它的一个特性,因此上面的公式就可以写成:
|ax ay az| |1 0 0| |px|
p’ = |bx by bz| * |0 1 0| * |py|
|cx cy cz| |0 0 1| |pz|
这个时候p’就是p在abc坐标系下的坐标了。此时a已经是一个坐标轴了,我们可以用旋转矩阵来做。
p’’ = RotMatrix * p’
[1 0 0] |p’x|
p’’ = [0 cos60 -sin60] * |p’y|
[0 sin60 cos60] |p’z|
最后,我们把p’’再次变换回xyz坐标系,得到最终的p’’’
|a b c| * p’’ = |x y z| * p’’’
p’’’ = |x y z|^-1 * |a b c| * p’’
p’’’ = |a b c| * p’’
最后
p’’’ = |a b c| * RotMatrix * |a b c|^T * p = M * p
这样就得到了xyz坐标系中点p围绕a旋转60度后的点。
最后,我用Quake3引擎的相应函数(来自idSoftware ——quake3-1[1].32b-source——mathlib.c)来完成对这个算法的说明:
/*
===============
RotatePointAroundVector
dst是一个float[3],也就是p’’’
dir相当于a,point就是p,degrees是旋转度数
===============
*/
void RotatePointAroundVector( vec3_t dst, const vec3_t dir, const vec3_t point,
float degrees ) {
float m[3][3];
float im[3][3];
float zrot[3][3];
float tmpmat[3][3];
float rot[3][3];
int i;
vec3_t vr, vup, vf;
float rad;
vf[0] = dir[0];
vf[1] = dir[1];
vf[2] = dir[2];
// 首先通过dir得到一个和它垂直的vr
// PerpendicularVector()函数用于构造和dir垂直的向量
// 也就是我们上面的第1步
PerpendicularVector( vr, dir );
// 通过cross multiply得到vup
// 现在已经构造出坐标轴向量vr, vup, vf
CrossProduct( vr, vf, vup );
// 把这三个单位向量放入矩阵中
m[0][0] = vr[0];
m[1][0] = vr[1];
m[2][0] = vr[2];
m[0][1] = vup[0];
m[1][1] = vup[1];
m[2][1] = vup[2];
m[0][2] = vf[0];
m[1][2] = vf[1];
m[2][2] = vf[2];
// 产生转置矩阵im
memcpy( im, m, sizeof( im ) );
im[0][1] = m[1][0];
im[0][2] = m[2][0];
im[1][0] = m[0][1];
im[1][2] = m[2][1];
im[2][0] = m[0][2];
im[2][1] = m[1][2];
// 构造旋转矩阵zrot
memset( zrot, 0, sizeof( zrot ) );
zrot[0][0] = zrot[1][1] = zrot[2][2] = 1.0F;
rad = DEG2RAD( degrees );
zrot[0][0] = cos( rad );
zrot[0][1] = sin( rad );
zrot[1][0] = -sin( rad );
zrot[1][1] = cos( rad );
// 开始构造变换矩阵M
// tmpmat = m * zrot
MatrixMultiply( m, zrot, tmpmat );
// rot = m * zrot * im
MatrixMultiply( tmpmat, im, rot );
// 则 rot = m * zrot * im 和我们上面推出的
// M = |a b c| * RotMatrix * |a b c|^T 一致
// 变换point这个点
// p’’’ = M * p
for ( i = 0; i < 3; i++ ) {
dst[i] = rot[i][0] * point[0] + rot[i][1] * point[1] + rot[i][2] * point[2];
}
}
四、世界空间到相机空间的变换
空间坐标系XYZ,相机坐标系UVN。这时候相机空间的基(以下简称相机)在空间坐标系中围绕各个坐标轴旋转了一定角度<a,b,c>,然后移动了<x,y,z>。对于模型我们可以看作相对于相机的逆运动,即模型旋转了一定角度<-a,-b,-c>,然后移动了<-x,-y,-z>,可以把相机和物体的运动看成两个互逆的变换。这样,可以通过对相机的变换矩阵求逆来得到模型的变换矩阵。下面来具体看一下,如何得到相机变换矩阵,并且求得它的逆矩阵。
首先声明一下,对于一个模型的变换,我们可以给模型矩阵左乘变换矩阵:
M x P = P'
| A B C D | | x | | Ax + By + Cz + D |
| E F G H | | y | | Ex + Fy + Gz + H |
x =
| I J K L | | z | | Ix + Jy + Kz + L |
| M N O P | | 1 | | Mx + Ny + Oz + P |
也可以右乘变换矩阵:
P^T x M^T = P'^T
| A E I M |
| B F J N |
| x y z 1| x = |Ax+By+Cz+D Ex+Fy+Gz+H Ix+Jy+Kz+L Mx+Ny+Oz+P|
| C G K O |
| D H L P |
可以看出两种变换方式是一个转置关系,结果只是形式上的不同,但这里我们使用后者,即右乘变换矩阵,因为比较普遍。
很显然,相机的变换可以分成两个阶段:旋转和平移。我们先来看旋转。
在空间坐标系中,相机旋转之前世界坐标系xyz和相机坐标系u0v0n0的各个轴向量的方向相同,有关系:
| u0 | | x |
P = |Pu0 Pv0 Pn0| x | v0 | = |Px Py Pz| x | y |
| n0 | | z |
这里P是空间坐标系中的一个向量。|u0 v0 n0|^T是相机基矩阵,|Pu0 Pv0 Pn0|是P在相机基矩阵下的坐标。|x y z|^T是
世界基矩阵,|Px Py Pz|是P在它下面的坐标。有Pu0 = Px, Pv0 =Py, Pn0 = Pz。
相机和向量P都旋转之后,有关系:
| u | | x |
P' = |Pu0 Pv0 Pn0| x | v | = |Px' Py' Pz'| x | y |
| n | | z |
P'是P同相机一起旋转后的向量。|u v n|^T是相机旋转后的基矩阵,|Pu0 Pv0 Pn0|是P'在它下面的坐标,因为P是和相机一起旋转的,所以坐标不变。|x y z|^T仍为世界基矩阵,|Px' Py' Pz'|是P'在它下面的坐标。
现在看
| u | | x |
|Pu0 Pv0 Pn0| x | v | = |Px' Py' Pz'| x | y |
| n | | z |
因为|x y z|^T为一个单位阵,且Pu0 = Px, Pv0 =Py, Pn0 = Pz。 所以得到
| u |
|Px Py Pz| x | v | = |Px' Py' Pz'|
| n |
即|Px Py Pz|和相机一起旋转后变成|Px' Py' Pz'|,即P x R = P',而旋转变换矩阵R就是:
| u |
| v |
| n |
写成标准4x4矩阵:
| ux uy uz 0|
| vx vy vz 0|
| nx ny nz 0|
| 0 0 0 1|
平移矩阵T很简单:
| 1 0 0 0 |
| 0 1 0 0 |
| 0 0 1 0 |
| x y z 1 |
则相机矩阵就是:
| ux uy uz 0 | | 1 0 0 0 |
| vx vy vz 0 | | 0 1 0 0 |
C = R x T = x
| nx ny nz 0 | | 0 0 1 0 |
| 0 0 0 1 | | x y z 1 |
它的逆矩阵,即相机的逆变换矩阵为
| 1 0 0 0 | | ux vx nx 0 | | ux vx nx 0 |
| 0 1 0 0 | | uy vy ny 0 | | uy vy ny 0 |
C^-1 = T^-1 x R^-1 = x =
| 0 0 1 0 | | uz nz nz 0 | | uz vz nz 0 |
| -x -y -z 1 | | 0 0 0 1 | |-T.u -T.v -T.n 1 |