OpenGL ES 从零开始系列9(完结)：四元数

在进入下一篇关于骨骼动画的文章之前，让我们先花点时间来了解一个马上会使用到的新数据类型：四元数[译者注:关于四元数的概念可以参考这个链接:点我]。我们用四元数存储单一骨骼在3个轴线上的旋转信息，换句话说，存储的是骨骼指向的方向。在下一部分介绍的仿真骨骼动画中，你将会看到，模型的顶点是同一个或多个骨骼相关联的，当骨骼移动时它们也会随之变化。相对于将欧拉角信息存储在3个GLfloats变量或一个 Vector3D 变量里来说, 使用四元数有2个优点：

1.四元数不会造成万向节死锁(gimbal lock)，但是欧拉角容易造成万向节死锁，使用四元数能够让我们的3D模型能够全方位的移动。
2.相比于给每个欧拉角做矩阵旋转转换计算，使用四元数结合多角度旋转可以显著的减少计算量。

从某些方面来看，四元数极其复杂且难于理解。它们是高级数学：完全疯狂的符咒。幸运的是，你不需要完全理解它们背后的数学含义。但是，我们现在需要使用它们来完成骨骼动画，所以还是值得我们花费些时间来讨论下它们的概念和怎么使用它们。

Discovery探索
从数学上讲,四元数是复数的一个扩展延伸，于1843年由Sir William Rowan Hamilton 发现。技术上讲，四元数表现为实数之上的4维正规可除代数。Zoiks!更简单的讲，四元数被认为是第四维度用来计算笛卡尔坐标中的3个坐标值。好吧，一切可能不那么简单，对吧？

先别怕，如果你不精通高等数学，四元数可能会让你头疼。但是，如我之前所说，如果你只是使用它们，完全不必深入了解。这玩意和你见过的一些概念是非常类似的。不知你是否还能想起我们在3维空间里涉及到的4X4矩阵的矩阵转换。当我们使用已转换的数据的时候，忽略了第4个值。我们可以把这里的第四个值当成四元数，为计算提供了一个位置。数学范畴内，请不要跟我说——过度简化有助于凡人在四元数世界里占有一席之地，有所作为。

四元数在探索时代里被认为是相当创新的，但最繁荣的时期却如此短暂。在1880中期，向量微积分开始在计算领域取代四元数理论，因为它用了一种更为容易理解和描述的概念描述了同样的现象。

Not Quite Dead Yet!虽死犹生
但在20世纪，四元数又重新获宠。正如我们在part 7里讨论的，有一个被称为gimbal lock 的现象，当你在每个轴线单独做旋转转换的时候就会发生，此现象的危害就是可能导致在三个轴中的一个轴上停止旋转。

尽管事实是四元数源于复数和理论数学，但它们都有实际应用。其中一个实际应用是三轴线上旋转角的展现。由于四元数用四个维度展示了笛卡尔（或三轴）旋转，此展现不会导致gimbal lock,而且你可以在四元数和旋转矩阵之间，四元数和欧拉角之间进行无损转换。这使得存储某些对象的旋转信息相当完美，比如。。。骨骼框架中的单独骨骼？不需要存贮3轴的角信息，而是存储一个单独的四元数。

四元数和矩阵一样，可以相乘，且存储于不同四元数中的旋转值通过相乘来合并计算。四元数乘积和2个旋转矩阵乘积的结果是完全一样的，考虑到减少计算量，这意味着除了要避免gimbal lock,还要减少每次程序循环运行的FLOPS(每秒浮点运算次数)。和矩阵乘法相比，四元数乘法不仅步骤少，而且可以通过一个四元数表达3轴所有数据。如果通过Vector3D 或3个GLfloats来存储旋转信息，我们经常不得不做3次矩阵乘法——每轴都要算一次。结论是，通过把存储旋转的独立角信息存为四元数，可以带来可观的性能提升。

 

The Quaternion Struct 四元数结构体
从数据上看来，四元数只不过是比Vector 3D多加了一个GLfloat，经常把它当成w字段。所以对我们来说一个四元数就象这样：
typedef struct {
   GLfloat x;
   GLfloat y;
   GLfloat z;
   GLfloat w;
} Quaternion3D;

Normalizing a Quaternion 四元数归一化
这个十分简单。四元数代表空间里的一个方向，就像Vector3Ds，实际距离的值并不在意，并且在进行一些计算之前／之后使他们正常下降到1.0。完成这些我们可以这样做：
static inline void Quaternion3DNormalize(Quaternion3D *quaternion)
{
   GLfloat magnitude;
  
   magnitude = sqrtf((quaternion->x * quaternion->x) +
                     (quaternion->y * quaternion->y) +
                     (quaternion->z * quaternion->z) +
                     (quaternion->w * quaternion->w));
  
   quaternion->x /= magnitude;
   quaternion->y /= magnitude;
   quaternion->z /= magnitude;
   quaternion->w /= magnitude;
}

Creating a Quaternion from a Rotation matrix 从一个旋转矩阵中创建一个四元数 
如果我们没法在四元数同其他的对象转换过程中保存旋转角度信息的话，四元数对我们来说还是没用的。首先我们从一个旋转矩阵中创建一个四元数。代码如下： 
static inline Quaternion3D Quaternion3DMakeWithMatrix3D(Matrix3D matrix)
{
   Quaternion3D quat;
   GLfloat trace, s;
  
   trace = matrix[0] + matrix[5] + matrix[10];
   if (trace > 0.0f)
   {
       s = sqrtf(trace + 1.0f);
       quat.w = s * 0.5f;
       s = 0.5f / s;
      
       quat.x = (matrix[9] - matrix[6]) * s;
       quat.y = (matrix[2] - matrix[8]) * s;
       quat.z = (matrix[4] - matrix[1]) * s;
   }
   else
   {
       NSInteger biggest;
       enum      {A,E,I};
       if (matrix[0] > matrix[5])
           if (matrix[10] > matrix[0])
               biggest = I;   
           else
               biggest = A;
       else
           if (matrix[10] > matrix[0])
               biggest = I;
           else
               biggest = E;
      
       switch (biggest)
       {
           case A:
               s = sqrtf(matrix[0] - (matrix[5] + matrix[10]) + 1.0f);
               if (s > QUATERNION_TRACE_ZERO_TOLERANCE)
               {
                   quat.x = s * 0.5f;
                   s = 0.5f / s;
                   quat.w = (matrix[9] - matrix[6]) * s;
                   quat.y = (matrix[1] + matrix[4]) * s;
                   quat.z = (matrix[2] + matrix[8]) * s;
                   break;
               }
               s = sqrtf(matrix[10] - (matrix[0] + matrix[5]) + 1.0f);
               if (s > QUATERNION_TRACE_ZERO_TOLERANCE)
               {
                   quat.z = s * 0.5f;
                   s = 0.5f / s;
                   quat.w = (matrix[4] - matrix[1]) * s;
                   quat.x = (matrix[8] + matrix[2]) * s;
                   quat.y = (matrix[9] + matrix[6]) * s;
                   break;
               }
               s = sqrtf(matrix[5] - (matrix[10] + matrix[0]) + 1.0f);
               if (s > QUATERNION_TRACE_ZERO_TOLERANCE)
               {
                   quat.y = s * 0.5f;
                   s = 0.5f / s;
                   quat.w = (matrix[2] - matrix[8]) * s;
                   quat.z = (matrix[6] + matrix[9]) * s;
                   quat.x = (matrix[4] + matrix[1]) * s;
                   break;
               }
               break;
              
           case E:
               s = sqrtf(matrix[5] - (matrix[10] + matrix[0]) + 1.0f);
               if (s > QUATERNION_TRACE_ZERO_TOLERANCE)
               {
                   quat.y = s * 0.5f;
                   s = 0.5f / s;
                   quat.w = (matrix[2] - matrix[8]) * s;
                   quat.z = (matrix[6] + matrix[9]) * s;
                   quat.x = (matrix[4] + matrix[1]) * s;
                   break;
               }
               s = sqrtf(matrix[10] - (matrix[0] + matrix[5]) + 1.0f);
               if (s > QUATERNION_TRACE_ZERO_TOLERANCE)
               {
                   quat.z = s * 0.5f;
                   s = 0.5f / s;
                   quat.w = (matrix[4] - matrix[1]) * s;
                   quat.x = (matrix[8] + matrix[2]) * s;
                   quat.y = (matrix[9] + matrix[6]) * s;
                   break;
               }
               s = sqrtf(matrix[0] - (matrix[5] + matrix[10]) + 1.0f);
               if (s > QUATERNION_TRACE_ZERO_TOLERANCE)
               {
                   quat.x = s * 0.5f;
                   s = 0.5f / s;
                   quat.w = (matrix[9] - matrix[6]) * s;
                   quat.y = (matrix[1] + matrix[4]) * s;
                   quat.z = (matrix[2] + matrix[8]) * s;
                   break;
               }
               break;
              
           case I:
               s = sqrtf(matrix[10] - (matrix[0] + matrix[5]) + 1.0f);
               if (s > QUATERNION_TRACE_ZERO_TOLERANCE)
               {
                   quat.z = s * 0.5f;
                   s = 0.5f / s;
                   quat.w = (matrix[4] - matrix[1]) * s;
                   quat.x = (matrix[8] + matrix[2]) * s;
                   quat.y = (matrix[9] + matrix[6]) * s;
                   break;
               }
               s = sqrtf(matrix[0] - (matrix[5] + matrix[10]) + 1.0f);
               if (s > QUATERNION_TRACE_ZERO_TOLERANCE)
               {
                   quat.x = s * 0.5f;
                   s = 0.5f / s;
                   quat.w = (matrix[9] - matrix[6]) * s;
                   quat.y = (matrix[1] + matrix[4]) * s;
                   quat.z = (matrix[2] + matrix[8]) * s;
                   break;
               }
               s = sqrtf(matrix[5] - (matrix[10] + matrix[0]) + 1.0f);
               if (s > QUATERNION_TRACE_ZERO_TOLERANCE)
               {
                   quat.y = s * 0.5f;
                   s = 0.5f / s;
                   quat.w = (matrix[2] - matrix[8]) * s;
                   quat.z = (matrix[6] + matrix[9]) * s;
                   quat.x = (matrix[4] + matrix[1]) * s;
                   break;
               }
               break;
              
           default:
               break;
       }
   }
   return quat;
}

好的，如果你想真正知道这里是怎么工作的。你需要了解矩阵运算、欧拉旋转理论、特征值和纹理。更不用说理解旋转在矩阵里的表示方法。当你完成你的数学博士学位，你可以回过头用它来解释剩余的usl。你会发现这个函数中使用的算法都是Matrix的FAQ上用伪代码列出来的。

 

Creating a Rotation Matrix from a Quaternion 从一个四元数中创建旋转矩阵
这另外的一个方法相对简单些。并且这个基本算法来自于Matrix FAQ，虽然我需要把它转换成行优先的顺序。
static inline void Matrix3DSetUsingQuaternion3D(Matrix3D matrix, Quaternion3D quat)
{
   matrix[0]  = (1.0f - (2.0f * ((quat.y * quat.y) + (quat.z * quat.z))));
   matrix[1]  = (2.0f * ((quat.x * quat.y) - (quat.z * quat.w)));
   matrix[2]  = (2.0f * ((quat.x * quat.z) + (quat.y * quat.w)));
   matrix[3] = 0.0f;
   matrix[4]  = (2.0f * ((quat.x * quat.y) + (quat.z * quat.w)));
   matrix[5]  = (1.0f - (2.0f * ((quat.x * quat.x) + (quat.z * quat.z))));
   matrix[6]  = (2.0f * ((quat.y * quat.z) - (quat.x * quat.w)));
   matrix[7] = 0.0f;
   matrix[8]  = (2.0f * ((quat.x * quat.z) - (quat.y * quat.w)));
   matrix[9]  = (2.0f * ((quat.y * quat.z) + (quat.x * quat.w)));
   matrix[10] = (1.0f - (2.0f * ((quat.x * quat.x) + (quat.y * quat.y))));
   matrix[11] = 0.0f;
   matrix[12]  = 0.0f;
   matrix[13]  = 0.0f;
   matrix[14] = 0.0f;
   matrix[15] = 1.0f;
}

Converting an Angle and Axis of Rotation to a Quaternion 把一个角度和旋转轴转换成一个四元数

四元数可以做的另外一种转换是，表示成在一个Vector3D表示的轴线上进行旋转。这在骨骼动画里面是非常有用的，因为这种表现形式通过矩阵是很难做到的。创建一个基于角度和轴旋转得四元数，我们可以这样做：
static inline Quaternion3D Quaternion3DMakeWithAxisAndAngle(Vector3D axis, GLfloat angle)
{
   Quaternion3D quat;
   GLfloat sinAngle;
  
   angle *= 0.5f;
    Vector3DNormalize(&axis);
   sinAngle = sinf(angle);
   quat.x = (axis.x * sinAngle);
   quat.y = (axis.y * sinAngle);
   quat.z = (axis.z * sinAngle);
   quat.w = cos(angle);
  
   return quat;
}

Extracting an Angle and Axis of Rotation from a Quaternion 从一个四元数中检测角度和轴得旋转
反过来，我们也可以从四元数中取得旋转的数据，包括旋转角度和深度，就像这样
static inline void Quaternion3DExtractAxisAndAngle(Quaternion3D quat, Vector3D *axis, GLfloat *angle)
{
   GLfloat s;
    Quaternion3DNormalize(&quat);
   s = sqrtf(1.0f - (quat.w * quat.w));
 
   if (fabs(s) < 0.0005f) s = 1.0f;
  
   if (axis != NULL)
   {
       axis->x = (quat.x / s);
       axis->y = (quat.y / s);
       axis->z = (quat.z / s);
   }
  
   if (angle != NULL)
       *angle = (acosf(quat.w) * 2.0f);
}

Quaternion Multiplication 四元数乘法
为了合并两种不同形式的四元数中得3D旋转信息。我们只需要让他们彼此相乘。好了继续我们得代码
static inline void Quaternion3DMultiply(Quaternion3D *quat1, Quaternion3D *quat2)
{
   Vector3D v1, v2, cp;
   float angle;
  
   v1.x = quat1->x;
   v1.y = quat1->y;
   v1.z = quat1->z;
   v2.x = quat2->x;
   v2.y = quat2->y;
   v2.z = quat2->z;
   angle = (quat1->w * quat2->w) - Vector3DDotProduct(v1, v2);
  
   cp = Vector3DCrossProduct(v1, v2);
   v1.x *= quat2->w;
   v1.y *= quat2->w;
   v1.z *= quat2->w;
   v2.x *= quat1->w;
   v2.y *= quat1->w;
   v2.z *= quat1->w;
  
   quat1->x = v1.x + v2.x + cp.x;
   quat1->y = v1.y + v2.y + cp.y;
   quat1->z = v1.z + v2.z + cp.z;
   quat1->w = angle;
}

Inverting a Quaternion 四元数转置
我们通过做一个四元数的共轭运算来取得四元数的转置。四元数做共轭运算其实就是将四元数中表示向量(x,y,z)的值取反。在这里的实现中，我们把它[四元数转置计算]作为四元数标准计算的一部分，而不是一个独立的步骤：
static inline void Quaternion3DInvert(Quaternion3D  *quat)
{
   GLfloat length = 1.0f / ((quat->x * quat->x) +
                    (quat->y * quat->y) +
                    (quat->z * quat->z) +
                    (quat->w * quat->w));
   quat->x *= -length;
   quat->y *= -length;
   quat->z *= -length;
   quat->w *= length;
}

Creating a Quaternion from Euler Angles 从欧拉角中创建四元数
前面我说过在旋转中最好不要使用欧拉角，但是有时候我们需要将欧拉角转换成四元数，比如说用户输入的信息是欧拉角信息。转换的步骤是，将欧拉轴用Vector3D表示出来，然后将Vector3D的值转换成四元数，最后将四元数相乘来得到结果：
static inline Quaternion3D Quaternion3DMakeWithEulerAngles(GLfloat x, GLfloat y, GLfloat z)
{
   Vector3D vx = Vector3DMake(1.f, 0.f, 0.f);
   Vector3D vy = Vector3DMake(0.f, 1.f, 0.f);
   Vector3D vz = Vector3DMake(0.f, 0.f, 1.f);
  
   Quaternion3D qx = Quaternion3DMakeWithAxisAndAngle(vx, x);
   Quaternion3D qy = Quaternion3DMakeWithAxisAndAngle(vy, y);
   Quaternion3D qz = Quaternion3DMakeWithAxisAndAngle(vz, z);


    Quaternion3DMultiply(&qx, &qy );
    Quaternion3DMultiply(&qx, &qz );
   return qx;
}

SLERPs and NLERPS

最后，最重要的部分来了：SLERPS 和 NLERPS。SLERPS 是 球面线性插值 的缩写，NLERPS是 归一化线性插值 的缩写，还记得我们在“Part 9a”是怎样计算 代数的 线性插值 的吗，我们由此计算出了每一个依赖于 过去了多长时间的点的准确位置。但是，在插值四元法中 通过这样的 计算我们得不到想要的结果。想象一个球和一个窝关节图：

知道哪里用得到这样的球窝关节吗？在你的  胳膊肘 的地方。当然这只是为了描叙方便。站在镜子的前面，将你的胳膊举过头顶，并保持你的胳膊竖直，然后将胳膊放下，放到你身体的一侧。
你的胳膊不是直线运动的，对吧？你的手是弧线运动的，并且比你的肘关节和你的手臂运动的要快，很有可能你的手并不是以一个恒定的速度运动。这就是我们在做动画旋转角 时想要模拟的基本的转动方式。

球面线性插值 和 归一化线性插值 是两种最通用的方法，球面线性插值 更常用，并且对现实提供了一个更好的模拟，但是它更慢，需要的配置高。归一化线性插值 更快，因为它和我们以前用过的线性插值本质上是一样的，仅仅是插值后作归一化，这个完全可以从“归一化线性插值”这个名字中体味出来。

当动画效果比对速度的要求更严格时，使用球面线性插值，否则就用 归一化线性插值，这两种方法完全相同的参数，在这两种方法之间转换应该没什么问题，你可以尝试看一下是否球面线性插值 要更多处理器消耗。
这里是一个简单的 归一化线性插值的例子：
static inline Quaternion3D Quaternion3DMakeWithNLERP(Quaternion3D *start, Quaternion3D *finish, GLclampf progress)
{
   Quaternion3D ret;
   GLfloat inverseProgress = 1.0f - progress;
   ret.x = (start->x * inverseProgress) + (finish->x * progress);    
   ret.y = (start->y * inverseProgress) + (finish->y * progress);
   ret.z = (start->z * inverseProgress) + (finish->z * progress);
   ret.w = (start->w * inverseProgress) + (finish->w * progress);
    Quaternion3DNormalize(&ret);
   return ret;
}

现在，球面线性插值的实现就有一点复杂了，它要能够 计算 四元数的点积 ，这跟计算一个向量的点积 是一样的：

static inline Quaternion3D Quaternion3DMakeWithSLERP(Quaternion3D *start, Quaternion3D *finish, GLclampf progress)
{
   GLfloat startWeight, finishWeight, difference;
   Quaternion3D ret;
  
   difference = ((start->x * finish->x) + (start->y * finish->y) + (start->z * finish->z) + (start->w * finish->w));
   if ((1.f - fabs(difference)) > .01f)
   {
       GLfloat theta, oneOverSinTheta;
      
       theta = acosf(fabsf(difference));
       oneOverSinTheta = (1.f / sinf(theta));
       startWeight = (sinf(theta * (1.f - progress)) * oneOverSinTheta);
       finishWeight = (sinf(theta * progress) * oneOverSinTheta);
       if (difference < 0.f)
           startWeight = -startWeight;
   } else
   {
       startWeight = (1.f - progress);
       finishWeight = progress;
   }
   ret.x = (start->x * startWeight) + (finish->x * finishWeight);
   ret.y = (start->y * startWeight) + (finish->y * finishWeight);
   ret.z = (start->z * startWeight) + (finish->z * finishWeight);
   ret.w = (start->w * startWeight) + (finish->w * finishWeight);
    Quaternion3DNormalize(&ret);
  
   return ret;
}

Finish Line

结束语：
现在我们有能力做一些仿真骨骼动画了，下次我们还会讲解这个东西，我已经用这些新的函数和数据类型更新了Xcode工程模板。
注：1、因为我们的程序使用了浮点数，实际上我们会看到由于“有效位丢失”造成轻微变化的想象，这不是因为我们使用四元法，是因为我们使用了浮点数。

 

原文：iphonedevelopment，OpenGL ES from the Ground Up 9 (intermission) Quaternions
iTyran翻译讨论地址：http://ityran.com/forum-36-1.html