标签：
概述
首先简要介绍一下应用场景：给定一个空间的Box范围和其中一系列空间曲面(使用三角网格Mesh表达的)，要对这个范围做一个四面体剖分(Tetrahedralization)，四面体们按一个合理的密度组成Box，并且曲面上的三角片一定是最后其中某些四面体的一个面。其效果图类似于下图左图到右图：
 |
 |
| box+给定曲面 | box的四面体剖分 |
用通俗的说法再解释一编就是，如何实现将这个box分成无数小三棱锥或者说四面体，并且要求四面体的面严格包含给定的三角片。这样的剖分结果有什么用处目前不是完全明白，目前知道的用途是做有限元分析，或者是做形态变化的演化等等，总之具体的用处可以不必关心。针对这个场景下的问题，目前通过调查研究，最后问题锁定到一个研究邻域上，叫做限定四面体剖分(Constrained Tetrahedralization)，特此写本文展示一步一步解决这个问题的过程，描述其中遇到的问题和最终的解决方案。
概念：三角剖分/Delaunay三角剖分/四面体剖分
为了研究本文的问题，从最原始最简单的一个基本问题点开始，那就是是三角剖分(Triangulation)，三角剖分有自己形式化的定义，但本文旨在形象的介绍这些概念。三角剖分简而言之就是将二维平面中的点集连成多个三角形，三角形之间没有重叠部分。对于一个固定的点集，三角剖分是不唯一的。三角剖分中最著名的算法是Delaunay三角剖分，这个算法简单理解可以认为输出一种最优的剖分，对给定的点集，Delaunay三角剖分的结果是唯一的，其因为具有很多优良特性因此被广泛使用。关于Delaunay三角剖分，本文不详细讲述，网上有很多资料可以查询。下图是一个Delaunay三角剖分示意：
四面体剖分是三角剖分扩展到高维度的情况，输入是点集，输出是空间四面体剖分，即空间区域以四面体为基本划分单位。下图是一个四面体剖分。
如何表示三角剖分结果/Mesh数据结构
虽然Mesh更多是用作空间形状的表示，但也可以用于表示三角剖分的结果，表示方法仍然是点集+三角片集合。关于Mesh的基础可以参考文章—"三角网格数据结构"和"三角网格数据结构2"。
如何表示四面体剖分结果/Tetrahedral Mesh的设计与实现
三角网格采用点集+三角片集合的形式来表示，推而广之，四面体剖分由于基本单元变成了四面体，故采用点集+四面体集合的方式表示，其中四面体采用4个点索引值表示，其定义如下：
struct Tetrahedra { int PIndex[4]; Tetrahedra(int p0,int p1,int p2,int p3) { PIndex[0]=p0; PIndex[1]=p1; PIndex[2]=p2; PIndex[3]=p3; } Tetrahedra() { PIndex[0]=-1; PIndex[1]=-1; PIndex[2]=-1; PIndex[3]=-1; } };
由于四面体由三角片构成，故设计一个三角片集合用于存储所有不重合三角片是有必要的。类似于三角网，四面体网格也存在各种几何图元的拓扑关系，这些关系在后续处理中也很重要，故给这些关系也设计相应的存储结构。新的四面体网格类型TetraMesh类可以继承自Mesh类，设计如下，其中相应的一些基本结构如Point等的定义可在"三角网格数据结构"和"三角网格数据结构2"里获知：
class Mesh { public: std::vector<Point3d> Vertices; std::vector<Triangle> Faces; std::vector<Vector> VertexNormals; public: Mesh(); ~Mesh(); MESH_INDEXTYPE AddVertex(Point3d& toAdd); MESH_INDEXTYPE AddFace(Triangle& triangle); }; class TetraMesh : public Mesh { public: std::vector<Triangle> InnerTriangles;//所有剖分的三角片，Triangles成员默认存外表面的三角片 std::vector<Tetrahedra> Tedtrahedras;//所有四面体 std::vector<Edge> Edges;//所有边，暂时用不上 std::vector<std::vector<int> > NeighborTetraPerTriangle;//每个三角片的邻接四面体集合 std::vector<std::vector<int> > NeighborTetraPerVertex;//每个点的邻接四面体集合 std::vector<std::vector<int> > NeighborTrianglePerEdge;//每个边的邻接三角片集合，用不上 std::vector<std::vector<int> > NeighborVertexPerVertex;//每个点的邻接点集合 public: MESH_INDEXTYPE AddEdge(Edge& e); MESH_INDEXTYPE AddTetra(Tetrahedra& tet); MESH_INDEXTYPE AddInnerTriangle(Triangle& t); void InitNeighborTetraPerVertex();//初始化VT关系 void InitInnerTriangles();//初始化三角片和FT关系 void InitSurfaceEdges();//初始化边和EF关系 void InitNeighborVertexPerVertex();//初始化VV关系 };
其中点集和四面体集合是必要信息，其余的拓扑信息可以从这两个信息中计算而来，下一节重点解决这一问题。
在Tetrahedral Mesh中求取相关邻接数据信息
首先在几何处理领域，以这几个简单的字母来指代相应几何元素，或者几何元素的集合：
然后用这些字母组合来指代在三角网、三角剖分和四面体剖分这样的数据中几何元素之间的关系，例如点与点(VT)，点与边(VE)等，其中在本文的范畴内需要使用的关系描述如下：
 |
 |
 |
 |
如何获得点体关系(VT)在本文具体实现中，这些都是关系采用vector<vector<int>>的方式实现的，例如VF关系是设置一个与V数组等大的数组每个元素对应相应的V的邻接面的索引，这样这个数组的每个元素也是数组，即vector<int>。
最基本的四面体数据结构只包含两个集合：点集V和四面体集合T，那么如何进一步从这两个集合推导出其他需要的信息呢？下面罗列出一些具体问题的解决方法：
1) 如何获取VT关系？
这个比较简单，参考博文"三角网格数据结构2"中三角网格获取邻接面信息的方法。简单点说就是：遍历四面体集合，对每个四面体T，将T中每个点的邻接四面体集合中加入T。
2) 如何获取三角片集合？
从TetraMesh数据结构中的V和T信息得出F信息可以采用这样的思路实现：首先我们知道四面体集合T，遍历所有的四面体，对每个四面体我们都可以获得4个三角片组合，这些三角片组合的全集经过去重之后就形成了三角片集合F。具体地，去重的思路是排序去重，下面的实现定义了一个重载了运算符的结构，用来排序去重，最后大致的思路如函数InitInnnerTriangles所示，在此过程中，NeighborTetraPerTriangle成员也被处理完成，这样既获得了所有三角片的集合也获得所有三角片相邻的四面体信息：
struct IntIndexTriple { int X; int Y; int Z; int Index; IntIndexTriple(int x,int y,int z):X(x),Y(y),Z(z){} void Sort() { int max,min; if(X>Y) { max=X; min=Y; } else { max=Y; min=X; } if(Z>max) { X=min; Y=max; } else if(Z<min) { X=Z; Y=min; Z=max; } else { X=min; Y=Z; Z=max; } } bool operator<(const IntIndexTriple &right) const { if(X < right.X) return true; else if(X>right.X) return false; else if(Y < right.Y) return true; else if(Y >right.Y) return false; else if(Z<right.Z) return true; else if(Z>right.Z) return false; else return false; } bool operator==(const IntIndexTriple& right) const { return X==right.X&&Y==right.Y&&Z==right.Z; } }; void InitInnerTriangles() { if(Tedtrahedras.size()==0) return; std::vector<IntIndexTriple> templist; for(size_t i=0;i<Tedtrahedras.size();i++) { Tetrahedra& te=Tedtrahedras[i]; IntIndexTriple t1(te.PIndex[0],te.PIndex[2],te.PIndex[1]); IntIndexTriple t2(te.PIndex[0],te.PIndex[1],te.PIndex[3]); IntIndexTriple t3(te.PIndex[0],te.PIndex[3],te.PIndex[2]); IntIndexTriple t4(te.PIndex[1],te.PIndex[2],te.PIndex[3]); t1.Sort(); t2.Sort(); t3.Sort(); t4.Sort(); t1.Index=i; t2.Index=i; t3.Index=i; t4.Index=i; templist.push_back(t1); templist.push_back(t2); templist.push_back(t3); templist.push_back(t4); } std::sort(templist.begin(),templist.end()); int uniquesize=1; Triangle t0(templist[0].X,templist[0].Y,templist[0].Z); InnerTriangles.push_back(t0); std::vector<int> vec0;vec0.reserve(2); vec0.push_back(templist[0].Index); NeighborTetraPerTriangle.push_back(vec0); for(size_t i=1;i<templist.size();i++) { if(!(templist[i]==templist[i-1])) { uniquesize++; Triangle t(templist[i].X,templist[i].Y,templist[i].Z); InnerTriangles.push_back(t); std::vector<int> vec;vec.reserve(2); vec.push_back(templist[i].Index); NeighborTetraPerTriangle.push_back(vec); } else { NeighborTetraPerTriangle[uniquesize-1].push_back(templist[i].Index); } } }
3) 如何获得面体关系
上述求三角片集合的方法中已经顺便获得此关系。
4) 如何获得体体关系？
基于面体关系，方法描述：遍历三角片集合，对每个三角片两侧的四面体T1、T2，T1的邻接体集合加入T2，T2的邻接体集合加入T1。
5) 其他邻接关系与信息？
可以依次类推，比较容易，例如求E集合和EF关系，可类比上面求F集合和FT关系的方式，如果理解了就不难实现。
概念：限定三角剖分/限定四面体剖分
所谓限定三角剖分，是指给定一系列点集和一系列边，要对点集进行三角剖分并能保证剖分的结果边集合里包含所有给定的边。下图是限定三角剖分示意图。可以看出限定剖分与Delaunay三角剖分的区别。后者的输入只需点集即可，而前者需要输入点集+限定边集。详细的理论介绍可见tetgen文档。
 |
→ |  |
| 点集输入 | 正常Delaunay三角剖分输出 | |
 |
→ |  |
| 点集+限定边输入 | 限定剖分输出，限定边附近三角片不满足Delaunay准则 |
限定四面体剖分是限定三角剖分扩展到高维度的一个情况，其输入为空间点集+限定三角片集合。要求是在剖分的结果三角片中一定存在限定三角片。详细的理论介绍可见可见tetgen文档，下图是一个限定四面体剖分，在不同颜色的交界处存在一个给定曲面，而这个四面体剖分是严格遵守了在给定曲面处的限定特性的。
方法设计：如何使用tetgen库来实现限定四面体剖分
对本场景应用的需求细化描述应该是如下的几条：
tetgen库是一个做四面体剖分的开源库，基于C/C++，网上介绍实在太多，这里不多说，贴上两个链接：
为了方便下文的描述，这里先介绍一下其中用于测试的输入数据，如下表所示：
| 数据名 | 描述 | 截图 |
| sample3.ply |
人工数据其box范围 为range：(0,0,0, 49,49,49) |
 |
| Yhs_test_1.dat |
地层数据，range: (0,0,0, 7161.497070, 10037.639648, 6000.000000)
|
 |
| Yhs_test_2.dat |
地层数据 range: 0.000000, 2219.188300, 1088.959961, 9523.590300, 8527.667400, 3450.000000 |
 |
数据的下载地址在此。其中dat文件的格式需要注意是它不是一个正规的网格文件格式，而是类似这样的一坨东西：
range代表box的范围(xmin,ymin,zmin,xmax,ymax,zmax)，nface指其中的曲面个数，设为n，然后下面一行是每个曲面的名字和三角片个数m，然后在下面是m行，每行一个9元组表示一个三角片，这个三角片集合需要恢复拓扑成mesh之后才能做限定剖分，那么如何恢复这个拓扑呢，参见子问题1。
子问题1-输入三角网的构建/如何从三角片集合恢复三角网拓扑结构？
由于本应用中输入的限定曲面的形式不是类似ply这种格式的mesh，而是一系列散乱三角片构成的N行，故在使用tetgen前应恢复拓扑构建三角网。这个所谓散乱的三角片指的是三角片由一个九元组表示，九个double数分别对应其三个顶点的XYZ坐标值。这样的九元组组成的集合又被称作triangle soup，使用triangle soup表示的三角网格是没有拓扑信息的，即无法获知顶点的共用情况，通常需要将其转为具有拓扑信息的Mesh数据结构，本博客之前写过一篇关于顶点焊接(Vertex Welding)的文章简要讲述了三维数据场里的顶点焊接，但其局限于三维数据场中设计hash表，不太具有普遍性。这里就再介绍一个非hash的简单一点的焊接顶点恢复拓扑的方法。
方法的思路其实非常简单，triangle soup中，三角片都是唯一的，只是各自的顶点没有归一化处理，那么最简单的想法是将这些N个三角片的3N个顶点取出来到数组里，然后给这个数组进行去重操作，需要做的是排序+遍历一次进行归一，由于Mesh的三角片需要记录顶点的索引，故采用一些辅助结构记录下排序去重后的顶点索引即可。具体实现的代码如下c++类所示：
struct OringinalTriangle { Point3d P[3]; OringinalTriangle(double x0,double y0,double z0,double x1,double y1,double z1,double x2,double y2,double z2) { P[0].X=x0; P[0].Y=y0; P[0].Z=z0; P[1].X=x1; P[1].Y=y1; P[1].Z=z1; P[2].X=x2; P[2].Y=y2; P[2].Z=z2; } }; class OriginalTriangleWelder { public: static Mesh* WeldingTrianglesToMesh(std::vector<OringinalTriangle>& inputTriangles) { OriginalTriangleWelder otw(inputTriangles); return otw.GetWeldedMesh(); }//输入OringinalTriangle集合，输出Mesh private: struct LFloatTriple { double X; double Y; double Z; int Index; bool operator<(const LFloatTriple &right) const { if(X < right.X) return true; else if(X>right.X) return false; else if(Y < right.Y) return true; else if(Y >right.Y) return false; else if(Z<right.Z) return true; else if(Z>right.Z) return false; else return false; } };//存储点以及旧索引位置，因为要排序，所以做了个运算符重载 std::vector<OringinalTriangle>& otriangles; public: OriginalTriangleWelder(std::vector<OringinalTriangle>& inputTriangles):otriangles(inputTriangles) { } ~OriginalTriangleWelder() { } Mesh* GetWeldedMesh() { Mesh *retMesh=new Mesh(); std::vector<LFloatTriple> vec; std::vector<Triangle> triangles; triangles.reserve(otriangles.size()); vec.reserve(otriangles.size()*3); for(size_t i=0;i<otriangles.size();i++) { LFloatTriple ft0,ft1,ft2; ft0.X=otriangles[i].P[0].X; ft0.Y=otriangles[i].P[0].Y; ft0.Z=otriangles[i].P[0].Z; ft1.X=otriangles[i].P[1].X; ft1.Y=otriangles[i].P[1].Y; ft1.Z=otriangles[i].P[1].Z; ft2.X=otriangles[i].P[2].X; ft2.Y=otriangles[i].P[2].Y; ft2.Z=otriangles[i].P[2].Z; int count=vec.size(); ft0.Index=count; ft1.Index=count+1; ft2.Index=count+2; vec.push_back(ft0); vec.push_back(ft1); vec.push_back(ft2); Triangle t(count,count+1,count+2); triangles.push_back(t); }//建好点集与三角形集合 std::vector<int> map1,map2; //对原位置为i的顶点，通过map1[i]为排序后的新位置 //对排序后的顶点i map2[i]为焊接聚合后的新索引位置 map1.resize(otriangles.size()*3); map2.resize(otriangles.size()*3); std::sort(vec.begin(),vec.end()); for(size_t i=0;i<vec.size();i++) map1[vec[i].Index]=i; int index=0; std::vector<Point3d> points; Point3d p0(vec[0].X,vec[0].Y,vec[0].Z); points.push_back(p0); for(size_t i=1;i<vec.size();i++) { if(!(vec[i]<vec[i-1])&&!(vec[i-1]<vec[i])) { map2[i]=index; } else { index++; map2[i]=index; Point3d pi(vec[i].X,vec[i].Y,vec[i].Z); points.push_back(pi); } } for(size_t i=0;i<triangles.size();i++) { triangles[i].P0Index=map2[map1[triangles[i].P0Index]]; triangles[i].P1Index=map2[map1[triangles[i].P1Index]]; triangles[i].P2Index=map2[map1[triangles[i].P2Index]]; } retMesh->Vertices.swap(points); retMesh->Faces.swap(triangles); return retMesh; } };//类，完成从9元组三角形集合，转化为标准Mesh
子问题2-输入适配/如何做一个PLC结构作为tetgen的合法输入？
现在有了表示限定曲面的Mesh和Box，那么就可以开始用tetgen来做限定剖分了。一个重要的问题是输入，输入问题确实在一开始搞的笔者很囧，一直没明白怎么把复杂的限定曲面输入进去，一开始总是各种不合法，导致程序强退，后来终于试出了规律，且看下文分解。
 |
 |
| 限定剖分 | 限定剖分结果 |
tetgen的一个重要函数是tetrahedralize函数，其原型声明如下：
// tetrahedralize() Interface for using TetGen‘s library to generate // // Delaunay tetrahedralizations, constrained Delaunay// // tetrahedralizations, quality tetrahedral meshes.// // ‘in‘ is an object of ‘tetgenio‘ which contains a PLC you want to tetrahed-// // ralize or a previously generated tetrahedral mesh you want to refine. It // // must not be a NULL. ‘out‘ is another object of ‘tetgenio‘ for storing the // // generated tetrahedral mesh. It can be a NULL. If so, the output will be // // saved to file(s). If ‘bgmin‘ != NULL, it contains a background mesh which // // defines a mesh size function. // void tetrahedralize(tetgenbehavior *b, tetgenio *in, tetgenio *out, tetgenio *addin = NULL, tetgenio *bgmin = NULL);
其中我们需要用到的参数为tetgenbehavior *b, tetgenio *in, tetgenio *out, tetgenio *addin这四个，这四个参数的含义分别指代限定剖分选项参数结构、输入网格结构，输出剖分结构，附加点集输入结构。从上文的描述可知，限定剖分其实就是输入Box+Mesh，输出TetraMesh的过程，下面进一步解释如何将输入的用于执行限定剖分的Box+Mesh转为tetgenio结构。
首先对tetgen里需要常使用到的成员数据结构做个介绍：
| 数据结构 | 部分重要成员 | 类型 | 含义 |
| tetgenio | pointlist | double* | 存储点集，每三个double表示一个点，其中XYZ坐标依次输入数组。 |
| numberofpoints | int | 存储点集的大小，点集数组一般是这个数的3倍大。 | |
| facetlist | facet* | 面片集合 | |
| numberoffacets | int | 面片个数 | |
| facet | polygonlist | polygon* | 一个面片的多边形集合 |
| numberofpolygons | int | 多边形个数 | |
| polygon | vertexlist | int* | 一个多边形的顶点列表 |
| numberofvertices | int | 顶点列表长度 |
根据上表我们把我们自己定义的数据结构Box+Mesh的信息输入到这些成员上。我们分两个部分来说明即：输入点集和输入面片集合：
1.输入点集：
点集的输入包括三部分点的输入，这里必须假定三部分点集不包含有位置相同的点，这三部分分别是：①限定Mesh中的点；②Box的八个顶点；③后续插入的自定义点集；其中第1、2部分要按顺序依次输入参数in的pointlist结构中，pointlist是一个double数组，每3个数表示一个点的XYZ坐标。故相应的函数如下，其中points数组来自限定Mesh：
static void InitPoints(tetgenio& in,std::vector<Point3d>& points,Box3Double box) { in.numberofpoints=points.size()+8; in.pointlist=new REAL[in.numberofpoints*3]; for(size_t i=0;i<points.size();i++) { in.pointlist[3*i]=points[i].X; in.pointlist[3*i+1]=points[i].Y; in.pointlist[3*i+2]=points[i].Z; } InitBoxPoints(in,points,box); }//将限定曲面点和box顶点加入剖分
第3部分附加的点集需要设置在addin的pointlist里（tetgen的tetrahedralize函数在内部实现中，是先将BOX+限定Mesh进行限定四面体剖分之后，再采用addin里的点击对剖分结果进行插点与再剖分）。
static void InitAddin(tetgenio& addin,std::vector<Point3d>& insertedPoints) { addin.numberofpoints=insertedPoints.size(); addin.pointlist=new REAL[insertedPoints.size()*3]; for(size_t i=0;i<insertedPoints.size();i++) { addin.pointlist[3*i]=insertedPoints[i].X; addin.pointlist[3*i+1]=insertedPoints[i].Y; addin.pointlist[3*i+2]=insertedPoints[i].Z; } }//设置addin插入点
2.输入三角片：
输入三角片的函数很容易用错，这里重点讲述这一问题。限定曲面的facetlist由两大部分组成：①限定Mesh中的三角形；②Box的6个面。
2.1输入限定Mesh上的三角片
输入限定Mesh的三角片比较容易，即每一个facet对象对应一个限定三角片，故其输入方式如下：
in.firstnumber = 0; in.numberoffacets = faces.size()+6; in.facetlist = new tetgenio::facet[in.numberoffacets]; in.facetmarkerlist = new int[in.numberoffacets]; memset(in.facetmarkerlist,0,sizeof(int)*in.numberoffacets); tetgenio::facet *f;tetgenio::polygon *p; for(size_t i=0;i<faces.size();i++) { f = &in.facetlist[i]; f->numberofpolygons = 1; f->polygonlist = new tetgenio::polygon[f->numberofpolygons]; f->numberofholes = 0; f->holelist = NULL; p = &f->polygonlist[0]; p->numberofvertices = 3; p->vertexlist = new int[p->numberofvertices]; p->vertexlist[0] = faces[i].P0Index; p->vertexlist[1] = faces[i].P1Index; p->vertexlist[2] = faces[i].P2Index; }
需要解释的是facet具有一个polygonlist成员和一个holelist成员，后者在本文中暂时用不上，对于前者我们认为每个facet的polygon数为1，这个polygon的顶点数为3，代表一个三角片，然后在这个polygon的vertexlist成员里依次输入三个顶点索引即可。
2.2 输入Box六个面的限定facet
Box六个面的输入需要注意一下，首先若限定Mesh无点落在box上（包括棱上和面上），则可以将6个面作为6个facet输入，每个facet包含一个polygon，顶点数为4。若限定Mesh与Box有相交的情况，如下图，则上述输入方式会导致剖分失败（程序强行退出或者输出无结果）。
 |
 |
| 限定曲面与Box相交的情况 | |
经过笔者的各种调试分析，后来终于理解了为何这种情况会出问题，归根结底是因为输入的facet应该符合tetgen文档里对PLC的定义，通俗的讲，当限定Mesh与Box的表面长方形有相交部分的时候(相交的含义是有三角片的1~3个顶点正好落在这个长方形上)则不可以简单的把这个面当成一个长方形facet。那么如何输入呢，这就需要对限定Mesh和Box的六个面的相交情况进行计算。
对Box的六个面中的任意一个面F(V0,V1,V2,V3)，设其与限定mesh的相交情况如下，则其对应一个facet f，其polygonlist的长数为1+n，其中1为F边构成的外圈多边形，n为限定曲面在F上的相交线段(每个线段是一个polygon)。在下图中，由于限定Mesh在F上相交得<I1,I8>,<I8,I9>,<I8,I10>,<I10,I4>,<I9,I6> 5条线段，外加外圈多边形为<V0,I1,V1,I4,V2,I6,V3> (注意不再是<V0,V1,V2,V3>)，故这个f的numberofpolygons设为6，第一个polygon的点集输入对应的V0,I1,V1,I4,V2,I6,V3，点数为7；之后的polygon依次输入5个线段的点，点数为2。
值得注意的是：由F的边构成的第一个polygon顶点数不是4了，而是4+m，其中m为限定曲面的三角片在F的边上相交的点数。但这个点输入涉及一个顺序问题，一般是按时针顺序输入，即V0,I1,V1,I4,V2,I6,V3，或者V3,V0,I1,V1,I4,V2,I6 ，但不可不按时针顺序，故在求得了f上的所有交点后还需要进行时针排序。所谓时针排序其实可以简单的按照点的三个坐标之和来排。有的棱上坐标和顺序排，有的棱上坐标和逆序排，最后可以找出一个合法的点集顺序。
相应的代码实现中，遍历三角片集合找出与Box相交的所有线段和位于Box的棱上的点，其实现方式如下属代码所示。
static void InitFacets(tetgenio& in,std::vector<Point3d>& points,std::vector<Triangle>& faces,Box3Double box) { in.firstnumber = 0; in.numberoffacets = faces.size()+6; in.facetlist = new tetgenio::facet[in.numberoffacets]; in.facetmarkerlist = new int[in.numberoffacets]; memset(in.facetmarkerlist,0,sizeof(int)*in.numberoffacets); tetgenio::facet *f;tetgenio::polygon *p; for(size_t i=0;i<faces.size();i++) { f = &in.facetlist[i]; f->numberofpolygons = 1; f->polygonlist = new tetgenio::polygon[f->numberofpolygons]; f->numberofholes = 0; f->holelist = NULL; p = &f->polygonlist[0]; p->numberofvertices = 3; p->vertexlist = new int[p->numberofvertices]; p->vertexlist[0] = faces[i].P0Index; p->vertexlist[1] = faces[i].P1Index; p->vertexlist[2] = faces[i].P2Index; } std::vector<std::vector<Int16Double> > segmentsInbound; segmentsInbound.resize(6); std::vector<std::vector<int> > pointsInbound; pointsInbound.resize(12); InitBoxEdgePoints(in,points,faces,box,pointsInbound); InitBoxSegments(in,points,faces,box,segmentsInbound); InitBoxFacets(in,points,faces,box,segmentsInbound,pointsInbound); } static void InitBoxEdgePoints(tetgenio& in,std::vector<Point3d>& points,std::vector<Triangle>& faces,Box3Double box,std::vector<std::vector<int> >& boundPoints) { int ret[12]; for(size_t i=0;i<points.size();i++) { if(!InBox(points[i],box)) { printf("%f %f %f",points[i].X,points[i].Y,points[i].Z); //throw new exception(); } CheckPointAlignBoxEdge(points[i],box,ret); for(int j=0;j<12;j++) { if(ret[j]==1) { boundPoints[j].push_back(i); } } } }//检测是否有顶点落在box棱上 找到这些点 加入boundPoints的对应集合 static void InitBoxSegments(tetgenio& in,std::vector<Point3d>& points,std::vector<Triangle>& faces,Box3Double box,std::vector<std::vector<Int16Double> >& boundSegs) { int ret[6]; for(size_t i=0;i<faces.size();i++) { Triangle& t=faces[i]; Point3d& p0=points[t.P0Index]; Point3d& p1=points[t.P1Index]; Point3d& p2=points[t.P2Index]; if(!InBox(p0,box)||!InBox(p1,box)||!InBox(p2,box)) throw std::exception(); CheckSegAlignBoxFace(p0,p1,box,ret); for(int j=0;j<6;j++) { if(ret[j]==1) { Int16Double seg(t.P0Index,t.P1Index); boundSegs[j].push_back(seg); } } CheckSegAlignBoxFace(p1,p2,box,ret); for(int j=0;j<6;j++) { if(ret[j]==1) { Int16Double seg(t.P1Index,t.P2Index); boundSegs[j].push_back(seg); } } CheckSegAlignBoxFace(p2,p0,box,ret); for(int j=0;j<6;j++) { if(ret[j]==1) { Int16Double seg(t.P2Index,t.P0Index); boundSegs[j].push_back(seg); } } } }//检测是否有三角片的某些边落在box面上 找到这些线段 加入boundSegs的对应集合 static void InitBoxFacets(tetgenio& in,std::vector<Point3d>& points,std::vector<Triangle>& faces,Box3Double box,std::vector<std::vector<Int16Double> >& segInBox,std::vector<std::vector<int> >& pb) { tetgenio::facet *f; tetgenio::polygon *p; std::vector<Int16Double>* segs; int count=faces.size(); in.facetmarkerlist[count]=-1; in.facetmarkerlist[count+1]=-2; in.facetmarkerlist[count+2]=-3; in.facetmarkerlist[count+3]=-4; in.facetmarkerlist[count+4]=-5; in.facetmarkerlist[count+5]=-6; int p0=in.numberofpoints-8; int p1=in.numberofpoints-8+1; int p2=in.numberofpoints-8+2; int p3=in.numberofpoints-8+3; int p4=in.numberofpoints-8+4; int p5=in.numberofpoints-8+5; int p6=in.numberofpoints-8+6; int p7=in.numberofpoints-8+7; //if(p.X==x0&&p.Y==y0) 1 5 pointsInbound[0] //if(p.X==x1&&p.Y==y0) 2 6 pointsInbound[1] //if(p.X==x1&&p.Y==y1) 3 7 pointsInbound[2] //if(p.X==x0&&p.Y==y1) 4 8 pointsInbound[3] // //if(p.Y==y0&&p.Z==z0) 1 2 pointsInbound[4] //if(p.Y==y1&&p.Z==z0) 3 4 pointsInbound[5] //if(p.Y==y1&&p.Z==z1) 7 8 pointsInbound[6] //if(p.Y==y0&&p.Z==z1) 5 6 pointsInbound[7] // //if(p.X==x0&&p.Z==z0) 1 4 pointsInbound[8] //if(p.X==x1&&p.Z==z0) 2 3 pointsInbound[9] //if(p.X==x1&&p.Z==z1) 6 7 pointsInbound[10] //if(p.X==x0&&p.Z==z1) 5 8 pointsInbound[11] // Facet x=x0 f = &in.facetlist[count]; segs=&(segInBox[0]); f->numberofpolygons =1+segs->size(); f->polygonlist = new tetgenio::polygon[f->numberofpolygons]; f->numberofholes = 0; f->holelist = NULL; p = &f->polygonlist[0]; ProcessBoxIntersectedPoints(p0,pb[0],true,p4,pb[11],true,p7,pb[3],false,p3,pb[8],false,p,points); ProcessBoxSegments(f,p,segs); // Facet x=x2 f = &in.facetlist[count+1]; segs=&(segInBox[1]); f->numberofpolygons =1+segs->size(); f->polygonlist = new tetgenio::polygon[f->numberofpolygons]; f->numberofholes = 0; f->holelist = NULL; p = &f->polygonlist[0]; ProcessBoxIntersectedPoints(p1,pb[9],true,p2,pb[2],true,p6,pb[10],false,p5,pb[1],false,p,points); ProcessBoxSegments(f,p,segs); // Facet y=y0 f = &in.facetlist[count+2]; segs=&(segInBox[2]); f->numberofpolygons =1+segs->size(); f->polygonlist = new tetgenio::polygon[f->numberofpolygons]; f->numberofholes = 0; f->holelist = NULL; p = &f->polygonlist[0]; ProcessBoxIntersectedPoints(p0,pb[4],true,p1,pb[1],true,p5,pb[7],false,p4,pb[0],false,p,points); ProcessBoxSegments(f,p,segs); // Facet y=y2 f = &in.facetlist[count+3]; segs=&(segInBox[3]); f->numberofpolygons =1+segs->size(); f->polygonlist = new tetgenio::polygon[f->numberofpolygons]; f->numberofholes = 0; f->holelist = NULL; p = &f->polygonlist[0]; ProcessBoxIntersectedPoints(p2,pb[5],false,p3,pb[3],true,p7,pb[6],true,p6,pb[2],false,p,points); ProcessBoxSegments(f,p,segs); // Facet z=z0 f = &in.facetlist[count+4]; segs=&(segInBox[4]); f->numberofpolygons =1+segs->size(); f->polygonlist = new tetgenio::polygon[f->numberofpolygons]; f->numberofholes = 0; f->holelist = NULL; p = &f->polygonlist[0]; ProcessBoxIntersectedPoints(p0,pb[4],true,p1,pb[9],true,p2,pb[5],false,p3,pb[8],false,p,points); ProcessBoxSegments(f,p,segs); // Facet z=z2 f = &in.facetlist[count+5]; segs=&(segInBox[5]); f->numberofpolygons =1+segs->size(); f->polygonlist = new tetgenio::polygon[f->numberofpolygons]; f->numberofholes = 0; f->holelist = NULL; p = &f->polygonlist[0]; ProcessBoxIntersectedPoints(p4,pb[7],true,p5,pb[10],true,p6,pb[6],false,p7,pb[11],false,p,points); ProcessBoxSegments(f,p,segs); }
这一步看似复杂其实很好理解。综上所述，合法的使用tetgen进行限定四面体剖分需要有合法的输入，在我们这个应用场景下，除了输入限定曲面的每个三角片之外，还需要输入Box的6个面。需要注意的是当限定曲面与Box有相交的时候，注意不能再把Box的面当长方形输入，要寻找到棱上的交点和面上的线段再进行输入。
当成功输入限定曲面之后，若不插点进行四面体剖分，调用tetrahedralize函数的结果如下：
| sample3.ply |  |
| yhs_test_1.dat |  |
| yhs_test_2.dat |  |
子问题3-如何向剖分结果插入特定点集?/生成自定义插入点集并使用addin参数
生成插入点集
生成插入点集的要求是均匀，这样本人简单设计了一个插入规则点阵列的方法，设置一个可变的参数unitLen，表示插入点之间的横纵距离，然后按这个参数在体表面和内部均匀插点。
插点具体分三种场合，在box内部插的点，在box棱上插的点，在box面上插的点，这三种插点方法用程序很容易实现，以在边上插点为例。
|
+ |  |
= |  |
| 在限定剖分的结果中插入自定义点集，tetgen的主函数里addin参数的作用即在于此 | ||||
设置点的自由度
由于在后面去这些插入点可能会自适应移动，但并非所有的插入点都能在xyz三个方向自由移动（例如棱上插的点要是移动了可能破坏box形状）所以需要给插入的点区分不同的自由度。每个点都有自己的自由度，用一个三元整数元组表示这个点在xyz轴方向上的自由度，例如（1，1，0）表示这个点的x，y坐标可以变，但z坐标不能变。根据以上的插值方式，在对棱，面，box内部进行插点，特此规定如下：
综上所述插点逻辑比较简单，代码如下所示：
class RandomVertexGenerator { //1.向BOX内部按unitLen均匀插点 //2. 向box的面上按unitLen均匀插点 //3. 向box的棱上按unitLen均匀插点 //这样做的初衷会使得插得点分布均匀，并不会重合， 改进思路可以考虑消灭随机性，直接用统一策略均匀插点 // 此类一开始是想在box插点的时候随机一把，不能摆的太整齐 // 后来觉得没必要 就整齐的生成点阵列插进剖分后，内部平滑几下就会一定程度上自适应分布了 private: Box3Double box; double unitLen; public: RandomVertexGenerator(double unitLength,Box3Double box) { this->unitLen=unitLength; this->box=box; } ~RandomVertexGenerator(){} public: int GenerateInnerPoints(std::vector<Point3d>& plist,std::vector<Int16Triple>& typelist,Int16Triple v) { int count=0; int xlenCount=(int)(box.GetXLength()/unitLen); int ylenCount=(int)(box.GetYLength()/unitLen); int zlenCount=(int)(box.GetZLength()/unitLen); double xlen=box.GetXLength()/xlenCount; double ylen=box.GetYLength()/ylenCount; double zlen=box.GetZLength()/zlenCount; for(int i=0;i<xlenCount-1;i++) { double xmin=(i+0.5f)*xlen+box.Min3[0]; double xmax=((xmin+xlen)>box.Max3[0]?(box.Max3[0]):(xmin+xlen)); for(int j=0;j<ylenCount-1;j++) { double ymin=(j+0.5f)*ylen+box.Min3[1]; double ymax=((ymin+ylen)>box.Max3[1]?(box.Max3[1]):(ymin+ylen)); for (int k=0;k<zlenCount-1;k++) { double zmin=(k+0.5f)*zlen+box.Min3[2]; double zmax=((zmin+zlen)>box.Max3[2]?(box.Max3[2]):(zmin+zlen)); Point3d p=GetClusterRandomPoints(xmin,xmax,ymin,ymax,zmin,zmax); plist.push_back(p); typelist.push_back(v); count++; } } } return count; }//内部插点，插点的时候点不会落到面和棱上，严格在内部 int GeneratePlanePoints(std::vector<Point3d>& plist,int dem,double value,std::vector<Int16Triple>& typelist,Int16Triple v) { int count=0; int xlenCount=0,ylenCount=0; double xlen=0,ylen=0,xst=0,xed=0,yst=0,yed=0; if(dem==2) { xlenCount=(int)(box.GetXLength()/unitLen); ylenCount=(int)(box.GetYLength()/unitLen); xlen=box.GetXLength()/xlenCount; ylen=box.GetYLength()/ylenCount; xst=box.Min3[0]; xed=box.Max3[0]; yst=box.Min3[1]; yed=box.Max3[1]; } if(dem==1) { xlenCount=(int)(box.GetXLength()/unitLen); ylenCount=(int)(box.GetZLength()/unitLen); xlen=box.GetXLength()/xlenCount; ylen=box.GetZLength()/ylenCount; xst=box.Min3[0]; xed=box.Max3[0]; yst=box.Min3[2]; yed=box.Max3[2]; } if(dem==0) { xlenCount=(int)(box.GetYLength()/unitLen); ylenCount=(int)(box.GetZLength()/unitLen); xlen=box.GetYLength()/xlenCount; ylen=box.GetZLength()/ylenCount; xst=box.Min3[1]; xed=box.Max3[1]; yst=box.Min3[2]; yed=box.Max3[2]; } for(int i=0;i<xlenCount-1;i++) { double xmin=(i+0.5f)*xlen+xst; double xmax=((xmin+xlen)>xed?(xed):(xmin+xlen)); for(int j=0;j<ylenCount-1;j++) { double ymin=(j+0.5f)*ylen+yst; double ymax=((ymin+ylen)>yed?(yed):(ymin+ylen)); double x=GetClusterRandom(xmin,xmax); double y=GetClusterRandom(ymin,ymax); if(dem==2) { Point3d p(x,y,value);plist.push_back(p); typelist.push_back(v); count++; } if(dem==1) { Point3d p(x,value,y);plist.push_back(p); typelist.push_back(v); count++; } if(dem==0) { Point3d p(value,x,y);plist.push_back(p); typelist.push_back(v); count++; } } } return count; }//平面内部插点，插点的时候点不会落到棱上，严格在面内部 int GenerateBoxPoints(std::vector<Point3d>& plist,std::vector<Int16Triple>& typelist,Int16Triple v) { //double xlen=box.GetXLength(); double x0=box.Min3[0],x1=box.Max3[0]; //double ylen=box.GetYLength(); double y0=box.Min3[1],y1=box.Max3[1]; //double zlen=box.GetZLength(); double z0=box.Min3[2],z1=box.Max3[2]; Point3d p0(x0,y0,z0); Point3d p1(x1,y0,z0); Point3d p2(x0,y1,z0); Point3d p3(x1,y1,z0); Point3d p4(x0,y0,z1); Point3d p5(x1,y0,z1); Point3d p6(x0,y1,z1); Point3d p7(x1,y1,z1); plist.push_back(p0); plist.push_back(p1); plist.push_back(p2); plist.push_back(p3); plist.push_back(p4); plist.push_back(p5); plist.push_back(p6); plist.push_back(p7); for(int i=0;i<8;i++) typelist.push_back(v); return 8; }//获取BOX的八个顶点，不属于插点，别处会用到 int GenerateLinePoints(std::vector<Point3d>& plist,int demParalTo,double d1,double d2,std::vector<Int16Triple>& typelist,Int16Triple v) { int count=0; int xlenCount=0; double xlen=0,xst=0,xed=0; if(demParalTo==0) { xlenCount=(int)(box.GetXLength()/unitLen); xlen=box.GetXLength()/xlenCount; xst=box.Min3[0]; xed=box.Max3[0]; } if(demParalTo==1) { xlenCount=(int)(box.GetYLength()/unitLen); xlen=box.GetYLength()/xlenCount; xst=box.Min3[1]; xed=box.Max3[1]; } if(demParalTo==2) { xlenCount=(int)(box.GetZLength()/unitLen); xlen=box.GetZLength()/xlenCount; xst=box.Min3[2]; xed=box.Max3[2]; } for(int i=0;i<xlenCount-1;i++) { double xmin=(i+0.5f)*xlen+xst; double xmax=((xmin+xlen)>xed?(xed):(xmin+xlen)); if(demParalTo==0) { Point3d p(GetClusterRandom(xmin,xmax),d1,d2); plist.push_back(p); typelist.push_back(v); count++; } if(demParalTo==1) { Point3d p(d1,GetClusterRandom(xmin,xmax),d2); plist.push_back(p); typelist.push_back(v); count++; } if(demParalTo==2) { Point3d p(d1,d2,GetClusterRandom(xmin,xmax)); plist.push_back(p); typelist.push_back(v); count++; } } return count; }//棱内部插点，插点的时候点不会落到棱端点上，严格在棱内部 private: static double GetRandom(double st,double ed) { /*int r=rand()%1000; double x=r*(ed-st)/1000+st; return x;*/ return (st+ed)/2.0f; }//获得st到ed间的随机值 static double GetClusterRandom(double st,double ed,double clusterRate=0.35f) { if(clusterRate<0||clusterRate>=0.49f) return (st+ed)/2; else { double len=ed-st; st+=clusterRate*len; ed-=clusterRate*len; return GetRandom(st,ed); } }//获得st+clusterRate*len到ed-clusterRate*len间的随机值，做这个函数是想让随机值别离st和ed太近 static Point3d GetClusterRandomPoints(double xmin,double xmax,double ymin,double ymax,double zmin,double zmax) { return Point3d(GetClusterRandom(xmin,xmax), GetClusterRandom(ymin,ymax), GetClusterRandom(zmin,zmax) ); }//在参数指定的小box内获取一个随机点 };//插点类，根据需求进行插点，原始需求为向BOX内以尺度unitLen均匀插点，为保证插点质量，将逻辑分裂为下列3项，
当采用addin参数，按照上文点集的InitAddin函数设置插入点之后，产生的效果如下：
| yhs_test_1.dat |  |
| yhs_test_2.dat |  |
子问题4-如何使得插入点自适应分布?/内部平滑的实现
插入的点阵是一个规则的点阵，在与限定曲面的点何在一起显得极为不自然，因此需要自适应分布，其示意图大致如下图所示：
|
 |
| 内部平滑输入 | 内部平滑输出 |
自适应分布可以采用类似网格平滑的思路来解决，只是输入不再是表征表面的三角网格，而是四面体网格。四面体网格中其实仍然有邻接点的概念，所以我们不妨对这些点也采用拉普拉斯平滑算子进行处理。其处理方式即为：将每个点移动到其邻接点的平均位置。
其实现方式如下：
static void ExecuteInnerSmoothing(TetraMesh& mesh,Mesh& surface,int iteration,std::vector<Int16Triple>& optypes) { std::vector<Int16Triple> alltypes; size_t fixcount=surface.Vertices.size()+8; if(fixcount+optypes.size()==mesh.Vertices.size()) { alltypes.resize(mesh.Vertices.size()); for(size_t i=0;i<mesh.Vertices.size();i++) { if(i<fixcount) alltypes[i]=Int16Triple(0,0,0); else alltypes[i]=optypes[i-fixcount]; } } else { printf("woca!\n"); throw std::exception(); } std::vector<double> vd1array; std::vector<double> vd2array; std::vector<double> vd3array; std::vector<short> numbers; size_t vcount=mesh.Vertices.size(); size_t tcount=mesh.Tedtrahedras.size(); vd1array.resize(mesh.Vertices.size(),0.0f); vd2array.resize(mesh.Vertices.size(),0.0f); vd3array.resize(mesh.Vertices.size(),0.0f); numbers.resize(mesh.Vertices.size(),0); for(int c=0;c<iteration;c++) { for(size_t i=0;i<tcount;i++) { Tetrahedra& t=mesh.Tedtrahedras[i]; Point3d p0=mesh.Vertices[t.PIndex[0]]; Point3d p1=mesh.Vertices[t.PIndex[1]]; Point3d p2=mesh.Vertices[t.PIndex[2]]; Point3d p3=mesh.Vertices[t.PIndex[3]]; int P0Index=t.PIndex[0]; int P1Index=t.PIndex[1]; int P2Index=t.PIndex[2]; int P3Index=t.PIndex[3]; vd1array[P0Index]+=p1.X;vd2array[P0Index]+=p1.Y;vd3array[P0Index]+=p1.Z; vd1array[P0Index]+=p2.X;vd2array[P0Index]+=p2.Y;vd3array[P0Index]+=p2.Z; vd1array[P0Index]+=p3.X;vd2array[P0Index]+=p3.Y;vd3array[P0Index]+=p3.Z; vd1array[P1Index]+=p0.X;vd2array[P1Index]+=p0.Y;vd3array[P1Index]+=p0.Z; vd1array[P1Index]+=p2.X;vd2array[P1Index]+=p2.Y;vd3array[P1Index]+=p2.Z; vd1array[P1Index]+=p3.X;vd2array[P1Index]+=p3.Y;vd3array[P1Index]+=p3.Z; vd1array[P2Index]+=p1.X;vd2array[P2Index]+=p1.Y;vd3array[P2Index]+=p1.Z; vd1array[P2Index]+=p0.X;vd2array[P2Index]+=p0.Y;vd3array[P2Index]+=p0.Z; vd1array[P2Index]+=p3.X;vd2array[P2Index]+=p3.Y;vd3array[P2Index]+=p3.Z; vd1array[P3Index]+=p1.X;vd2array[P3Index]+=p1.Y;vd3array[P3Index]+=p1.Z; vd1array[P3Index]+=p0.X;vd2array[P3Index]+=p0.Y;vd3array[P3Index]+=p0.Z; vd1array[P3Index]+=p2.X;vd2array[P3Index]+=p2.Y;vd3array[P3Index]+=p2.Z; numbers[P0Index]+=3; numbers[P1Index]+=3; numbers[P2Index]+=3; numbers[P3Index]+=3; } for(size_t i=0;i<vcount;i++) { if(numbers[i]!=0) { if(alltypes[i].X==1) mesh.Vertices[i].X=vd1array[i]/numbers[i]; if(alltypes[i].Y==1) mesh.Vertices[i].Y=vd2array[i]/numbers[i]; if(alltypes[i].Z==1) mesh.Vertices[i].Z=vd3array[i]/numbers[i]; } vd1array[i]=0.0f; vd2array[i]=0.0f; vd3array[i]=0.0f; numbers[i]=0; } } //内部平滑实现，原理同拉普拉斯平滑，把每个点往邻接点的平均位 //置移动，不同的是根据alltypes，不同点允许改的坐标不都一样 // Box顶点 和限定曲面点完全不能动 // 棱插点可以在棱上滑动 // 面插点可在面上动 // 内部插点可随便动 }
结合拉普拉斯平滑算法的博客"几种网格平滑算法的实现"，本步骤不难理解。经过此步骤后，四面体剖分对测试数据的处理结果如下：
子问题5-拓扑分离/如何将限定曲面处的共享顶点分离开？
拓扑分离不属于四面体剖分的范畴，但在本文的应用场景下有此需求，所以花了一点心思编写了这部分代码，将限定曲面两侧的四面体的顶点从共享状态弄成分开状态。本来针对这一问题没有思路，之后在与mgm同学讨论的过程中收到启发，实现了这一功能。
实现拓扑分离的思路其实就是对限定曲面上的每个顶点进行逐点处理的过程，可以按如下的图标所示的思路进行处理（图是二维示意图，但三维情形完全类似）：
| 步骤 | 示意图 | 步骤说明 |
| 1 |  |
左边这是剖分结果。 |
| 2 |  |
首先对剖分结果中限定曲面上的点进行标记 对每个限定曲面上的点，求出相邻四面体。 |
| 3 |  |
根据体体关系（被限定曲面分开的强制不相邻）， 对每个顶点的四面体集合分成若干组。 |
| 4 |  |
根据组数求出每个顶点需要的额外副本个数。 |
| 5 |  |
补充副本个数个顶点追加到点集末尾。 |
| 6 |  |
将不同组的四面体的改点指向新的副本点。 |
其实现的相应代码如下：
class VertexSpliter3d { //拓扑分离策略，根据MGM的思路实现 // 在每个限定顶点点处，找到限定三角片作为边界，然后找出这些边界区分出的n个区域，然后在点集中补上n-1个点 // ，让n-1个区域的四面体们的这个点索引指向那些补上的对应的新点 class MassiveIndexedHash { //一种自定义的特殊hash表，使用场合为key的范围为[0-比较大但内存能接受的n], // 但哈希set的次数不多，此表又需要经常清空以支持大量复用，所以清空所有key要快 // 设计思路是做个int[n]存一堆默认值，即table ，再附带一个记录key的records容器 // 每次set 除了设table[key]=v之外，同时把key push到records里 // 这样每次清表的时候，按records里的key把table里的值恢复成默认再清空records就行 // 不设records记录key的话，每次clear遍历n大的数组搞毛 private: std::vector<int> records; std::vector<int> table; int defaultValue; public: MassiveIndexedHash(int size,int v) { table.resize(size,v); this->defaultValue=v; } ~MassiveIndexedHash(){} void SetValue(int key,int v) { if(key<0||(size_t)key>table.size()) throw std::exception(); if(defaultValue!=v) { records.push_back(key); table[key]=v; } } int GetValue(int key) { if(key<0||(size_t)key>table.size()) throw std::exception(); return table[key]; } void ClearTable() { for(size_t i=0;i<records.size();i++) table[records[i]]=defaultValue; records.clear(); } int GetCount() { return (int)records.size(); } }; struct ConnectedVF { int pcenterIndex; std::vector<int> teds; int targetpIndex; ConnectedVF() { pcenterIndex=-1; targetpIndex=-1; } };//结构，存储点索引和其所在的部分体的集合，用于表示被限定面分开的一个联通区域 private: static inline int Max(int a,int b,int c) { int ab=(a>b?a:b); return ab>c?ab:c; } static inline int Min(int a,int b,int c) { int ab=(a<b?a:b); return ab<c?ab:c; } static void AppendPoints(std::vector<Point3d>&points,int count,Point3d v) { if(count==0) return; for(int i=0;i<count;i++) { points.push_back(v); } }//向数组末尾push count个一样的v static void ChangeTetraPoints(Tetrahedra* t,int orindex,int newindex) { for(int i=0;i<4;i++) { if(t->PIndex[i]==orindex) { t->PIndex[i]=newindex; return; } } printf("error!\n"); throw std::exception(); }//将Tetrahedra中一个顶点的索引变成新值 static bool TriangleEqual(Triangle& t1,Triangle&t2) { int sum1=t1.P0Index+t1.P1Index+t1.P2Index; int sum2=t2.P0Index+t2.P1Index+t2.P2Index; if(sum1!=sum2) return false; int max1=Max(t1.P0Index,t1.P1Index,t1.P2Index); int max2=Max(t2.P0Index,t2.P1Index,t2.P2Index); if(max1!=max2) return false; int min1=Min(t1.P0Index,t1.P1Index,t1.P2Index); int min2=Min(t2.P0Index,t2.P1Index,t2.P2Index); if(min1!=min2) return false; return true; }//检查两个Triangle对象是否相等 private: TetraMesh& mesh; std::vector<Triangle>& Constrains; std::vector<int> ConstrainsMapToTriangleIndex;//第i个Constrains集合里的三角形在三角形集合里的索引 std::vector<int> TriangleMapToConstrainsIndex;//第i个三角形集合里的三角形在Constains里的索引 std::vector<bool> VBoundFlag; std::vector<int> VTypeFlag;//记录每个顶点在多少个限定曲面上 std::vector<int> VRegionRes; std::vector<std::vector<int> > VFAdjs;//点面关系 std::vector<std::vector<int> > TTAdjs;//体体关系 std::vector<Triangle> InnerTriangles; std::vector<std::vector<int> > FTAdjs;//面体关系 std::vector<std::vector<int> > VTAdjs;//点体关系 private: std::vector<bool> tempmark;//循环中复用 std::vector<ConnectedVF> regions;//循环中复用 std::queue<int> queue;//循环中复用 MassiveIndexedHash hashtable;//循环中复用 public: VertexSpliter3d(TetraMesh& m,Mesh& constrains):mesh(m),Constrains(constrains.Faces),hashtable(m.Tedtrahedras.size(),-1) { } ~VertexSpliter3d() { } private: void InitVT() { VTAdjs.resize(mesh.Vertices.size()); for(size_t i=0;i<mesh.Tedtrahedras.size();i++) { Tetrahedra& t=mesh.Tedtrahedras[i]; if(VTypeFlag[t.PIndex[0]]>0) VTAdjs[t.PIndex[0]].push_back(i); if(VTypeFlag[t.PIndex[1]]>0) VTAdjs[t.PIndex[1]].push_back(i); if(VTypeFlag[t.PIndex[2]]>0) VTAdjs[t.PIndex[2]].push_back(i); if(VTypeFlag[t.PIndex[3]]>0) VTAdjs[t.PIndex[3]].push_back(i); } }//初始化点体关系 void InitTypeFlags() { VTypeFlag.resize(mesh.Vertices.size(),0); for(size_t i=0;i<Constrains.size();i++) { VTypeFlag[Constrains[i].P0Index]++; VTypeFlag[Constrains[i].P1Index]++; VTypeFlag[Constrains[i].P2Index]++; } VRegionRes.resize(mesh.Vertices.size(),0); }//初始化点与限定曲面的关系，见VTypeFlag定义 void InitBoundFlags() { VBoundFlag.resize(mesh.Vertices.size(),false); for(size_t i=0;i<FTAdjs.size();i++) { if(FTAdjs[i].size()==1) { Triangle &t=InnerTriangles[i]; VBoundFlag[t.P0Index]=true; VBoundFlag[t.P1Index]=true; VBoundFlag[t.P2Index]=true; } } }//初始化外表面三角片顶点标记 void InitInnerTriangle_FT() { if(mesh.InnerTriangles.size()!=0&&mesh.NeighborTetraPerTriangle.size()!=0) { this->FTAdjs.swap(mesh.NeighborTetraPerTriangle); this->InnerTriangles.swap(mesh.InnerTriangles); } else { mesh.InnerTriangles.clear(); mesh.NeighborTetraPerTriangle.clear(); mesh.InitInnerTriangles(); this->FTAdjs.swap(mesh.NeighborTetraPerTriangle); this->InnerTriangles.swap(mesh.InnerTriangles); } }//根据V与T 求出F集合 参见InitInnerTriangles函数 void InitVF() { VFAdjs.resize(mesh.Vertices.size()); for(size_t i=0;i<InnerTriangles.size();i++) { VFAdjs[InnerTriangles[i].P0Index].push_back(i); VFAdjs[InnerTriangles[i].P1Index].push_back(i); VFAdjs[InnerTriangles[i].P2Index].push_back(i); } }//初始化点三角片关系 void InitConstrainsMap() { TriangleMapToConstrainsIndex.resize(InnerTriangles.size(),-1); for(size_t i=0;i<Constrains.size();i++) { Triangle& con=Constrains[i]; std::vector<int>& adjs=VFAdjs[con.P0Index]; for(size_t j=0;j<adjs.size();j++) { int findex=adjs[j]; Triangle &f=InnerTriangles[findex]; if(TriangleEqual(con,f)) { TriangleMapToConstrainsIndex[findex]=i; } } } ConstrainsMapToTriangleIndex.resize(Constrains.size(),-1); for(size_t i=0;i<TriangleMapToConstrainsIndex.size();i++) { int indexInCon=TriangleMapToConstrainsIndex[i]; if(indexInCon!=-1) { ConstrainsMapToTriangleIndex[indexInCon]=i; } } } void InitTT() { TTAdjs.resize(mesh.Tedtrahedras.size()); for(size_t i=0;i<InnerTriangles.size();i++) { std::vector<int>& fts=FTAdjs[i]; if(fts.size()==2) { int ti1=fts[0]; int ti2=fts[1]; TTAdjs[ti1].push_back(ti2); TTAdjs[ti2].push_back(ti1); } if(fts.size()==0||fts.size()>2) throw std::exception(); } }//初始化TT关系，根据FT关系推出 void InitTT_Type() { TTAdjs.resize(mesh.Tedtrahedras.size()); for(size_t i=0;i<InnerTriangles.size();i++) { if(TriangleMapToConstrainsIndex[i]!=-1) continue; else { std::vector<int>& fts=FTAdjs[i]; if(fts.size()==2) { int ti1=fts[0]; int ti2=fts[1]; TTAdjs[ti1].push_back(ti2); TTAdjs[ti2].push_back(ti1); } if(fts.size()==0||fts.size()>2) throw std::exception(); } } }//初始化TT关系，根据FT关系推出，不同之处在于这个函数考虑限定曲面，认为限定曲面F两侧的T不相邻 int GetInsertedNum() { int sum=0; for(size_t i=0;i<VTypeFlag.size();i++) { if(!VBoundFlag[i]) sum+=(VTypeFlag[i]-1); else sum+=VTypeFlag[i]; } return sum; }//计算出拓扑分离总共需要增加多少副本点 void InitRegions(int i,bool check=false) { regions.clear();regions.reserve(VTypeFlag[i]); std::vector<int>& allTetra=VTAdjs[i]; tempmark.clear();tempmark.resize(allTetra.size(),false); hashtable.ClearTable(); for(size_t j=0;j<allTetra.size();j++) hashtable.SetValue(allTetra[j],j); for(size_t j=0;j<allTetra.size();j++) { if(!tempmark[j]) { ConnectedVF region; region.pcenterIndex=i; region.teds.push_back(allTetra[j]); queue.push(j); tempmark[j]=true; while(!queue.empty()) { int k=queue.front(); queue.pop(); int tindex=allTetra[k]; std::vector<int>& TTNeighb=TTAdjs[tindex]; for(size_t n=0;n<TTNeighb.size();n++) { int nj=hashtable.GetValue(TTNeighb[n]); if(nj!=-1) { if(!tempmark[nj]) { tempmark[nj]=true; region.teds.push_back(allTetra[nj]); queue.push(nj); } } } } regions.push_back(region); } } if(check) { size_t sum=0; for(size_t j=0;j<regions.size();j++) { sum+=regions[j].teds.size(); } if(sum!=VTAdjs[i].size()) { printf("region fail!\n"); throw std::exception(); } }//check 是debug用的参数 }//对顶点i，求取限定曲面分开的各个邻接四面体子集，使用了广度遍历算法 void MarshalNewPoints(int index) { int orsize=mesh.Vertices.size()-regions.size()+1; regions[0].targetpIndex=index; for(size_t i=1;i<regions.size();i++) { regions[i].targetpIndex=orsize+i-1; std::vector<int>& facesInthisregion=regions[i].teds; for(size_t j=0;j<facesInthisregion.size();j++) { ChangeTetraPoints(&(mesh.Tedtrahedras[facesInthisregion[j]]),regions[i].pcenterIndex,regions[i].targetpIndex); } } }//将顶点i处限定曲面分开的不同类的四面体的顶点指向i的副本 void ProcessSingleVertex(int i) { InitRegions(i,true);//求取顶点i的邻接四面体分类后集合 VRegionRes[i]=regions.size(); AppendPoints(mesh.Vertices,regions.size()-1,mesh.Vertices[i]);//将副本点加入mesh点集 MarshalNewPoints(i);//将顶点i的各类四面体的顶点指向副本点 }//对每个点执行拓扑分离 public: void ProcessSpliting() { InitTypeFlags(); InitVT(); InitInnerTriangle_FT(); InitBoundFlags(); InitVF(); InitConstrainsMap(); InitTT_Type(); //以上为初始化工作 mesh.Vertices.reserve(mesh.Vertices.size()+GetInsertedNum()); size_t orcount=mesh.Vertices.size(); for(size_t i=0;i<orcount;i++) { if(VTypeFlag[i]>0) { ProcessSingleVertex(i); } } //以上为对每个点拓扑分离 std::vector<std::vector<int> >().swap(TTAdjs); std::vector<bool>().swap(VBoundFlag); std::vector<int>().swap(VTypeFlag); std::vector<int>().swap(VRegionRes); std::vector<std::vector<int> >().swap(VFAdjs); std::vector<std::vector<int> >().swap(TTAdjs); std::vector<Triangle>().swap(InnerTriangles); std::vector<std::vector<int> >().swap(FTAdjs ); std::vector<std::vector<int> >().swap(VTAdjs); std::vector<Triangle>().swap(mesh.InnerTriangles); std::vector<std::vector<int> >().swap(mesh.NeighborTetraPerTriangle); //以上为清除vector空间 mesh.InitInnerTriangles();//重新计算三角片集合 } };
经过拓扑分离的处理，若限定曲面能将Box的区域完全分隔开，则通过找联通区域的算法可以给不同的区域染上不同的颜色，测试结果的展示如下图：
 |
 |
 |
 |
其他需要注意的问题
项目工程下载：
本文的项目下载地址为：https://github.com/chnhideyoshi/ConstrainedTetrahedralization
原文链接为：http://www.cnblogs.com/chnhideyoshi/p/ConstrainedTetrahedralization.html 转载请注明出处。
标签：
原文地址：http://www.cnblogs.com/chnhideyoshi/p/ConstrainedTetrahedralization.html
