HTML clipboard
HTML clipboard
四元数的对数、指数和标量乘运算
首先，让我们重写四元数的定义，引入一个新的变量α，等于半角θ/2：
α = θ/2
|| n || = 1
q = [cosα n sinα] = [cosα xsinα
ysinα zsinα]
q 的对数定义为公式10.15：
log q = log([cosα n sinα]) ≡ [0 α
n ]
公式10.15 四元数的对数
≡表示"恒等于"，注意log q 的结果，它一般不是单位四元数。
指数以严格相反的方式定义，首先，设四元数 p 的形式为[0 α n ]，
n 为单位向量：
p = [0 α n ] = [0 (αx αy αz)]
|| n || = 1
接着，指数定义为公式10.16:
exp p = exp([0 α n ]) = [cosα n
sinα]
公式10.16 四元数的指数
根据定义，exp p 总是返回单位四元数。
四元数的对数和指数类似于它们的标量形式，回忆一下，对于标量α，有下列关系成立:
同样，四元数指数运算为四元数对数运算的逆运算：
exp(log q ) = q
最后，四元数能与一个标量相乘。其计算方法非常直接：每个分量都乘以这个标量，给定标量k和四元数 q ，有公式 10.17：
k q = k[w v ] = [kw k v
] = k[w (x y z)] = [kw kx ky kz]
公式10.17 四元数和标量相乘
一般不会得到单位四元数，这也是为什么在表达角位移的场合中标量乘不是那么有用的原因。
四元数求幂
四元数能作为底数，记作 qt
（不要和指数运算混淆，指数运算只接受一个四元数作为参数，而四元数求幂有两个参数 ---- 四元数和指数）。四元数求幂的意义类似于实数求幂。回忆一下，a0
= 1， a1 = a，a为非零标量。当t从0变到1时，at从1到a。四元数求幂有类似的结论：当t从0变到1，
qt从[1, 0 ]到 q 。
这对四元数求幂非常有用，因为它可以从角位移中抽取"一部分"。例如，四元数 q
代表一个角位移，现在想要得到代表1/3这个角位移的四元数，可以这样计算： q1/3。
指数超出[0, 1]范围外的几何行为和预期的一样（但有一个重要的注意事项）。例如， q2代表的角位移是
q 的两倍。假设 q 代表绕x轴顺时针旋转30度，那么 q2代表绕x轴顺时针旋转60度，
q-1/3代表绕x轴逆时针旋转10度。
上面提到的注意事项是，四元数表达角位移时使用最短圆弧，不能"绕圈"。继续上面的例子， q4不是预期的绕x轴顺时针旋转240度，而是逆时针80度。显然，向一个方向旋转240度等价于向相反的方向旋转80度，都能得到正确的"最终结果"。但是，在此基础上的进一步运算，产生的就可能不是预期的结果了。例如，(q4)1/2不是 q 2，尽管我们感觉应该是这样。一般来说，凡是涉及到指数运算的代数公式，如(as)t
= a(st)，对四元数都不适用。
现在，我们已经理解四元数求幂可以为我们做什么了。让我们看看它的数学定义，四元数求幂定义在前一节讨论的"有用"运算上，定义如公式10.18：
注意，对于标量求幂，也有类似结论：
不难理解为什么当t从0变到1时 q'从单位四元数变到 q
。注意到对数运算只是提取了轴 n
和角度θ；接着，和指数t进行标量乘时，结果是θ乘以t；最后，指数运算"撤销"了对数运算，从tθ和 n 重新计算w和
v 。上面给出的定义就是标准数学定义，在理论上非常完美，但直接转换到代码却是很复杂的。程序清单10.1所示的代码展示了怎样计算
q'的值。
Listing 10.1: Code to raise a quaternion to a power
// Quaternion (input and output)
float w,x,y,z;
// Input exponent
float exponent;
// Check for the case of an identity quaternion.
// This will protect against divide by zero
if (fabs(w) < .9999f)
{
// Extract the half angle alpha (alpha = theta/2)
float alpha = acos(w);
// Compute new alpha value
float newAlpha = alpha * exponent;
// Compute new w value
w = cos(newAlpha);
// Compute new xyz values
float mult = sin(newAlpha) / sin(alpha);
x *= mult;
y *= mult;
z *= mult;
}
关于这些代码，需要注意的地方有:
（1）有必要做单位四元数的检查。因为w=+(-)1会导致mult的计算中出现除零现象。单位四元数的任意次方还是单位四元数。因此，如果检测到输入是单位四元数，忽略指数直接返回原四元数即可。
（2）计算alpha时，使用了acos函数，它的返回值是正的角度。这并不会违背一般性，任何四元数都能解释成有正方向的旋转角度，因为绕某轴的负旋转等价于绕指向相反方向的轴的正旋转。
四元数插值 ---- "slerp"
当今3D数学中四元数存在的理由是由于一种称作slerp的运算，它是球面线性插值的缩写（Spherical Linear
Interpolation）。slerp运算非常有用，因为它可以在两个四元数间平滑插值。slerp运算避免了欧拉角插值的所有问题。
slerp是一种三元运算，这意味着它有三个操作数。前两个操作数是两个四元数，将在它们中间插值。设这两个"开始"和"结束"四元数分别为 q 0和
q 1。插值参数设为变量t，t在0到1之间变化，slerp函数：slerp( q 0,
q 1, t)，将返回 q 0到 q 1之间的插值方位。
能否利用现有的数学工具推导出slerp公式呢？如果是在两个标量a0和a1间插值，我们会使用下面的标准线性插值公式：
Δa = a1 - a0
lerp(a0, a1, t) = a0 + tΔa
标准线性插值公式从a0开始，并加上a0和a1差的t倍，有三个基本步骤：
（1）计算两个值的差。
（2）取得差的一部分。
（3）在初始值上加上差的一部分。
可以使用同样的步骤在四元数间插值：
（1）计算两个值的差， q 0到 q 1的角位移由Δ q
= q 0-1 q 1给出。
（2）计算差的一部分，四元数求幂可以做到，差的一部分由(Δ q )t给出。
（3）在开始值上加上差的一部分，方法是用四元数乘法来组合角位移： q 0(Δ q )t
这样，得到slerp的公式如公式10.19所示：
这是理论上的slerp计算过程，实践中，将使用一种更加有效的方法。
我们在4D空间中解释四元数，因为所有我们感兴趣的四元数都是单位四元数，所以它们都"存在"于一个 4D"球面"上。
slerp的基本思想是沿着4D球面上连接两个四元数的弧插值(这就是球面线性插值这个名称的由来)。
可以把这种思想表现在平面上，设两个2D向量 v 0和 v 1，都是单位向量。我们要计算
v 1，它是沿 v 0到 v 1弧的平滑插值。设w是
v 0到 v 1弧所截的角，那么 v t就是绕
v 1沿弧旋转tw的结果，如图10.10所示：
将 v t表达成 v 0和 v 1的线性组合，从另一方面说，存在两个非零常数k0和k1，使得：
v t = k0 v 0
+ k1 v 1
可以用基本几何学求出k0和k1，图10.11展示了计算的方法：
对以k1 v 1为斜边的直角三角形应用三角公式得：
这里有两点需要考虑。第一，四元数 q 和- q
代表相同的方位，但它们作为slerp的参数时可能导致不一样的结果，这是因为4D球面不是欧式空间的直接扩展。而这种现象在2D和3D中不会发生。解决方法是选择
q 0和 q 1的符号使得点乘 q 0 . q
1的结果是非负。第二个要考虑的是如果 q 0和 q
1非常接近，sinθ会非常小，这时除法可能会出现问题。为了避免这样的问题，当sinθ非常小时使用简单的线性插值。程序清单10.2把所有的建议都应用到了计算四元数的slerp中：
Listing 10.2: How slerp is computed in practice
// The two input quaternions
float w0,x0,y0,z0;
float w1,x1,y1,z1;
// The interpolation parameter
float t;
// The output quaternion will be computed here
float w,x,y,z;
// Compute the "cosine of the angle" between the
// quaternions, using the dot product
float cosOmega = w0*w1 + x0*x1 + y0*y1 + z0*z1;
// If negative dot, negate one of the input
// quaternions to take the shorter 4D "arc"
if (cosOmega < 0.0f) {
w1 = –w1;
x1 = –x1;
y1 = –y1;
z1 = –z1;
cosOmega = –cosOmega;
}
// Check if they are very close together to protect
// against divide-by-zero
float k0, k1;
if (cosOmega > 0.9999f) {
// Very close - just use linear interpolation
k0 = 1.0f–t;
k1 = t;
} else {
// Compute the sin of the angle using the
// trig identity sin^2(omega) + cos^2(omega) = 1
float sinOmega = sqrt(1.0f – cosOmega*cosOmega);
// Compute the angle from its sin and cosine
float omega = atan2(sinOmega, cosOmega);
// Compute inverse of denominator, so we only have
// to divide once
float oneOverSinOmega = 1.0f / sinOmega;
// Compute interpolation parameters
k0 = sin((1.0f – t) * omega) * oneOverSinOmega;
k1 = sin(t * omega) * oneOverSinOmega;
}
// Interpolate
w = w0*k0 + w1*k1;
x = x0*k0 + x1*k1;
y = y0*k0 + y1*k1;
z = z0*k0 + z1*k1;