本文参考文献:
http://blog.sina.com.cn/s/blog_557d254601018dfv.html
http://blog.csdn.net/candycat1992/article/details/41254799
http://blog.csdn.net/cppyin/article/details/6177742
一、四元数概念及运算
1. 四元数引入
将实数域扩充到复数域,并用复数来表示平面向量,用复数的加、乘运算表示平面向量的合成、伸缩和旋转变换,这些观念已经在中学课程中学过了。那么,很自然的问题就是,在三维,或更高维空间中是否也有复数的类似物?也就是说,像扩充实数那样,在复数域的基础上添加一个或几个新的元素,并且让它们跟原来的复数做加减乘除,是否就可以得到一个新的数集,并且其中的元素还可以像复数域那样做加、减、乘、除运算,并满足通常复数的那些运算律,包括加法和乘法的交换律与结合律、乘法对加法的分配律等待?更进一步,我们是否可以期望用这样的数来表示三维或更高维空间中的伸缩和旋转,就像用复数表示平面向量的伸缩旋转那样方便?
历史上有很多数学家试图寻找过三维的复数,但后来证明这样的三维复数是不存在的 ,即使不考虑空间旋转,只从代数角度来说,三维的复数域作为普通复数域的扩张域也是不存在的。
知道了复数不能推广到三维,我们把目光移向四维复数,即四元数。四元数是由爱尔兰数学家威廉·卢云·哈密顿在1843年发现的。复数推广到四元数,必须牺牲掉数域的某一条或几条性质,哈密尔顿抛弃了乘法交换律。
2. 四元数定义
四元数都是 1、i、j 和 k 的线性组合,一般可表示为 d + ai + bj + ck, a、b、c、d是实数。 如把四元数的集合考虑成多维实数空间的话,四元数就代表着一个四维空间,相对于复数为二维空间。
3.加乘运算
要把两个四元数相加只需将相类的系数加起来就可以,就像复数一样。至于乘法则可跟随以下的乘数表:
以上表格中,最左边的列表示被乘数,最上面行表示乘数。
即 i^2=j^2=k^2=-1,ij=k, ji=-k, jk=i, kj=-i, ki=j, ik=-j
可以立即验证加法交换律、结合律,以及等式 p+0=0+p=p,方程 p+x=0 恒有解,乘法结合律,还有乘法对加法的分配律都是成立的,只不过没有乘法交换律。
例如:
假设:x = 3 + i, y = 5i + j - 2k
那么:
x + y = 3 + 6i + j - 2k
xy =( {3 + i} )( {5i + j - 2k} ) = 15i + 3j - 6k + 5i^2 + ij - 2ik
= 15i + 3j - 6k - 5 + k + 2j = - 5 + 15i + 5j - 5k
四元数也可以表示为: q=[w,v], 其中v=(x,y,z)是矢量,w是标量。
如果向量部分外积为零,两个四元数的乘积就可交换了。
4.模的定义
四元数也是可以归一化的,四元数的单位化与Vector类似,首先(a^2+b^2+c^2+d^2)的平方根 称为四元数的模,即|q| = Norm(q)=sqrt(w^2 + x^2 + y^2 + z^2)。所以Normlize(q)=q/Norm(q)=q / sqrt(w2 + x2 + y2 + z2).
5.共轭定义
若p=w+xi+yj+zk=w+v,则p*=w-xi-yj-zk=w-v 表示与p实部相等,向量部分相反的四元数,称为p的共轭。
p.p* = p*.p = w^2 + x^2 + y^2 + z^2 = |p|^2
6.逆的定义
四元数的转置通过 被定义。
q×q-1 = 1
两边同时乘以q*:
q×q-1×q* = q*
因为:q×q* = |q|2
所以q-1=q* / |q|2.
可以注意到,如果q是一个单位四元数的话,那么q的倒数就等于q的共轭。即:q-1=q*
二、四元数与三维空间的旋转
我们要关心的是三维空间上任意的伸缩旋转变换是否可用四元数的乘积来表示,而这一点对四元数来说是完全能够胜任的。
如果已知一个三维空间的伸缩旋转的转轴方向、旋转角度和伸缩比例,来求相应的四元数,是比较容易的。
特别地,单位化的四元数用来描述旋转:
以原点为旋转中心,旋转的轴是(α, β, γ)( α^2 + β^2 + γ^2 = 1),
(右手系的坐标定义的话,望向向量(α, β, γ)的前进方向反时针) 转θ角的旋转,用四元数表示就是,Q = (cos(θ/2); α * sin(θ/2), β * sin(θ/2), γ * sin(θ/2))
四元数的乘法的意义类似于Matrix的乘法-可以将两个旋转合并,例如:
Q=Q1*Q2 表示Q的是先做Q2的旋转,再做Q1的旋转的结果,而多个四元数的旋转也是可以合并的,当有多次旋转操作时,使用四元数可以获得更高的计算效率。
补充:
一个四元数可以表示为q = w + xi + yj + zk,现在就来回答这样一个简单的式子是怎么和三维旋转结合在一起的。为了方便,我们下面使用q = ((x, y, z),w) = (v, w),其中v是向量,w是实数,这样的式子来表示一个四元数。
我们先来看问题的答案。我们可以使用一个四元数q=((x,y,z)sinθ2, cosθ2) 来执行一个旋转。具体来说,如果我们想要把空间的一个点P绕着单位向量轴u = (x, y, z)表示的旋转轴旋转θ角度,我们首先把点P扩展到四元数空间,即四元数p = (P, 0)。那么,旋转后新的点对应的四元数(当然这个计算而得的四元数的实部为0,虚部系数就是新的坐标)为:
其中,
我们举个最简单的例子:把点P(1, 0, 1)绕旋转轴u = (0, 1, 0)旋转90°,求旋转后的顶点坐标。首先将P扩充到四元数,即p = (P, 0)。而q = (u*sin45°, cos45°)。求的值。建议大家一定要在纸上计算一边,这样才能加深印象,连笔都懒得动的人还是不要往下看了。最后的结果p’ = ((1, 0, -1), 0),即旋转后的顶点位置是(1, 0, -1)。
如果想要得到复合旋转,只需类似复合矩阵那样左乘新的四元数,再进行运算即可。
我们来总结下四元数旋转的几个需要注意的地方:
1.用于旋转的四元数,每个分量的范围都在(-1,1);
2.每一次旋转实际上需要两个四元数的参与,即q和q*;
3.所有用于旋转的四元数都是单位四元数,即它们的模是1;
// Quaternion.cpp
#include <math.h>
#include <iostream.h>
/// Define Data type
typedef struct
{
double t; // real-component
double x; // x-component
double y; // y-component
double z; // z-component
} quaternion;
Bill 注:Kakezan 在日语里是 “乘法”的意思
quaternion Kakezan(quaternion left, quaternion right)
{
quaternion ans;
double d1, d2, d3, d4;
d1 = left.t * right.t;
d2 = -left.x * right.x;
d3 = -left.y * right.y;
d4 = -left.z * right.z;
ans.t = d1+ d2+ d3+ d4;
d1 = left.t * right.x;
d2 = right.t * left.x;
d3 = left.y * right.z;
d4 = -left.z * right.y;
ans.x = d1+ d2+ d3+ d4;
d1 = left.t * right.y;
d2 = right.t * left.y;
d3 = left.z * right.x;
d4 = -left.x * right.z;
ans.y = d1+ d2+ d3+ d4;
d1 = left.t * right.z;
d2 = right.t * left.z;
d3 = left.x * right.y;
d4 = -left.y * right.x;
ans.z = d1+ d2+ d3+ d4;
return ans;
}
Make Rotational quaternion
quaternion MakeRotationalQuaternion(double radian, double AxisX, double AxisY, double AxisZ)
{
quaternion ans;
double norm;
double ccc, sss;
ans.t = ans.x = ans.y = ans.z = 0.0;
norm = AxisX * AxisX + AxisY * AxisY + AxisZ * AxisZ;
if(norm <= 0.0) return ans;
norm = 1.0 / sqrt(norm);
AxisX *= norm;
AxisY *= norm;
AxisZ *= norm;
ccc = cos(0.5 * radian);
sss = sin(0.5 * radian);
ans.t = ccc;
ans.x = sss * AxisX;
ans.y = sss * AxisY;
ans.z = sss * AxisZ;
return ans;
}
Put XYZ into quaternion
quaternion PutXYZToQuaternion(double PosX, double PosY, double PosZ)
{
quaternion ans;
ans.t = 0.0;
ans.x = PosX;
ans.y = PosY;
ans.z = PosZ;
return ans;
}
/ main
int main()
{
double px, py, pz;
double ax, ay, az, th;
quaternion ppp, qqq, rrr;
cout << "Point Position (x, y, z) " << endl;
cout << " x = ";
cin >> px;
cout << " y = ";
cin >> py;
cout << " z = ";
cin >> pz;
ppp = PutXYZToQuaternion(px, py, pz);
while(1) {
cout << "\nRotation Degree ? (Enter 0 to Quit) " << endl;
cout << " angle = ";
cin >> th;
if(th == 0.0) break;
cout << "Rotation Axis Direction ? (x, y, z) " << endl;
cout << " x = ";
cin >> ax;
cout << " y = ";
cin >> ay;
cout << " z = ";
cin >> az;
th *= 3.1415926535897932384626433832795 / 180.0; /// Degree -> radian;
qqq = MakeRotationalQuaternion(th, ax, ay, az);
rrr = MakeRotationalQuaternion(-th, ax, ay, az);
ppp = Kakezan(rrr, ppp);
ppp = Kakezan(ppp, qqq);
cout << "\nAnser X = " << ppp.x
<< "\n Y = " << ppp.y
<< "\n Z = " << ppp.z << endl;
}
return 0;
}
三、关于插值
使用四元数的原因就是在于它非常适合插值,这是因为他是一个可以规格化的4维向量,最简单的插值算法就是线性插值,公式如:
q(t)=(1-t)q1+t q2
q(t)=((1-t)q1+tq2 )/ |(1-t)q1+tq2|
这样,插值向量v(t)的端点就会沿着v1,v2端点构成的圆弧行进。因为v1,v2是等长的,这个圆弧实际上是位于v1,v2构成的球面上的一段,所以又叫球面线性插值。
参考郑军“四元数插值算法实现游戏角色平滑旋转”一文:
// 插值计算
double x0,y0,z0,wO;
double xl,yl,zl,wl;
//参数t为插值变量
double t;
//在参数t下输出一个四元数,得到一个插值方位
double x,y,z,w;
//计算开始方位和结束方位的点积
double cosO=xO*xl+y0*yl+zO*zl+WO*wl;
if(cosO<0){
xO=一xO;
yo=一yO;
zO=一zO;
wO=一wo;
cosO=一eosO;
}
//进行插值运算
double kO,kl;
if(cosO>0.9999){
kO=1一t;
kl=t;
}else{
double sinO=Mathf.Sqrt(1一Mathf.Pow(cosO,2));
double O=Mathf.Atan2(sinO,cosO);
float one=1.O/sinO;
kO=Mathf.Sin((1一t) * O)*one;
kl=Mathf.Sin(t * O) * one;
}
//输出在t参数下的插值结果
x=xO * kO+xl * kl:
Z=zO * kO+zl * kl:
w=wo * kO+wl * kl: