【图形学】Loop细分算法及半边结构实现(C++)

本文原理参考： GAMES101课程 ,基于Unity3D的Loop Subdivision 网格细分算法

本文复现参考：Loop subdivision,三角网格下的Half-Edge数据结构实现方法
本文代码链接：https://pan.baidu.com/s/1R_Wn4iqBadGKxjFT8NgpsA 提取码：9x9i

文章目录

• Loop 细分
• • 算法介绍
  • • 增加新顶点
    • 更新旧顶点
  • 算法实现
• 半边结构
• • 半边结构的表示方法
  • 半边结构的构造方法
  • • 创建所有顶点
    • 创建所有面及半边
  • 半边结构的使用

Loop 细分

顾名思义，网格细分是将粗糙网格精细化的过程，如下动图。Loop细分是众多网格细分算法的一种，Loop细分仅仅对三角形网格模型有效。值得注意的是，虽然叫Loop细分，但是不能理解为“循环”细分，叫这个名字是因为作者的名字是这个而已。

算法介绍

Loop细分算法主要分两步进行

1. 增加新顶点：增加三角形三个边的中点，将一个三角形分成四个三角形
   

2. 移动新的三角形和老的三角形顶点让细分之后的结果更加光滑

在介绍具体的算法过程之前，先介绍一些概念：

边界边：处于三角形网格边界的边，只被网格中的一个三角形占用。

内部边：处于三角形网格内部的边，被网格中的两个三角形（最多只有两个）占用。

边界点：边界边的两个顶点为边界点。

内部点：除了边界点的所有点为内部点。

下面一图涵盖所有定义：绿色框和红色框分别表示边界点和内部点，绿色标记和红色标记分别表示边界边和内部边

增加新顶点

通过在每个边上取中点，从而增加顶点，根据边的属性（边界，非边界）有如下两个计算新增顶点位置的规则：

1. 内部边上增加点

计算规则： v = 3 / 8 ∗ ( v 0 + v 1 ) + 1 / 8 ∗ ( v 2 + v 3 ) v=3 / 8^{*}\left(v_{0}+v_{1}\right)+1 / 8^{*}\left(v_{2}+v_{3}\right) v=3/8∗(v0​+v1​)+1/8∗(v2​+v3​)

2. 边界边上增加点
   
   计算规则： v = 1 / 2 ∗ ( v 0 + v 1 ) v=1 / 2^{*}\left(v_{0}+v_{1}\right) v=1/2∗(v0​+v1​)

更新旧顶点

1. 内部点

   

   更新规则： v = ( 1 − n β ) v 0 + β ∑ i = 1 n v i v=(1-n \beta) v_{0}+\beta \sum_{i=1}^{n} v_{i} v=(1−nβ)v0​+β∑i=1n​vi​，n是顶点的度，$\beta$是一个与n有关的数， β = 1 n [ 5 8 − ( 3 8 + 1 4 cos ⁡ 2 π n ) 2 ] \beta=\frac{1}{n}\left[\frac{5}{8}-\left(\frac{3}{8}+\frac{1}{4} \cos \frac{2 \pi}{n}\right)^{2}\right] β=n1​[85​−(83​+41​cosn2π​)2]

   这个公式也比较好理解， 更新一个旧的顶点位置应该受与旧顶点相连的其他点的位置的影响以及旧的顶点自身位置的影响，通过对这两个部分进行加权从而得到一个旧顶点的位置更新。加权的含义在于，如果一个三角形连接和很多其他点，即度比较大，那么该点原始位置就不太重要，占有的比重应该比较小，而其他周围的点影响该点的程度比较高，占有的比重会比较大，反之亦然。

2. 边界点

   影响边界点的只有周围的边界点和自身，计算规则为： v = 3 / 4 ∗ v 0 + 1 / 8 ∗ ( v 1 + v 2 ) \mathrm{v}=3 / 4^{*} \mathrm{v}_{0}+1 / 8^{*}\left(\mathrm{v}_{1}+\mathrm{v}_{2}\right) v=3/4∗v0​+1/8∗(v1​+v2​)
   

算法实现

本文算法是通过半边结构实现的，因为半边结构能够快速的得到我们期望的结果，如，获得一个顶点所有相邻的点这种操作在半边结构中是非常迅速的，如果不了解半边结构可以看一下本文的附加内容。

代码中TriMesh类是以半边结构组织的网格模型的数据和拓扑信息。

算法实现流程：

1. 保存更新旧顶点之后的位置
2. 用未更新旧顶点的模型计算出新增顶点的位置，通过扫描所有边，插入新的顶点
3. 建立新增顶点和更新之后的旧顶点之间的拓扑关系，连接成网格，通过扫描所有面建立关系。

TriMesh* create_subDivisionMesh(TriMesh* oriMesh)
{//创建一个新的TriMesh用于存储细分之后的结果TriMesh* resultMesh = new TriMesh();auto &oriVertexs = oriMesh->Vertexs();//建立旧顶点到新顶点的映射（根据id，newVertexs与旧顶点一一对应）std::vector 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 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）都已经被创建，就不用重新创建一次了。用于这种查询可以使用哈希表进行加速，以两个顶点作为查询依据，返回的是构建的半边。

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

3. 获得一个顶点的相邻顶点

   std::vector MeshLib::TriMesh::getNeighborVertexs(const Vertex* vertex)
   {std::vector 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;
   }
   

本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！