四元数

为什么使用四元数

为了回答这个问题，先来看看一般关于旋转（面向）的描述方法－欧拉描述法。它使用最简单的x,y,z值来分别表示在x，y，z轴上的旋转角度，其取值为0-360(或者0-2pi），一般使用roll，pitch，yaw来表示这些分量的旋转值。需要注意的是，这里的旋转是针对世界坐标系说的，这意味着第一次的旋转不会影响第二、三次的转轴，简单的说，三角度系统无法表现任意轴的旋转，只要一开始旋转，物体本身就失去了任意轴的自主性，这也就导致了万向轴锁（Gimbal Lock）的问题。

还有一种是轴角的描述方法（即我一直以为的四元数的表示法），这种方法比欧拉描述要好，它避免了Gimbal Lock，它使用一个3维向量表示转轴和一个角度分量表示绕此转轴的旋转角度，即(x,y,z,angle)，一般表示为(x,y,z,w)或者(v,w)。但这种描述法却不适合插值。

那到底什么是Gimbal Lock呢？正如前面所说，因为欧拉描述中针对x,y,z的旋转描述是世界坐标系下的值，所以当任意一轴旋转90°的时候会导致该轴同其他轴重合，此时旋转被重合的轴可能没有任何效果，这就是Gimbal Lock，这里有个例子演示了Gimbal Lock，点击这里下载。运行这个例子，使用左右箭头改变yaw为90°，此时不管是使用上下箭头还是Insert、Page Up键都无法改变Pitch，而都是改变了模型的roll。

那么轴、角的描述方法又有什么问题呢？虽然轴、角的描述解决了Gimbal Lock，但这样的描述方法会导致差值不平滑，差值结果可能跳跃，欧拉描述同样有这样的问题。



什么是四元数

四元数一般定义如下：

q=w+xi+yj+zk

其中w是实数，x,y,z是虚数，其中:

i*i=-1

j*j=-1

k*k=-1

也可以表示为：

q=[w,v]

其中v=(x,y,z)是矢量，w是标量，虽然v是矢量，但不能简单的理解为3D空间的矢量，它是4维空间中的的矢量，也是非常不容易想像的。

四元数也是可以归一化的，并且只有单位化的四元数才用来描述旋转（面向），四元数的单位化与Vector类似，

首先||q|| = Norm(q)=sqrt(w2 + x2 + y2 + z2)

因为w2 + x2 + y2 + z2=1

所以Normlize(q)=q/Norm(q)=q / sqrt(w2 + x2 + y2 + z2)

说了这么多，那么四元数与旋转到底有什么关系？我以前一直认为轴、角的描述就是四元数，如果是那样其与旋转的关系也不言而喻，但并不是这么简单，轴、角描述到四元数的转化：

w = cos(theta/2)

x = ax * sin(theta/2)

y = ay * sin(theta/2)

z = az * sin(theta/2)

其中(ax,ay,az)表示轴的矢量，theta表示绕此轴的旋转角度，为什么是这样？和轴、角描述到底有什么不同？这是因为轴角描述的“四元组”并不是一个空间下的东西，首先（ax,ay,az）是一个3维坐标下的矢量，而theta则是级坐标下的角度，简单的将他们组合到一起并不能保证他们插值结果的稳定性，因为他们无法归一化，所以不能保证最终插值后得到的矢量长度（经过旋转变换后两点之间的距离）相等，而四元数在是在一个统一的4维空间中，方便归一化来插值，又能方便的得到轴、角这样用于3D图像的信息数据，所以用四元数再合适不过了。



关于四元数的运算法则和推导这里有篇详细的文章介绍，重要的是一点，类似与Matrix的四元数的乘法是不可交换的，四元数的乘法的意义也类似于Matrix的乘法－可以将两个旋转合并，例如：

Q=Q1*Q2

表示Q的是先做Q2的旋转，再做Q1的旋转的结果，而多个四元数的旋转也是可以合并的，根据四元数乘法的定义，可以算出两个四元数做一次乘法需要16次乘法和加法，而3x3的矩阵则需要27运算，所以当有多次旋转操作时，使用四元数可以获得更高的计算效率。



为什么四元数可以避免Gimbal Lock

在欧拉描述中，之所以会产生Gimbal Lock是因为使用的三角度系统是依次、顺序变换的，如果在OGL中，代码可能这样：

glRotatef( angleX, 1, 0, 0)

glRotatef( angleY, 0, 1, 0)

glRotatef( angleZ, 0, 0, 1)



注意：以上代码是顺序执行，而使用的又是统一的世界坐标，这样当首先旋转了Y轴后，Z轴将不再是原来的Z轴，而可能变成X轴，这样针对Z的变化可能失效。

而四元数描述的旋转代码可能是这样：

TempQ = From Eula(x,y,z)

FinalQ =CameraQ * NewQ

theta, ax, ay, az = From (FinalQ)

glRotatef(theta, ax, ay, az);

其中(ax,ay,az)描述一条任意轴，theta描述了绕此任意轴旋转的角度，而所有的参数都来自于所有描述旋转的四元数做乘法之后得到的值，可以看出这样一次性的旋转不会带来问题。这里有个例子演示了使用四元数不会产生Gimbal Lock的问题。



关于插值

使用四元数的原因就是在于它非常适合插值，这是因为他是一个可以规格化的4维向量，最简单的插值算法就是线性插值，公式如：

q(t)=(1-t)q1+t q2

但这个结果是需要规格化的，否则q(t)的单位长度会发生变化，所以

q(t)=(1-t)q1+t q2 / || (1-t)q1+t q2 ||

如图:

 

尽管线性插值很有效，但不能以恒定的速率描述q1到q2之间的曲线，这也是其弊端，我们需要找到一种插值方法使得q1->q(t)之间的夹角θ是线性的，即θ(t)=(1-t)θ1+t*θ2，这样我们得到了球形线性插值函数q(t)，如下：

q(t)=q1 * sinθ(1-t)/sinθ + q2 * sinθt/sineθ

如果使用D3D，可以直接使用D3DXQuaternionSlerp函数就可以完成这个插值过程。





第二篇:
四元数入门
四元数常常可以在3D的书上看到。
但我的那本3D图形学书上，在没讲四元数是干什么的之前，就列了几张纸的公式，
大概因为自己还在上高中，不知道的太多，看了半天没看懂。。。
终于，在gameres上看到了某强人翻译的一个“4元数宝典 ”（原文是日本人写的。。。），感觉很好，分享下。
★旋转篇：
　我将说明使用了四元数（si yuan shu, quaternion）的旋转的操作步骤
（１）四元数的虚部，实部和写法
所谓四元数，就是把4个实数组合起来的东西。
4个元素中，一个是实部，其余3个是虚部。
比如，叫做Q的四元数，实部t而虚部是x,y,z构成，则像下面这样写。
Q = (t; x, y, z)
又，使用向量 V=(x,y,z)，
Q = (t; V)
也可以这么写。

正规地用虚数单位i,j,k的写法的话，
Q = t + xi + yj + zk
也这样写，不过，我不大使用

（２）四元数之间的乘法
虚数单位之间的乘法
ii = -1, ij = -ji = k (其他的组合也是循环地以下同文)
有这么一种规则。（我总觉得，这就像是向量积（外积），对吧）
用这个规则一点点地计算很麻烦，所以请用像下面这样的公式计算。

A = (a; U)
B = (b; V)
AB = (ab - U•V; aV + bU + U×V)
不过，“U•V”是内积，「U×V」是外积的意思。
注意：一般AB<>BA所以乘法的左右要注意！

（3）3次元的坐标的四元数表示
如要将某坐标(x,y,z)用四元数表示，
P = (0; x, y, z)
则要这么写。
 
另外，即使实部是零以外的值，下文的结果也一样。用零的话省事所以我推荐。

（４）旋转的四元数表示
以原点为旋转中心，旋转的轴是(α, β, γ)
（但 α^2 + β^2 + γ^2 = 1），
（右手系的坐标定义的话，望向向量(α, β, γ)的前进方向反时针地）
转θ角的旋转，用四元数表示就是，
Q = (cos(θ/2); α sin(θ/2), β sin(θ/2), γ sin(θ/2))
R = (cos(θ/2); -α sin(θ/2), -β sin(θ/2), -γ sin(θ/2))
(另外R 叫 Q 的共轭四元数。）

那么，如要实行旋转，
则 R P Q = (0; 答案)

请像这样三明治式地计算。这个值的虚部就是旋转之后的点的坐标值。
 （另外，实部应该为零。请验算看看）

例子代码
/// Quaternion.cpp
/// (C) Toru Nakata, [e-mail]
/// 2004 Dec 29
  
#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;
}





又一篇

在以前涉及的程序中，处理物体的旋转通常是用的矩阵的形式。由于硬件在纹理映射和光栅化上的加强，程序员可以将更多的CPU周期用于物理模拟等工作，这样将使得程序更为逼真。
一，原文出处： http://www.gamasutra.com/features/19980703/quaternions_01.htm
二，摘录：
有三种办法表示旋转：矩阵表示，欧拉角表示，以及四元组表示。矩阵，欧拉角表示法在处理插值的时候会遇到麻烦。
There are many ways to represent the orientation of an object. Most programmers use 3x3 rotation matrices or three Euler angles to store this information. Each of these solutions works fine until you try to smoothly interpolate between two orientations of an object.

矩阵表示法不适合于进行插值(在转动的前后两个朝向之间取得瞬时朝向)。矩阵有9个自由度（3×3矩阵），而实际上表示一个旋转只需要3个自由度（旋转轴3）。
Rotations involve only three degrees of freedom (DOF), around the x, y, and z coordinate axes. However, nine DOF (assuming 3x3 matrices) are required to constrain the rotation - clearly more than we need.
Another shortcoming of rotation matrices is that they are extremely hard to use for interpolating rotations between two orientations. The resulting interpolations are also visually very jerky, which simply is not acceptable in games any more.
欧拉角表示法
You can also use angles to represent rotations around three coordinate axes. You can write this as (q, c, f); simply stated, "Transform a point by rotating it counterclockwise about the z axis by q degrees, followed by a rotation about the y axis by c degrees, followed by a rotation about the x axis by f degrees."
欧拉角表示法的缺陷：把一个旋转变成了一系列的旋转。并且，进行插值计算也不方便。
However, there is no easy way to represent a single rotation with Euler angles that corresponds to a series of concatenated rotations. Furthermore, the smooth interpolation between two orientations involves numerical integration

目录
[隐藏]
• 1 四元组表示法
• 2 四元组的基本运算法则
• 3 用四元组进行旋转的公式
• 4 旋转叠加性质

[编辑]
四元组表示法
[w, v]其中v是矢量，表示旋转轴。w标量，表示旋转角度。所以，一个四元组即表示一个完整的旋转。
There are several notations that we can use to represent quaternions. The two most popular notations are complex number notation (Eq. 1) and 4D vector notation (Eq. 2).
w + xi + yj + zk (where i2 = j2 = k2 = -1 and ij = k = -ji with real w, x, y, z) (Eq. 1)
[w, v] (where v = (x, y, z) is called a "vector" and w is called a "scalar") (Eq. 2)
4D空间以及单元四元组：
Each quaternion can be plotted in 4D space (since each quaternion is comprised of four parts), and this space is called quaternion space. Unit quaternions have the property that their magnitude is one and they form a subspace, S3, of the quaternion space. This subspace can be represented as a 4D sphere. (those that have a one-unit normal), since this reduces the number of necessary operations that you have to perform.

[编辑]
四元组的基本运算法则
Table 1. Basic operations using quaternions.
Addition: q + q´ = [w + w´, v + v´]
Multiplication: qq´ = [ww´ - v • v´, v x v´ + wv´ +w´v] (• is vector dot product and x is vector cross product); Note: qq´ ? q´q
//为什么是这样？－－定义成这样滴，木有道理可以讲
Conjugate: q* = [w, -v] 共轭
Norm: N(q) = w2 + x2 + y2 + z2 模
Inverse: q-1 = q* / N(q)
Unit Quaternion: q is a unit quaternion if N(q)= 1 and then q-1 = q*
Identity: [1, (0, 0, 0)] (when involving multiplication) and [0, (0, 0, 0)] (when involving addition)

[编辑]
用四元组进行旋转的公式
重要!!!!：只有单元四元组才表示旋转。为什么????－－已解决（shoemake有详细证明）
It is extremely important to note that only unit quaternions represent rotations, and you can assume that when I talk about quaternions, I'm talking about unit quaternions unless otherwise specified.
Since you've just seen how other methods represent rotations, let's see how we can specify rotations using quaternions. It can be proven (and the proof isn't that hard) that the rotation of a vector v by a unit quaternion q can be represented as
v´ = q v q-1 (where v = [0, v]) (Eq. 3)

////////////////////////////////////////////////////////////
四元组和旋转参数的变换一个四元组由(x,y,z,w)四个变量表达,假设给定一个旋转轴axis和角度angle,那么这个四元组的计算公式为:

void Quaternion::fromAxisAngle(const Vector3 &axis, Real angle)
{
       Vector3 u = axis;
       u.normalize();
       Real s = Math::rSin(angle/2.f);
       x = s*u.x;
       y = s*u.y;
       z = s*u.z;
       w = Math::rCos(angle/2.f);
}

由四元组反算回旋转轴axis和角度angle的公式为:
/
void Quaternion::toAxisAngle(Vector3 &axis, Real &angle) const
{
       angle = acos(w)*2
       axis.x = x
       axis.y = y
       axis.z = z
}