四元数定义三维旋转

        四元数（Quaternion）是由爱尔兰数学家 哈密顿(Hamilton）在1843年发明的概念。四元数的乘法不符合交换律（commutative law）。

        四元数 可以描述为 (x *
i, y * j,
z * k, w) ，其中（x,y,z）表示虚部，w表示实部。

        四元数 最大的价值在于 描述旋转，解决三维中的 万向节锁 的问题。

万向节锁：

       三维中的旋转通常用绕 XYZ 三个轴的旋转来表示，这个旋转角度即 欧拉角，也就是经常说的 Pitch、Yaw、Roll。


        采用欧拉角 描述旋转的 优点是直观简单（矩阵乘），容易理解，旋转过程分解为：按照一定的顺序（比如X、Y、Z） 依次独立地绕轴旋转。

        这种独立的旋转带来了一个问题，也就是 未考虑到旋转之间的关联性。

        假设 物体先绕转X轴旋转，然后再绕Y轴，最后绕Z轴选择，当绕Y轴旋转90度之后
万向节锁 的问题就出现了，因为X轴已经被求值了，它不再随同其他两个轴旋转，这样X轴与Z轴就指向同一个方向（它们相当于同一个轴了）。

欧拉角 与 四元数

        欧拉角可以和四元数进行相互转换，一个四元数可以 由欧拉角（对应 φ,θ,ψ）变换得到：

        


        同样，欧拉角 也可以通过四元数 变换得到：

        


        代码可以参考作者定义的 四元数 类（Quaternion.h）

/* 四元数定义 - 用于旋转
linolzhang, 2010.2.10
*/

#ifndef QUATERNION_H
#define QUATERNION_H

#include "BaseOperation.h"
#include "Matrix4.h"

template <typename T>
class Quaternion
{
public:
T x, y, z, w; // (x,y,z)(i,j,k) + w - 虚部（x,y,z），实部（w）

public:
// --------------------------------------------------------
// 构造函数
Quaternion()
{
x = y = z = 0; // 初始化成员变量
w = 1;
}
Quaternion(T _x,T _y,T _z,T _w)
{
x = _x; y = _y; z = _z;
w = _w;
}
Quaternion(const Quaternion<T> &q) // 拷贝构造函数
{
x = q.x; y = q.y; z = q.z;
w = q.w;
}
const Quaternion &operator=(const Quaternion<T> &q) // 赋值操作符
{
x = q.x; y = q.y; z = q.z;
w = q.w;
return *this;
}
// --------------------------------------------------------
// 构造函数 - 由其它转换
Quaternion(const Tuple3<T> &rAxisVector, float rAngle)  // 由 旋转轴向量，旋转角 得到
{
w = cos(0.5*rAngle);

float t = sin(0.5*rAngle);
x = rAxisVector.x * t;
y = rAxisVector.y * t;
z = rAxisVector.z * t;
}
Quaternion(const Matrix4<T> &mat) // 由 4*4 矩阵转换得到
{
T trace = mat.m[0] + mat.m[5] + mat.m[10]; // 用来保存矩阵的 迹(trace)
float s = 0;

if (trace > 0) // 如果矩阵的trace>0
{
float temp = sqrt(trace + 1);
s = 0.5 / temp;

w = 0.5 * temp;
x = (mat.m[9] - mat.m[6]) * s;
y = (mat.m[2] - mat.m[8]) * s;
z = (mat.m[4] - mat.m[1]) * s;
}
else // 矩阵的trace<0
{
T tq[3];

tq[0] = 1 + mat.m[0] - mat.m[5] - mat.m[10];
tq[1] = 1 - mat.m[0] + mat.m[5] - mat.m[10];
tq[2] = 1 - mat.m[0] - mat.m[5] + mat.m[10];

// 查找最大的
int maxPos = 0;
for(int i=1; i<3; i++)
maxPos = tq[i]>tq[maxPos] ? i : maxPos;

float temp = sqrt(tq[maxPos]);
if( temp!= 0)
s = 0.5 / temp;

// 检查对角线
if (maxPos == 0)
{
w = (mat.m[9] - mat.m[6]) * s;
x = 0.5 * temp;
y = (mat.m[4] + mat.m[1]) * s;
z = (mat.m[2] + mat.m[8]) * s;
}
else if (maxPos == 1)
{
w = (mat.m[2] - mat.m[8]) * s;
x = (mat.m[4] + mat.m[1]) * s;
y = 0.5 * temp;
z = (mat.m[9] + mat.m[6]) * s;
}
else /* if (maxPos == 2) */
{
w = (mat.m[4] - mat.m[1]) * s;
x = (mat.m[2] + mat.m[8]) * s;
y = (mat.m[9] + mat.m[6]) * s;
z = 0.5 * temp;
}
}
}

// --------------------------------------------------------
// 重载 +=, -=, *=
void operator +=(const Quaternion<T> &q)
{
w += q.w;
x += q.x;
y += q.y;
z += q.z;
}
void operator -=(const Quaternion<T> &q)
{
w -= q.w;
x -= q.x;
y -= q.y;
z -= q.z;
}
void operator *=(const Quaternion<T> &q)
{
T w1 = w*q.w - x*q.x - y*q.y - z*q.z;
T x1 = w*q.x + x*q.w + y*q.z - z*q.y;
T y1 = w*q.y + y*q.w + z*q.x - x*q.z;
T z1 = w*q.z + z*q.w + x*q.y - y*q.x;

w = w1;
w = x1; y = y1; z = z1;
}

// --------------------------------------------------------
// 转换为矩阵表示
const Matrix4<T> convertToMatrix4()const
{
//  注：
//     [ 1-2y2-2z2 , 2xy-2wz , 2xz+2wy ]
//     [ 2xy+2wz , 1-2x2-2z2 , 2yz-2wx ]
//     [ 2xz-2wy , 2yz+2wx , 1-2x2-2y2 ]

T xx = x*x;  T xy = x*y;  T xz = x*z;  T xw = x*w;
T yy = y*y;  T yz = y*z;  T yw = y*w;
T zz = z*z;  T zw = z*w;

return Matrix4<T>(  1-2*(yy+zz),  2*(xy-zw),    2*(xz+yw),    0,
2*(xy+zw),    1-2*(xx+zz),  2*(yz-xw),    0,
2*(xz-yw),    2*(yz+xw),    1-2*(xx+yy),  0,
0,            0,            0,            1  );

// 如果使用规范后的四元数 - 单位四元数
// 代入 w2 = 1-x2-y2-z2 得到
//     [ w2+x2-y2-z2 , 2xy-2wz , 2xz+2wy ]
//     [ 2xy+2wz , w2-x2+y2-z2 , 2yz-2wx ]
//     [ 2xz-2wy , 2yz+2wx , w2-x2-y2+z2 ]
}
const Quaternion getConjugate()const // 得到共轭四元数
{
return Quaternion<T> (-x, -y, -z, w);
}

// --------------------------------------------------------
// 运算操作
void setIdentity() // 设置初始化值
{
x = y = z = 0;
w = 1;
}
void normalize() // 单位化
{
float magnitude = sqrtf(x*x + y*y + z*z + w*w);
if (magnitude != 0)
{
x /= magnitude;
y /= magnitude;
z /= magnitude;
w /= magnitude;
}
}
};

////////////////////////////////////////////////////////////////////////////////////////////////////////
// --------------------------------------------------------
// 外部运算，重载*
template <typename T>
inline const Quaternion<T> operator*(const Quaternion<T> &left,const Quaternion<T> &right)
{
return Quaternion<T>(   left.w * right.x + left.x * right.w + left.y * right.z - left.z * right.y,
left.w * right.y - left.x * right.z + left.y * right.w + left.z * right.x,
left.w * right.z + left.x * right.y - left.y * right.x + left.z * right.w,
left.w * right.w - left.x * right.x - left.y * right.y - left.z * right.z  );
}

// --------------------------------------------------------
// 外部运算 - 插值(q1->q2) t[0,1]   注：接受单位四元数
template<typename T>
const Quaternion<T> Slerp(const Quaternion<T> &q1, const Quaternion<T> &q2, float t)
{
// 四元数插值公式
// 线性插值 q(t) = q1 * (1-t) + q2 * t
// 球面插值 q(t) = q1 * sin( (1-t)*θ) ) /sinθ + q2 * sin( t*θ))/sinθ

float theta;           // 表示四元数之间的夹角（实际上应该是 theta）
float scale0, scale1;  // 缩放因子

float qi[4]; // 结果暂存

float cosTheta = q1.x*q2.x + q1.y*q2.y + q1.z*q2.z + q1.w*q2.w;  // 点乘 = cos(theta)

if (cosTheta < 0)
{
cosTheta = -cosTheta;
qi[0] = -q2.x;
qi[1] = -q2.y;
qi[2] = -q2.z;
qi[3] = -q2.w;
}
else
{
qi[0] = q2.x;
qi[1] = q2.y;
qi[2] = q2.z;
qi[3] = q2.w;
}

// If the quaternions are really close, do a simple linear interpolation.
if ( 1-cosTheta <= 0.0001 ) // 线性插值公式
{
scale0 = 1 - t;
scale1 = t;
}
else // 球面插值公式
{
theta  = (float) acos( cosTheta );
float temp  = 1.0f/sinf(theta);

scale0 = sinf((1 - t) * theta) * temp;
scale1 = sinf(t * theta)  * temp;
}

// 根据公式得到四元数结果
return Quaternion( scale0 * q1.x + scale1 * qi[0],
scale0 * q1.y + scale1 * qi[1],
scale0 * q1.z + scale1 * qi[2],
scale0 * q1.w + scale1 * qi[3] );
}

typedef Quaternion<float>  Quaternionf;
typedef Quaternion<double> Quaterniond;

#endif

        四元数 关键应用在于其表示旋转的插值，通过 四元数插值 能够有效进行平滑，在 OSG、OGRE、Unity、CryEngine、UnReal等引擎中早已成熟应用，其中的开源代码部分大家可以借鉴。