本文原理參考論文： GAMES101課程 ,基於Unity3D的Loop Subdivision 網格細分演算法

本文復現參考：Loop subdivision,三角網格下的Half-Edge資料結構實現方法
本文程式碼連結：https://pan.baidu.com/s/1R_Wn4iqBadGKxjFT8NgpsA 提取碼：9x9i

Loop 細分

顧名思義，網格細分是將粗糙網格精細化的過程，如下動圖。Loop細分是眾多網格細分演算法的一種，Loop細分僅僅對三角形網格模型有效。值得注意的是，雖然叫Loop細分，但是不能理解為“迴圈”細分，叫這個名字是因為作者的名字是這個而已。

演算法介紹

Loop細分演算法主要分兩步進行

1. 增加新頂點：增加三角形三個邊的中點，將一個三角形分成四個三角形
   

2. 移動新的三角形和老的三角形頂點讓細分之後的結果更加光滑

在介紹具體的演算法過程之前，先介紹一些概念：

邊界邊：處於三角形網格邊界的邊，只被網格中的一個三角形佔用。

內部邊：處於三角形網格內部的邊，被網格中的兩個三角形（最多隻有兩個）佔用。

邊界點：邊界邊的兩個頂點為邊界點。

內部點：除了邊界點的所有點為內部點。

下面一圖涵蓋所有定義：綠色框和紅色框分別表示邊界點和內部點，綠色標記和紅色標記分別表示邊界邊和內部邊

增加新頂點

通過在每個邊上取中點，從而增加頂點，根據邊的屬性（邊界，非邊界）有如下兩個計算新增頂點位置的規則：

1. 內部邊上增加點

計算規則：$v=3/8^{\ast }\left(v_0+v_1\right)+1/8^{\ast }\left(v_2+v_3\right)$

2. 邊界邊上增加點
   
   計算規則：$v=1/2^{\ast }\left(v_0+v_1\right)$

更新舊頂點

1. 內部點

   

   更新規則：$v=(1-n\beta )v_0+\beta {\sum }_{i=1}^n{v}_i$，n是頂點的度，$\beta$是一個與n有關的數，$\beta =\frac{1}{n}\left[\frac{5}{8}-{\left(\frac{3}{8}+\frac{1}{4}\cos \frac{2\pi }{n}\right)}^2\right]$

   這個公式也比較好理解， 更新一箇舊的頂點位置應該受與舊頂點相連的其他點的位置的影響以及舊的頂點自身位置的影響，通過對這兩個部分進行加權從而得到一箇舊頂點的位置更新。加權的含義在於，如果一個三角形連線和很多其他點，即度比較大，那麼該點原始位置就不太重要，佔有的比重應該比較小，而其他周圍的點影響該點的程度比較高，佔有的比重會比較大，反之亦然。

2. 邊界點

   影響邊界點的只有周圍的邊界點和自身，計算規則為：$\mathrm{v}=3/4^{\ast }{\mathrm{v}}_0+1/8^{\ast }\left({\mathrm{v}}_1+{\mathrm{v}}_2\right)$
   

演算法實現

本文演算法是通過半邊結構實現的，因為半邊結構能夠快速的得到我們期望的結果，如，獲得一個頂點所有相鄰的點這種操作在半邊結構中是非常迅速的，如果不瞭解半邊結構可以看一下本文的附加內容。

程式碼中TriMesh類是以半邊結構組織的網格模型的資料和拓撲資訊。

演算法實現流程：

1. 儲存更新舊頂點之後的位置
2. 用未更新舊頂點的模型計算出新增頂點的位置，通過掃描所有邊，插入新的頂點
3. 建立新增頂點和更新之後的舊頂點之間的拓撲關係，連線成網格，通過掃描所有面建立關係。

TriMesh* create_subDivisionMesh(TriMesh* oriMesh)
{
 	//建立一個新的TriMesh用於儲存細分之後的結果
    TriMesh* resultMesh = new TriMesh();
    auto &oriVertexs = oriMesh->Vertexs();
    //建立舊頂點到新頂點的對映（根據id，newVertexs與舊頂點一一對應）
    std::vector<Vertex* > newVertexs; 
    int vid = 0;
    //1. 儲存更新舊頂點之後的位置
    // 更新舊頂點位置，儲存在新的TriMesh和newVertexs陣列中
    for (int i = 0; i < oriVertexs.size(); i++)
    {
        vec3 newPosition(0,0,0);
        if(!oriVertexs[i]->isBoundary) //內部點處理
        {
            auto neighborVertexs = oriMesh->getNeighborVertexs(oriVertexs[i]);
            int n = neighborVertexs.size();//頂點的度
            float u = (5.0 / 8.0 - pow(3.0 / 8.0 + 1.0/4.0*std::cos(2 * PAI / n), 2)) / n;
            vec3 neighborPosition_sum(0,0,0);
            for (int j = 0; j < neighborVertexs.size(); j++)
            {
                auto neightbor = neighborVertexs[j];
                neighborPosition_sum += neightbor->vcoord;
            }
            newPosition = oriVertexs[i]->vcoord *(1 - n * u) + neighborPosition_sum*u;
        }
        else if (oriVertexs[i]->isBoundary) //邊界點處理
        {
            auto boundaryNeighborVert = oriMesh->getBoundaryNeighborVertexs(oriVertexs[i]);
            vec3 nvert_sum = boundaryNeighborVert[0]->vcoord + boundaryNeighborVert[1]->vcoord;
            newPosition = oriVertexs[i]->vcoord * 0.75 + nvert_sum * 0.125;
        }
        Vertex* newVertex = resultMesh->create_vertex(newPosition, vid++);
        newVertexs.push_back(newVertex);
    }
    //2. 用未更新舊頂點的模型計算出新增頂點的位置,通過掃描所有邊，插入新的頂點
    //建立邊到新插入點的對映
    std::unordered_map<EdgeKey, Vertex*, EdgeKeyHashFuc> map_ev; 
    auto oriEdges = oriMesh->HalfEges();
    int Count = 0;
    for (int i = 0; i < oriEdges.size(); i++)
    {
        HalfEdge* he = oriEdges[i];
        Vertex* v1 = he->next->next->vert;
        Vertex* v2 = he->vert;
        EdgeKey key(v1->id,v2->id);
        //半邊所在的邊已經建立了點
        if (map_ev.find(key) != map_ev.end())
        {
            continue;
        }
        if (he->isBoundary) //邊界邊插入點
        {
            vec3 newPos = (v1->vcoord + v2->vcoord) / 2.0;
            Vertex* newVertex = resultMesh->create_vertex(newPos, vid++);
            Count++;
            map_ev[EdgeKey(v1->id, v2->id)] = newVertex;
            map_ev[EdgeKey(v2->id, v1->id)] = newVertex;
        }
        else //內部邊插入點
        {
            Vertex* v3 = he->next->vert;
            Vertex* v4 = he->opposite->next->vert;
            vec3 newPos = (v1->vcoord + v2->vcoord) * (3.0 / 8.0) + (v3->vcoord + v4->vcoord) * (1.0 / 8.0);
            Vertex* newVertex = resultMesh->create_vertex(newPos, vid++);
            Count++;
            map_ev[EdgeKey(v1->id, v2->id)] = newVertex;
            map_ev[EdgeKey(v2->id, v1->id)] = newVertex;
        }
    }
    //3. 建立新增頂點和更新之後的舊頂點之間的拓撲關係，連線成網格
    // 使用舊的Trimesh的面中的頂點資訊建立面並儲存，新建立的面是原始TriMesh的四倍
    auto oriFaces = oriMesh->Faces();
    for (int i = 0; i < oriFaces.size(); i++)
    {
        Face* face = oriFaces[i];
        HalfEdge* fhe[3];
        fhe[0] = face->halfEdge;
        fhe[1] = fhe[0]->next;
        fhe[2] = fhe[1]->next;
        Vertex* center[3];
        // 根據線段找到剛剛新建立的點
        for (int i = 0; i < 3; i++) {
            auto key = getHalfEdgeKey(fhe[i]);
            center[i] = map_ev[key];
        }
        //建立以三個邊的中心點為中點的三角形面片
        resultMesh->create_face(center);
        Vertex* triVert[3];
        //建立另外三個三角形面片
        for (int i = 0; i < 3; i++)
        {
            int oriVertexId = fhe[i]->next->next->vert->id;
            triVert[0] = newVertexs[oriVertexId];
            triVert[1] = map_ev[getHalfEdgeKey(fhe[i])];
            triVert[2] = map_ev[getHalfEdgeKey(fhe[(i+2)%3])];
            resultMesh->create_face(triVert);
        }
    }
    //建立邊界標記
    resultMesh->creat_boundaryFlag();
    return resultMesh;
}

半邊結構

對三角形網格之進行搜尋是經常使用的，比如在平滑頂點法線的時候，需要搜尋到使用該頂點的所有面，然後對這些面法線進行加權平均。這樣的操作是非常繁瑣的，所以我們希望有一個資料結構能夠幫助我們快速的檢索三角面片資訊，半邊結構就是一個非常好的資料結構。

半邊結構的表示方法

半邊結構的主要特點在於將一條邊分成兩條方向相反的半邊，這就能夠很好的區分一條線段屬於那一條面，因為沒有一條半邊可以同時屬於兩個面。

半邊的資料結構如下：儲存了半邊指向的頂點，與半邊相反方向的對邊，下一條半邊，以及半邊所在面的資訊。

頂點的資料結構如下：儲存一個頂點相關資訊，以及從該點發出的半邊，雖然有很多個從該點發出的半邊，但是隻需要要記錄一個半邊即可，因為可以很快的檢索到其他相關的半邊。

**面的資料結構如下：**只需要儲存一個半邊，一個半邊可以不斷查詢下一個半邊以檢索到面中所有的頂點和半邊

半邊結構的構造方法

已知的資料有：已知頂點陣列 vector, 法線陣列vector可選，面陣列vector，這些都是可以從模型檔案(如obj檔案)中得到的。face中儲存三角形三個頂點的索引。

從已知資料獲得半邊結構的構造主要分兩步：建立所有頂點，建立所有面（同時建立邊）

建立所有頂點

建立頂點的過程中，只儲存頂點資料（位置，法線等），半邊資料設定為NULL，在之後建立面和半邊的過程中進行填充。讀取模型的過程可以得到三角網格的所有頂點的資訊，分別建立一個半邊結構中的Vertex物件，並儲存。

Vertex* create_vertex(vec3 point,int id)
{
    Vertex* vert = new Vertex(id,point);
    m_vertices.push_back(vert);
    return vert;
}

建立所有面及半邊

這一部分是半邊結構構造的核心過程，思想是然後遍歷每個三角形面，面中的每兩個頂點建立一個半邊。在建立半邊的時候可以對頂點的半邊資訊進行填充。

1. 建立面的過程如下，這一過程非常簡單，值得注意的是這個函式中僅能看到填充了半邊的下一個邊和麵的資訊。半邊的對邊資訊和頂點資訊在create_edge中進行填充。

   Face* create_face(Vertex* vertexs[3])
   {
       Face* face = new Face();
       HalfEdge* edges[3];
       //每兩個頂點構建一個邊
       for (int i = 0; i < 3; i++)
       {
           edges[i] = create_edge(vertexs[i % 3], vertexs[(i + 1) % 3]);
       }
       for (int i = 0; i < 3; i++)
       {
           edges[i]->next = edges[(i + 1) % 3];
           edges[i]->face = face;
           m_edges.push_back(edges[i]);
       }
       face->halfEdge = edges[0];
       m_faces.push_back(face);
       return face;
   }
   

2. 建立半邊：根據兩個頂點能夠構建出帶有頂點資料的半邊，半邊結構的精髓之處就是構建了半邊的對邊。如何構建半邊和半邊對邊之間的關係就非常重要了。具體的解決過程是：每當建立一個半邊的時候，都隨之建立一個該半邊的對邊。那麼，知道用於構建半邊的兩個頂點是否已經被用過了是一個重要的問題，如果這兩個頂點被用過了，說明這兩個頂點構建的兩個半邊（v1->v2,v2->v1）都已經被建立，就不用重新建立一次了。用於這種查詢可以使用雜湊表進行加速，以兩個頂點<v1,v2>作為查詢依據，返回的是構建的半邊。

HalfEdge* MeshLib::TriMesh::create_edge(Vertex* v1, Vertex* v2)
{
    if (v1 == NULL || v2 == NULL)
    {
        return NULL;
    }
    //首先查詢是否存在以這v1 v2建立的半邊
    EdgeKey key(v1->id,v2->id);
    if (m_hashmap_edge.find(key) != m_hashmap_edge.end()) //如果存在說明是之前建立好了的對邊
    {
        return m_hashmap_edge[key];
    }
    //不存在 則建立v1->v2的半邊以及其對邊( v2->v1的半邊),對邊儲存在雜湊表中以便於之後填充資料
    //建立時對邊的連線關係
    HalfEdge* he = new HalfEdge();
    HalfEdge* he_op = new HalfEdge();
    he->vert = v2;
    he->opposite = he_op;
    v1->halfEdge = he;
    he_op->vert = v1;
    he_op->opposite = he;
    //存入雜湊表
    m_hashmap_edge[EdgeKey(v1->id, v2->id)] = he;
    m_hashmap_edge[EdgeKey(v2->id, v1->id)] = he_op;
    return he;
}

半邊結構的使用

1. 獲得一個頂點發出的所有半邊

   //獲得一個頂點發出的所有半邊，半邊的上一個半邊的對邊和半邊的對邊的下一條邊是以該起點發出的半邊
   //如果沒有遇到邊界的邊則，只需要按照一個規則（預設搜尋半邊的對邊的下一條邊）不斷檢索下去
   //如果遇到半邊，則需要按照另一個規則（半邊的上一個半邊的對邊）檢索
   std::vector<HalfEdge*> getEdgesFromVertex(const Vertex* vertex)
   {
       std::vector<HalfEdge*> halfEdges;
       HalfEdge *he = vertex->halfEdge;
       HalfEdge *phe = he;
      do
      {
          if (phe->isBoundary) {
              halfEdges.push_back(phe);
              break;
          }
          halfEdges.push_back(phe);
          phe = phe->opposite->next;
      } while (phe!= he);
      if (phe->isBoundary)// 遇到邊界
      {
          phe = he->next->next;
          do
          {
              if (phe->isBoundary || phe==NULL) {     
                  break;
              }
              halfEdges.push_back(phe->opposite);
              phe = phe->opposite->next->next;
          } while (phe->vert->id== vertex->id);
      }
       return halfEdges;
   }
   

2. 獲得一個頂點相鄰的所有面

std::vector<Face*> MeshLib::TriMesh::getFacesFromVertex(const Vertex* vertex)
{
    std::vector<Face*> faces;
    auto edges = getEdgesFromVertex(vertex);
    for (int i = 0; i < edges.size(); i++) {
        faces.push_back(edges[i]->face);
    }
    return faces;
}

3. 獲得一個頂點的相鄰頂點

   std::vector<Vertex*> MeshLib::TriMesh::getNeighborVertexs(const Vertex* vertex)
   {
       std::vector<Vertex*> neighbors;
       auto faces = getFacesFromVertex(vertex);
       for (int i = 0; i < faces.size(); i++)
       {
           auto vertexs = getVertexsFromFace(faces[i]);
           for (int j = 0; j < vertexs.size(); j++)
           {
               bool isInsert = true;
               for( int t =0; t< neighbors.size(); t++)
                   if (vertexs[j]->id == neighbors[t]->id )
                   {
                       isInsert = false;
                   }
               isInsert = isInsert && (vertexs[j]->id != vertex->id);
               if(isInsert)
                   neighbors.push_back(vertexs[j]);
           }
       }
       return neighbors;
   }