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

04-03

在從零開始手敲次世代遊戲引擎（四十三）當中，我們根據場景物體的頂點數據生成了包裹場景物體的AABB碰撞盒。這個碰撞盒可以用於我們在從零開始手敲次世代遊戲引擎（四十二）介紹的BroadPhase。

同樣，在從零開始手敲次世代遊戲引擎（四十二）當中我們介紹了NarrowPhase，它使用凸包取代AABB碰撞盒，來獲得更為精確的碰撞檢測結果。接下來就讓我們看看如何通過場景物體的頂點數據來生成這個凸包。

從頂點數據生成凸包的研究在計算幾何學（圖形學）當中已經有至少幾十年的歷史，很多演算法被提出來。參考引用1當中列出了一些著名的演算法及其複雜度。這裡面Graham Scan是一個複雜度為 $O(n log n)$ 的比較經典的演算法（參考引用2，參考引用4），後面的複雜度為 $O(n log h)$ 的Chans演算法（參考引用3）也是基於它的。

Graham Scan演算法。圖片來源：Wikipedia

據Wikipedia，Graham Scan命名自美國數學家Ronald Graham，供職於加州大學伯克利分校。他老婆叫金芳蓉，出生於台灣，也是一位數學家，供職在賓夕法尼亞大學。（雖然這與我們要講的沒啥關係。。。）

Graham Scan演算法由3個步驟組成：

找到起點。這個點必須保證是在最終的凸包上；
從起點開始，對點集進行某種遞增的排序，以獲得一個比較初始的外殼檢索順序；
從起點開始按照第二步確定的檢索順序依次判定點是否在凸包上。

然而這個演算法我沒有能夠找到將其拓展到3D空間的相關的資料。從個人感覺上來說，因為這個演算法在2D空間是根據 $v_{n-2},v_{n-1},v_{n}$ 之間的折線「左拐」還是「右拐」來推斷 $v_{n-1}$ 是否位於凸包上的，那麼將其擴展到3D空間的話，應該是需要至少考察 $v_{n-3},v_{n-2},v_{n-1},v_{n}$ 這4個點，看看由 $v_{n-3},v_{n-2},v_{n-1}$ 這3個點所確定的平面與 $v_{n-2},v_{n-1},v_{n}$ 這3個點所確定的平面之間的夾角是「左拐」還是「右拐」。但是這裡面有幾個問題：

如何對3D的點集進行排序，使其滿足Graham Scan的要求？
當判定為「右拐」的時候，是應該捨棄 $v_{n-2}$ 還是 $v_{n-1}$ ？

等等。這些看來需要比較深入的思考。在參考引用4當中，也沒有Graham Scan拓展到3D空間的例子，只有Gift Wrapping拓展到3D空間的例子，看來這個拓展應該不是很容易，或者不存在。

因為我們的主要對象是3D場景，所以我這裡選用已知可以拓展到3D場景的Quickhall演算法。Quickhall演算法的工作原理如下圖所示：

2D Quickhull演算法工作流程演示。圖片來源：Wikipedia

3D空間當中的Quickhull演算法的工作原理與2D空間當中相同，但是基本型由三角形變成了四面體，因此稍微複雜一些：

我們首先分別在 $x,y,z$ 三個坐標軸方向上找到最大最小點，共3對6個點（這6個點必在凸包上）；
在這6個點當中，選取距離最大的兩個點為初始點；
在這6個點當中，選取距離第二步所選2個點連線最遠的點作為第3個點；
遍歷點集當中所有點，找到距離這3個點所確定的平面最遠的點作為第4個點；
這4個點組成了一個四面體，將所有位於四面體內的點排除在搜索對象外；
將所有位於4面體外的點按照四面體的四個面進行分組（分組的依據是這個點位於4面體的哪個面所在平面的前方，也就是這個點能夠看到四面體的哪個面。有一些位置上的點是可以同時看到四面體的至多3個面的。對於這種情況，將其放在可見的任意一個面的分組裡）
檢查每個面分組裡的點的個數。如果為零，這個面無需繼續處理。如果大於零，那麼這個面需要繼續處理，將該面分組壓入棧；
從棧中彈出面分組，遍歷分組當中的點，找到離開面距離最遠的點，然後從此點觀察棧中所有的面，剔除這些面，並將所形成的缺口上的點與本步驟當中找到的新點連線，形成一組新的面。將被剔除的那些面的面分組當中的點重新分配到新的面，捨棄位於所有新面背後的點。然後回到步驟7；
當所有面分組裡面都不再有新點的時候，我們就已經找到了點集所對應的凸包。

下面就是實際執行的結果。我們先在三維坐標系的單位空間當中隨機布置了100個點，然後用上面的演算法計算其凸包。為了展示這個計算過程，我們每次執行完第8步的時候就輸出圖形，並且等待按下鍵盤上的U鍵才回到第7步繼續迭代。也就是說，每按下一次「U」鍵，迭代一次，得到一個更為精確的凸包。事實上，在商業的遊戲物理引擎當中，也是採用這個迭代的方法，將凸包的計算量分散到多次的物理引擎基本循環當中，來避免對CPU在某個時間點產生大的負荷。況且，對於還比較遠的物體，我們只需要它的一個大概的包圍盒即可以，隨著物體的接近，我們再逐步迭代，完善其包圍盒。

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

因為本篇所用畫面是通過我們之前從文章從零開始手敲次世代遊戲引擎（四十一）到文章從零開始手敲次世代遊戲引擎（四十三）所寫的調試用繪圖介面繪製的，實際上並沒有場景物體，所以應用程序剛剛啟動的時候什麼都沒有，需要按下鍵盤上的「D」按鍵打開調試用圖形才能看到。

有了這個演算法之後，我們就可以對任意的三維模型生成它的凸包，也就是用於物理碰撞檢測的NarrowPhase的包圍盒。

下面是演算法的核心代碼。完整的代碼在netwarm007/GameEngineFromScratch當中。這個實現還有很多可以改進效率的地方，比如在將點進行分組的地方，現在是每次迭代都對所有剩下的點進行重新分組。然而事實上，我們只需要對那些在上一個迭代當中變化了的面的分組當中的點進行重新分組就可以了。

另外，我在寫這個代碼的時候其實前前後後零零碎碎也花了大約1周的時間。遇到的主要問題有下面這兩個：

用for( : )對std::set進行迭代，並在循環體當中對set的元素進行erase()操作導致for(:)永遠不終止或者拋出異常的情況。這裡需要特別注意的時候erase()會將傳遞給他的參數invalidate，當這個參數正好也是循環變數的時候就會導致這個問題。解決的方法是不要使用for(:)而是使用通常的迭代器（Iterator），然後用erase()的返回值更新這個iterator。由於erase()返回的iterator已經指向了被刪除元素後面的一個元素，所以使用for（auto it = begin(); it != end(); it++)這種形式也是不行的，會導致少循環一次；
為了減少內存消耗，減少複製從而提高性能，我在Point、Edge、Face等結構裡面存儲的都是指向圖元的指針，比如PointPtr、EdgePtr、FacePtr。比如一個Edge當中就包括了起點和終點的PointPtr。由於在編碼的時候就保證了每個圖元只有一個副本，加上所使用的std::shared_ptr重載了比較運算符，使得指向同一個圖元的不同的std::shared_ptr實例能夠被正確判定為相等，因此大多數情況當中沒有問題。然而，在上述步驟的第8步當中，對於移除表面所形成的空洞的邊緣（Edge）的計算，因為同一個Edge在不同的Face裡面方向不同（因為我們用的是右手坐標系，所以Edge必須按照逆時針方向組織到Face裡面，這個Face才是面向外面的），就導致同一條邊在不同Face裡面實際上是不同的實例。所以需要重載Edge的比較運算符，保證這些被不同Face共享的Edge能被正確識別出來。

void QuickHull::ComputeInitialTetrahydron(){ if(m_PointSet.size() < 4) { // too few points in the point set, nothing could be done return; } PointPtr ExtremePointXMin = make_shared<Point>(numeric_limits<float>::max()); PointPtr ExtremePointYMin = make_shared<Point>(numeric_limits<float>::max()); PointPtr ExtremePointZMin = make_shared<Point>(numeric_limits<float>::max()); PointPtr ExtremePointXMax = make_shared<Point>(numeric_limits<float>::min()); PointPtr ExtremePointYMax = make_shared<Point>(numeric_limits<float>::min()); PointPtr ExtremePointZMax = make_shared<Point>(numeric_limits<float>::min()); // copy the point set into temporary working set m_PointWaitProcess = m_PointSet; // finding the Extreme Points [O(n) complexity] for(auto point_ptr : m_PointWaitProcess) { if(point_ptr->x < ExtremePointXMin->x) ExtremePointXMin = point_ptr; if(point_ptr->y < ExtremePointYMin->y) ExtremePointYMin = point_ptr; if(point_ptr->z < ExtremePointZMin->z) ExtremePointZMin = point_ptr; if(point_ptr->x > ExtremePointXMax->x) ExtremePointXMax = point_ptr; if(point_ptr->y > ExtremePointYMax->y) ExtremePointYMax = point_ptr; if(point_ptr->z > ExtremePointZMax->z) ExtremePointZMax = point_ptr; } PointSet ExtremePoints; ExtremePoints.insert({ ExtremePointXMin, ExtremePointXMax, ExtremePointYMin, ExtremePointYMax, ExtremePointZMin, ExtremePointZMax }); // now we find the most distant pair among 6 extreme points { float distance_x, distance_y, distance_z; distance_x = Length(*ExtremePointXMax - *ExtremePointXMin); distance_y = Length(*ExtremePointYMax - *ExtremePointYMin); distance_z = Length(*ExtremePointZMax - *ExtremePointZMin); int max_distance_index = (distance_x < distance_y)? ( (distance_y < distance_z)?2:1 ) : ( (distance_x < distance_z)?2:0 ); PointPtr A, B; switch(max_distance_index) { case 0: A = ExtremePointXMin; B = ExtremePointXMax; break; case 1: A = ExtremePointYMin; B = ExtremePointYMax; break; case 2: A = ExtremePointZMin; B = ExtremePointZMax; break; default: assert(0); } ExtremePoints.erase(A); ExtremePoints.erase(B); // now we have 2 points on the convex hull, find the 3rd one among remain // extreme points PointPtr C; { double max_distance = 0.0f; for (auto point_ptr : ExtremePoints) { Vector3f numerator; Vector3f denominator; CrossProduct(numerator, *point_ptr - *A, *point_ptr - *B); denominator = *B - *A; double distance = Length(numerator) / Length(denominator); if (distance > max_distance) { C = point_ptr; max_distance = distance; } } } m_PointWaitProcess.erase(A); m_PointWaitProcess.erase(B); m_PointWaitProcess.erase(C); // now we find the 4th point to form a tetrahydron PointPtr D; { float max_distance = 0; for (auto point_ptr : m_PointWaitProcess) { auto distance = PointToPlaneDistance({A,B,C}, point_ptr); if (distance > max_distance) { D = point_ptr; max_distance = distance; } } } center_of_tetrahydron = make_shared<Point>((*A + *B + *C + *D) * 0.25f); m_ConvexHull.AddTetrahydron({A,B,C,D}); m_PointWaitProcess.erase(D); }}void QuickHull::IterateHull() { AssignPointsToFaces(); cerr << "remain point count: " << m_PointWaitProcess.size() << endl; auto pPoint = *m_PointWaitProcess.begin(); auto pFace = m_PointAboveWhichFacies.find(pPoint)->second; float max_distance = 0.0f; auto vertices = pFace->GetVertices(); auto range = m_PointsAboveFace.equal_range(pFace); PointPtr far_point = nullptr; for_each( range.first, range.second, [&](decltype(m_PointsAboveFace)::value_type x) { auto distance = PointToPlaneDistance(vertices, x.second); if (distance > max_distance) { far_point = x.second; max_distance = distance; } } ); if (far_point) { // remove all faces this point can see and // create new faces by connecting all vertices // on the border of hole to the new point set<Edge> edges_on_hole; auto range = m_PointAboveWhichFacies.equal_range(far_point); for_each( range.first, range.second, [&](decltype(m_PointAboveWhichFacies)::value_type x) { auto face_to_be_removed = x.second; for (auto edge : face_to_be_removed->Edges) { if (edges_on_hole.find(*edge) != edges_on_hole.end()) { // this edge is shared by faces going to be removed // so it is not on the border of hole, remove it edges_on_hole.erase(*edge); } else { // temporary add it edges_on_hole.insert(*edge); } } m_ConvexHull.Faces.erase(face_to_be_removed); } ); // now we have edges on the hole // so we create new faces by connecting // them with the new point //assert(edges_on_hole.size() >= 3); for (auto edge : edges_on_hole) { m_ConvexHull.AddFace({edge.first, edge.second, far_point}, center_of_tetrahydron); } // now the point has been proceeded // so remove it from the waiting queue m_PointWaitProcess.erase(far_point); }}void QuickHull::AssignPointsToFaces(){ m_PointsAboveFace.clear(); m_PointAboveWhichFacies.clear(); auto it = m_PointWaitProcess.begin(); while (it != m_PointWaitProcess.end()) { bool isInsideHull = true; for (auto pFace : m_ConvexHull.Faces) { if (isPointAbovePlane(pFace->GetVertices(), *it)) { m_PointsAboveFace.insert({pFace, *it}); // record all faces // the point can "see" in order to extrude the // convex hull face m_PointAboveWhichFacies.insert({*it, pFace}); isInsideHull = false; } } if (isInsideHull) { // return an iterator to the element that follows the last element removed // (or set::end, if the last element was removed). it = m_PointWaitProcess.erase(it); } else { it++; } }}