HTML clipboard
HTML clipboard
四元數(shù)的對(duì)數(shù)、指數(shù)和標(biāo)量乘運(yùn)算
首先，讓我們重寫四元數(shù)的定義，引入一個(gè)新的變量α，等于半角θ/2：
α = θ/2
|| n || = 1
q = [cosα n sinα] = [cosα xsinα
ysinα zsinα]
q 的對(duì)數(shù)定義為公式10.15：
log q = log([cosα n sinα]) ≡ [0 α
n ]
公式10.15 四元數(shù)的對(duì)數(shù)
≡表示"恒等于"，注意log q 的結(jié)果，它一般不是單位四元數(shù)。
指數(shù)以嚴(yán)格相反的方式定義，首先，設(shè)四元數(shù) p 的形式為[0 α n ]，
n 為單位向量：
p = [0 α n ] = [0 (αx αy αz)]
|| n || = 1
接著，指數(shù)定義為公式10.16:
exp p = exp([0 α n ]) = [cosα n
sinα]
公式10.16 四元數(shù)的指數(shù)
根據(jù)定義，exp p 總是返回單位四元數(shù)。
四元數(shù)的對(duì)數(shù)和指數(shù)類似于它們的標(biāo)量形式，回憶一下，對(duì)于標(biāo)量α，有下列關(guān)系成立:
同樣，四元數(shù)指數(shù)運(yùn)算為四元數(shù)對(duì)數(shù)運(yùn)算的逆運(yùn)算：
exp(log q ) = q
最后，四元數(shù)能與一個(gè)標(biāo)量相乘。其計(jì)算方法非常直接：每個(gè)分量都乘以這個(gè)標(biāo)量，給定標(biāo)量k和四元數(shù) q ，有公式 10.17：
k q = k[w v ] = [kw k v
] = k[w (x y z)] = [kw kx ky kz]
公式10.17 四元數(shù)和標(biāo)量相乘
一般不會(huì)得到單位四元數(shù)，這也是為什么在表達(dá)角位移的場(chǎng)合中標(biāo)量乘不是那么有用的原因。
四元數(shù)求冪
四元數(shù)能作為底數(shù)，記作 qt
（不要和指數(shù)運(yùn)算混淆，指數(shù)運(yùn)算只接受一個(gè)四元數(shù)作為參數(shù)，而四元數(shù)求冪有兩個(gè)參數(shù) ---- 四元數(shù)和指數(shù)）。四元數(shù)求冪的意義類似于實(shí)數(shù)求冪。回憶一下，a0
= 1， a1 = a，a為非零標(biāo)量。當(dāng)t從0變到1時(shí)，at從1到a。四元數(shù)求冪有類似的結(jié)論：當(dāng)t從0變到1，
qt從[1, 0 ]到 q 。
這對(duì)四元數(shù)求冪非常有用，因?yàn)樗梢詮慕俏灰浦谐槿?一部分"。例如，四元數(shù) q
代表一個(gè)角位移，現(xiàn)在想要得到代表1/3這個(gè)角位移的四元數(shù)，可以這樣計(jì)算： q1/3。
指數(shù)超出[0, 1]范圍外的幾何行為和預(yù)期的一樣（但有一個(gè)重要的注意事項(xiàng)）。例如， q2代表的角位移是
q 的兩倍。假設(shè) q 代表繞x軸順時(shí)針旋轉(zhuǎn)30度，那么 q2代表繞x軸順時(shí)針旋轉(zhuǎn)60度，
q-1/3代表繞x軸逆時(shí)針旋轉(zhuǎn)10度。
上面提到的注意事項(xiàng)是，四元數(shù)表達(dá)角位移時(shí)使用最短圓弧，不能"繞圈"。繼續(xù)上面的例子， q4不是預(yù)期的繞x軸順時(shí)針旋轉(zhuǎn)240度，而是逆時(shí)針80度。顯然，向一個(gè)方向旋轉(zhuǎn)240度等價(jià)于向相反的方向旋轉(zhuǎn)80度，都能得到正確的"最終結(jié)果"。但是，在此基礎(chǔ)上的進(jìn)一步運(yùn)算，產(chǎn)生的就可能不是預(yù)期的結(jié)果了。例如，(q4)1/2不是 q 2，盡管我們感覺(jué)應(yīng)該是這樣。一般來(lái)說(shuō)，凡是涉及到指數(shù)運(yùn)算的代數(shù)公式，如(as)t
= a(st)，對(duì)四元數(shù)都不適用。
現(xiàn)在，我們已經(jīng)理解四元數(shù)求冪可以為我們做什么了。讓我們看看它的數(shù)學(xué)定義，四元數(shù)求冪定義在前一節(jié)討論的"有用"運(yùn)算上，定義如公式10.18：
注意，對(duì)于標(biāo)量求冪，也有類似結(jié)論：
不難理解為什么當(dāng)t從0變到1時(shí) q'從單位四元數(shù)變到 q
。注意到對(duì)數(shù)運(yùn)算只是提取了軸 n
和角度θ；接著，和指數(shù)t進(jìn)行標(biāo)量乘時(shí)，結(jié)果是θ乘以t；最后，指數(shù)運(yùn)算"撤銷"了對(duì)數(shù)運(yùn)算，從tθ和 n 重新計(jì)算w和
v 。上面給出的定義就是標(biāo)準(zhǔn)數(shù)學(xué)定義，在理論上非常完美，但直接轉(zhuǎn)換到代碼卻是很復(fù)雜的。程序清單10.1所示的代碼展示了怎樣計(jì)算
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;
}
關(guān)于這些代碼，需要注意的地方有:
（1）有必要做單位四元數(shù)的檢查。因?yàn)閣=+(-)1會(huì)導(dǎo)致mult的計(jì)算中出現(xiàn)除零現(xiàn)象。單位四元數(shù)的任意次方還是單位四元數(shù)。因此，如果檢測(cè)到輸入是單位四元數(shù)，忽略指數(shù)直接返回原四元數(shù)即可。
（2）計(jì)算alpha時(shí)，使用了acos函數(shù)，它的返回值是正的角度。這并不會(huì)違背一般性，任何四元數(shù)都能解釋成有正方向的旋轉(zhuǎn)角度，因?yàn)槔@某軸的負(fù)旋轉(zhuǎn)等價(jià)于繞指向相反方向的軸的正旋轉(zhuǎn)。
四元數(shù)插值 ---- "slerp"
當(dāng)今3D數(shù)學(xué)中四元數(shù)存在的理由是由于一種稱作slerp的運(yùn)算，它是球面線性插值的縮寫（Spherical Linear
Interpolation）。slerp運(yùn)算非常有用，因?yàn)樗梢栽趦蓚€(gè)四元數(shù)間平滑插值。slerp運(yùn)算避免了歐拉角插值的所有問(wèn)題。
slerp是一種三元運(yùn)算，這意味著它有三個(gè)操作數(shù)。前兩個(gè)操作數(shù)是兩個(gè)四元數(shù)，將在它們中間插值。設(shè)這兩個(gè)"開(kāi)始"和"結(jié)束"四元數(shù)分別為 q 0和
q 1。插值參數(shù)設(shè)為變量t，t在0到1之間變化，slerp函數(shù)：slerp( q 0,
q 1, t)，將返回 q 0到 q 1之間的插值方位。
能否利用現(xiàn)有的數(shù)學(xué)工具推導(dǎo)出slerp公式呢？如果是在兩個(gè)標(biāo)量a0和a1間插值，我們會(huì)使用下面的標(biāo)準(zhǔn)線性插值公式：
Δa = a1 - a0
lerp(a0, a1, t) = a0 + tΔa
標(biāo)準(zhǔn)線性插值公式從a0開(kāi)始，并加上a0和a1差的t倍，有三個(gè)基本步驟：
（1）計(jì)算兩個(gè)值的差。
（2）取得差的一部分。
（3）在初始值上加上差的一部分。
可以使用同樣的步驟在四元數(shù)間插值：
（1）計(jì)算兩個(gè)值的差， q 0到 q 1的角位移由Δ q
= q 0-1 q 1給出。
（2）計(jì)算差的一部分，四元數(shù)求冪可以做到，差的一部分由(Δ q )t給出。
（3）在開(kāi)始值上加上差的一部分，方法是用四元數(shù)乘法來(lái)組合角位移： q 0(Δ q )t
這樣，得到slerp的公式如公式10.19所示：
這是理論上的slerp計(jì)算過(guò)程，實(shí)踐中，將使用一種更加有效的方法。
我們?cè)?D空間中解釋四元數(shù)，因?yàn)樗形覀兏信d趣的四元數(shù)都是單位四元數(shù)，所以它們都"存在"于一個(gè) 4D"球面"上。
slerp的基本思想是沿著4D球面上連接兩個(gè)四元數(shù)的弧插值(這就是球面線性插值這個(gè)名稱的由來(lái))。
可以把這種思想表現(xiàn)在平面上，設(shè)兩個(gè)2D向量 v 0和 v 1，都是單位向量。我們要計(jì)算
v 1，它是沿 v 0到 v 1弧的平滑插值。設(shè)w是
v 0到 v 1弧所截的角，那么 v t就是繞
v 1沿弧旋轉(zhuǎn)tw的結(jié)果，如圖10.10所示：
將 v t表達(dá)成 v 0和 v 1的線性組合，從另一方面說(shuō)，存在兩個(gè)非零常數(shù)k0和k1，使得：
v t = k0 v 0
+ k1 v 1
可以用基本幾何學(xué)求出k0和k1，圖10.11展示了計(jì)算的方法：
對(duì)以k1 v 1為斜邊的直角三角形應(yīng)用三角公式得：
這里有兩點(diǎn)需要考慮。第一，四元數(shù) q 和- q
代表相同的方位，但它們作為slerp的參數(shù)時(shí)可能導(dǎo)致不一樣的結(jié)果，這是因?yàn)?D球面不是歐式空間的直接擴(kuò)展。而這種現(xiàn)象在2D和3D中不會(huì)發(fā)生。解決方法是選擇
q 0和 q 1的符號(hào)使得點(diǎn)乘 q 0 . q
1的結(jié)果是非負(fù)。第二個(gè)要考慮的是如果 q 0和 q
1非常接近，sinθ會(huì)非常小，這時(shí)除法可能會(huì)出現(xiàn)問(wèn)題。為了避免這樣的問(wèn)題，當(dāng)sinθ非常小時(shí)使用簡(jiǎn)單的線性插值。程序清單10.2把所有的建議都應(yīng)用到了計(jì)算四元數(shù)的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;