從零開始手敲次世代遊戲引擎（四十八）

上一篇我們實現了一個由時間控制單變數形成的動畫。（雖然最終輸出的是一個4X4的矩陣，但是這個矩陣是由兩個獨立的旋轉：x軸旋轉、z軸旋轉合成的）

雖然對於平移和縮放，由於其線性特徵，我們可以將其按照坐標軸分解，使得每一步都是一個受時間控制的單變數過程（如：我們可以將「沿著向量{5, 3, 6}的平移」拆解為「沿x軸移動5，沿y軸移動3，沿z軸移動6」這3個步驟，且這3個步驟可以按照任何順序（包括同時進行）進行而不會影響最終結果），然而對於旋轉，我們則不能這麼處理了。（如：「繞著x軸旋轉90度，然後繞著z軸旋轉90度」與「繞著z軸旋轉90度，然後繞著x軸旋轉90度」所得到的結果是不一樣的，並且與「繞著x軸和z軸同時旋轉90度」也不是一個結果）

但是另外一方面，無論我們按照什麼順序如何進行旋轉，如果只看旋轉前的物體朝向（命名為v1）與旋轉之後的物體朝向（命名為v2），那麼我們都可以通過一次v1到v2的旋轉來得到同樣的結果。並且，這個v1到v2的旋轉可以看成是v1與v2所確定的平面（將v1與v2的起點重合，那麼這個起點和v1與v2的兩個終點一共3個點確定了一個平面）法線（可以通過v1叉乘v2得到）的旋轉。那麼如果記開始旋轉的時刻為t1，結束旋轉的時刻為t2，對於它們之間的任意時刻t

此時物體的朝向可以通過內插值v1與v2所形成的夾角得到。這個過程被稱為slerp過程。（參考引用1）

因此，對於平移矩陣和縮放矩陣，我們可以直接進行線性內插計算，或者使用從零開始手敲次世代遊戲引擎（四十七）所寫一元三次貝塞爾曲線進行內插計算，只不過需要將參數類型從浮點置換為向量類型就可以了。而對於旋轉矩陣，我們不能直接進行內插計算，而是需要將其轉換到球面坐標系（四元數）之後，再按照線性或者貝塞爾曲線方程來進行內插計算；或者，我們也可以計算出X、Y、Z方向的歐拉角，然後對這個角度進行線性內插計算，最後再用插補後的角度計算出旋轉矩陣。

下面的視頻給出了一個這樣的按照各個分量進行插值的動畫。

https://www.zhihu.com/video/965986560651460608

然而，還有一個更為麻煩的情況，那就是動畫參數是以一個4X4仿射矩陣（transform matrix）的形式存在，而不是按照縮放平移旋轉這樣的基本操作分別導出的。比如Blender當中的骨骼動畫，在導出為OGEX文件格式之後，就會發現它是以這種方式導出的。

首先，如果我們對於這個4X4矩陣直接進行內插值計算，雖然也可以得到動畫效果，但是這個動畫是會很奇怪的。動畫物體在兩個關鍵幀之間會發生很多形變和扭曲，整個動畫過程變得無法預料且不可控制。

這是因為一個4X4的仿射矩陣所能表達的幾何變化除了旋轉平移縮放這些之外，其實還包括歪斜（如同長方體變平行四邊形那樣）、透視、鏡面反射等情況。而且體現這些操作的矩陣當中的元素與旋轉和縮放是重疊的。

正確的做法是需要將這個矩陣進行分解，拆出平移、縮放、旋轉這些基本要素之後，再各自進行內插運算，再將它們合起來。然而，就如前面所說，在一個4X4的矩陣當中，左上方的3X3矩陣當中的元素包含了旋轉、縮放、歪斜、反射等多種情況，而且我們可以通過不同的步驟去改變一個物體，使其最終狀態是一致的。

這就是說，對於一個4X4矩陣的分解，理論上是可以有多數個可能性的。如果僅僅知道這個矩陣，要完全還原出各個基本操作，在理論上是不可能的。

但是對於動畫，一般來說我們只需要看起來合理就可以。所謂看起來合理，那就是符合人的經驗，至於是否完全還原製作時每一步的操作，其實是不重要，或者說不希望1:1還原的。

所以，我們這裡假定4X4矩陣是按照縮放（含反射）旋轉和平移這樣的步驟按順序完成的，並且每個步驟至多只有一次。也就是說，如果已知最終的變換矩陣M是由縮放矩陣S、旋轉矩陣R和平移矩陣T按順序相乘得到的，那麼是有辦法從M反推出S、R和T的。

這個過程稱為矩陣的分解（Matrix Decompose）。如上面所述，平移矩陣可以直接從4X4矩陣的第4行得到，所以接下來的就是分解4X4矩陣的左上角的3X3矩陣。在動畫當中常用的矩陣分解方法有QR分解法，SVD分解法和Polar分解法。

QR分解法（參考引用4）可以將任意一個方陣分解為一個純旋轉矩陣（Q）和一個上三角矩陣（R），但是問題在於分解出的Q並不是唯一的（也就是說不見得等同於我們原始的旋轉矩陣。事實上，當我們的縮放在各坐標軸上是同等的時候，QR分解所得矩陣就是我們的原始旋轉矩陣，然而當縮放不同等的時候，所得矩陣就不對了）。而且，QR分解所得的結果不能直接得到縮放矩陣。

SVD分解法(參考引用5)則可以將M分解為一個旋轉矩陣，一個按坐標軸方向的縮放矩陣，還有第二個旋轉矩陣。採用這個方法分解出的縮放矩陣雖然就是我們原始的縮放矩陣，但是旋轉部分被分解為兩個矩陣，且因為矩陣乘法不滿足交換律，我們不能簡單地將這兩個位於縮放矩陣一前一後的旋轉矩陣相乘來導出原始的旋轉矩陣。

Polar分解法（參考引用6）其實是SVD分解法的一個特例，然而當我們按照縮放->旋轉這樣的順序進行操作的時候，Polar分解法得到的兩個矩陣UP，U恰巧就是我們原始的旋轉矩陣，而P是一個對角矩陣，也就是我們原始的縮放矩陣。

（參考引用8）當中就對這三種分解法進行了理論上的比較，得出了Polar分解法最適合用在3D動畫矩陣插補當中的結論。

Polar分解法的理論基礎以及數值演算法基礎在（參考引用6）當中有詳細的說明，我這裡直接給出其演算法的實現：

template <typename T> inline void MatrixPolarDecompose(const Matrix<T, 3, 3>& in_matrix, Matrix<T, 3, 3>& U, Matrix<T, 3, 3>& P) { U = in_matrix; Matrix<T, 3, 3> U_inv; Matrix<T, 3, 3> U_pre; do { U_pre = U; U_inv = U; if (!InverseMatrix3X3f(U_inv)) assert(0); Matrix<T, 3, 3> U_inv_trans; Transpose(U_inv_trans, U_inv); U = (U + U_inv_trans) * (T)0.5; } while(U != U_pre); P = in_matrix * U_inv; }

演算法其實是從要分解的矩陣 作為起點，採用矩陣平方根的數值迭代演算法，計算出與 最為相似的單位陣 的平方根 。然後再將 乘以 ，得到P。

下面讓我們看看這個分解是否真的可以從 還原出

#include <iomanip>#include <iostream>#include <random>#include "MatrixComposeDecompose.hpp"using namespace My;using namespace std;int main(int , char** ){ default_random_engine generator; generator.seed(48); uniform_real_distribution<float> distribution_r(-1.0f * PI, 1.0f * PI); uniform_real_distribution<float> distribution_s(0.1f, 100.0f); uniform_real_distribution<float> distribution_t(-1000.0f, 1000.0f); auto dice_r = std::bind(distribution_r, generator); auto dice_s = std::bind(distribution_s, generator); auto dice_t = std::bind(distribution_t, generator); Vector3f translation ({dice_t(), dice_t(), dice_t()}); Vector3f scale ({dice_s(), dice_s(), dice_s()}); Vector3f rotation ({dice_r(), dice_r(), dice_r()}); Matrix4X4f matrix; Matrix4X4fCompose(matrix, rotation, scale, translation); cerr << "Scale: " << scale; Matrix3X3f A = {{ {matrix[0][0], matrix[0][1], matrix[0][2]}, {matrix[1][0], matrix[1][1], matrix[1][2]}, {matrix[2][0], matrix[2][1], matrix[2][2]} }}; Matrix3X3f U, P; MatrixPolarDecompose(A, U, P); cout.precision(4); cout.setf(ios::fixed); cout << "Polar Decompose of matrix A: " << endl; cout << A; cout << "U:" << endl; cout << U; cout << "Orthogonal: " << (U.isOrthogonal()?"true":"false") << endl; cout << "P:" << endl; cout << P; Matrix3X3f A_dash = P * U; cout << "U * P: " << A_dash; cout << "Error: " << A_dash - A; return 0;}

上面這個例子首先用偽隨機數生成了一個Transform矩陣，然後用Polar分解法對它進行分解。下面是執行的結果：

chenwenlideMBP:GameEngineFromScratch chenwenli$ build/Test/PolarDecomposeTestScale: ( 0.2078,8.2476,86.5060 )Polar Decompose of matrix A:( 0.1197,-0.1359,0.1019 )( 6.1709,5.4718,0.0487 )( -28.4736,31.4404,75.3926 )U:( 0.5761,-0.6540,0.4903 )( 0.7482,0.6634,0.0059 )( -0.3292,0.3634,0.8715 )Orthogonal: trueP:( 0.2078,0.0000,-0.0000 )( 0.0000,8.2476,0.0000 )( 0.0000,0.0000,86.5060 )U * P:( 0.1197,-0.1359,0.1019 )( 6.1709,5.4718,0.0487 )( -28.4736,31.4404,75.3926 )Error:( -0.0000,0.0000,-0.0000 )( 0.0000,0.0000,0.0000 )( 0.0000,0.0000,-0.0000 )chenwenlideMBP:GameEngineFromScratch chenwenli$

我們可以看到分解出的P是一個對角矩陣，其中的值正好與我們隨機產生的原始矩陣的縮放值相同。而U是旋轉矩陣，它的值其實與我們隨機產生的旋轉所對應的旋轉矩陣是一樣的。（這個例子里看不太出來）。關於這一點，我們可以在下面這個完整的對變換矩陣進行線性插值的例子當中驗證：

#include <iostream>#include <random>#include "Linear.hpp"#include "geommath.hpp"using namespace My;using namespace std;int main (int argc, char** argv){ int interpolate_count = 100; if (argc > 1) { interpolate_count = atoi(argv[1]); } default_random_engine generator; generator.seed(48); uniform_real_distribution<float> distribution_r(-1.0f * PI, 1.0f * PI); uniform_real_distribution<float> distribution_s(0.1f, 100.0f); uniform_real_distribution<float> distribution_t(-1000.0f, 1000.0f); auto dice_r = std::bind(distribution_r, generator); auto dice_s = std::bind(distribution_s, generator); auto dice_t = std::bind(distribution_t, generator); // generate start point matrix Vector3f translation_1 ({dice_t(), dice_t(), dice_t()}); Vector3f scale_1 ({dice_s(), dice_s(), dice_s()}); Vector3f rotation_1 ({dice_r(), dice_r(), dice_r()}); Matrix4X4f matrix_transform_1; Matrix4X4fCompose(matrix_transform_1, rotation_1, scale_1, translation_1); Vector3f v ({dice_t(), dice_t(), dice_t()}); // generate end point matrix Vector3f translation_2 ({dice_t(), dice_t(), dice_t()}); Vector3f scale_2 ({dice_s(), dice_s(), dice_s()}); Vector3f rotation_2 ({dice_r(), dice_r(), dice_r()}); Matrix4X4f matrix_transform_2; Matrix4X4fCompose(matrix_transform_2, rotation_2, scale_2, translation_2); auto v1 = v; TransformCoord(v1, matrix_transform_1); auto v2 = v; TransformCoord(v2, matrix_transform_2); cout << "Start Point:" << v1; cout << "End Point:" << v2; cout << "_________________" << endl; cout << "Start Translation: " << translation_1; cout << "End Translation: " << translation_2; cout << "_________________" << endl; cout << "Start Scalar: " << scale_1; cout << "End Scalar: " << scale_2; cout << "_________________" << endl; cout << "Start Rotation: " << rotation_1; cout << "End Rotation: " << rotation_2; cout << "=================" << endl; Linear<Matrix4X4f, float> linear_introplator({matrix_transform_1, matrix_transform_2}); cout << endl; cout << "Interpolate: " << endl; for (int i = 0; i <= interpolate_count; i++) { auto inter_matrix = linear_introplator.Interpolate(i * 1.0f / interpolate_count, 1); Vector3f rotation, scalar, translation; Matrix4X4fDecompose(inter_matrix, rotation, scalar, translation); auto v_inter = v; TransformCoord(v_inter, inter_matrix); cout << "#" << i << endl; cout << "_________________" << endl; cout << "Interpolated Position: " << v_inter; cout << "_________________" << endl; cout << "Rotation: " << rotation; cout << "Scalar: " << scalar; cout << "Translation: " << translation; cout << "=================" << endl; } return 0;}

上面這個例子隨機產生兩個變換矩陣和一個向量，並將向量分別乘以這兩個矩陣得到其初始位置（相當於動畫開始時的Pose）和最終位置（相當於動畫結束時的Pose）。然後我們對這兩個矩陣進行線性內插值，以便計算出向量在開始位置到終止位置的運動軌跡。我們可以在命令行輸入具體插值的點數（預設為插100個點）。為了方便起見，我們下面展現的是插入3個點的結果：

chenwenlideMBP:GameEngineFromScratch chenwenli$ build/Test/LinearInterpolateTest 3Start Point:( 21090.7422,-26166.0117,-58756.9805 )End Point:( -30576.4004,-36326.3008,10601.1221 )_________________Start Translation: ( -997.8421,-836.8849,729.8501 )End Translation: ( 437.6792,-791.1947,238.9918 )_________________Start Scalar: ( 0.2078,8.2476,86.5060 )End Scalar: ( 79.7544,46.2412,10.6004 )_________________Start Rotation: ( -3.1348,-2.6292,2.2929 )End Rotation: ( 1.8683,-0.2396,-2.4812 )=================Interpolate:#0_________________Interpolated Position: ( 21090.7422,-26166.0117,-58756.9766 )_________________Rotation: ( 0.0068,-0.5124,-0.8487 )Scalar: ( 0.2078,8.2476,86.5060 )Translation: ( -997.8421,-836.8849,729.8501 )=================#1_________________Interpolated Position: ( 31601.7500,-26435.5840,-29497.2344 )_________________Rotation: ( 0.6273,-0.4215,-1.3929 )Scalar: ( 26.7233,20.9121,61.2041 )Translation: ( -519.3350,-821.6548,566.2307 )=================#2_________________Interpolated Position: ( 12253.8975,-41567.6953,-133.1098 )_________________Rotation: ( 1.2478,-0.3305,-1.9370 )Scalar: ( 53.2389,33.5767,35.9023 )Translation: ( -40.8279,-806.4247,402.6112 )=================#3_________________Interpolated Position: ( -30576.3984,-36326.2969,10601.1201 )_________________Rotation: ( 1.8683,-0.2396,-2.4812 )Scalar: ( 79.7544,46.2412,10.6004 )Translation: ( 437.6792,-791.1947,238.9918 )=================chenwenlideMBP:GameEngineFromScratch chenwenli$

可以看到，我們精確地從兩個隨機生成的矩陣得到了原始的動畫參數，並且無論是向量的位置、中間變換矩陣的縮放、旋轉、位移都實現了線性的變化。

最終代碼在article_48當中。

參考引用