Tetgen学习，算法篇
Tetgen学习，算法篇
一、Delaunay 四面体化
tetgen使用了两种增量算法来生成3D PLC 的Delaunay constructing tetrahedralizations
the Bowyer-Watson (B-W) algorithm and the incremental flip algorithm of Edelsbrunner and Shah
函数incrementalflip()实现了[Edelsbrunner and Shah, 1996]的翻转算法。它翻转一个由本地非delaunay面组成的队列(以任意顺序)。当Delaunay四面化以每次增加一个顶点的增量方式构建时，成功得到了保证。
函数locate()在当前DT中查找包含新点的四面体。它采用了一种简单的随机行走算法:从DT中的任意一个四面体出发，每次访问一个四面体找到目的地，如果有多个选择，则随机选择一个四面体。由于Edelsbrunner的无环定理，该算法终止。
选择一个好的四面体对行走的速度至关重要。TetGen最初使用无需预处理的快速随机点定位。”第12届ACM计算几何研讨会论文集，. 它首先随机抽样DT中的几个四面体，然后选择最接近的一个开始行走。
上述算法会随着点数量的增加而大大减慢下载速度，如果对点进行预先排序，使空间中附近的点具有附近的指标，然后按此顺序将点相加，则可以在常数时间内确定点的位置。他们沿着三维希尔伯特曲线对这些点进行排序。
函数hilbert_sort3()沿着3D希尔伯特曲线对一组3D点进行排序。它根据映射到点集边界框的子框的希尔伯特指标递归地分割一个点集。希尔伯特指数是用巴兹算法在1971年计算出来的。在Hamilton, C.的论文“紧凑型希尔伯特指数”，技术报告CS-2006-07，计算机科学，Dalhousi大学，2006(第2节)中可以找到这个算法的一个很好的阐述。我的实现还引用了Stever Witham的“希尔伯特行走”的性能(希望它仍然可以在http://www上找到)。tiac。net/sw/2008/10/Hilbert/)TetGen使用Boissonnat, j . d论文中的方法对点进行排序。Devillers, 0。和Hornus,年代。“Delaunay三角剖分和Delaunay图的增量构造”，《第25届ACM计算几何研讨会论文集》，2009。首先使用Amenta et al 2003的有偏随机插入排序(BRIO)将点随机分组，然后沿三维希尔伯特曲线对每个子组中的点进行排序。以这种顺序插入点可以确保点在域中的随机“分散”，而对每个子集进行排序则提供了局部性。
二、四面体化的发展历程
step1:Chew and Ruppert 首次提出生成高质量网格，他们的工作保证了三角网格的最小角，并且在相同尺寸中生成的网格质量最好
step2：Shewchuk将Rupport的理论推广到三维空间，他的工作保证了生成的四面体网格的半径边缘比得到了限制，但是他却不能移除slivers（一种扁平的四面体，半径边缘比良好但是体积和最小二面角都很差），同时如果输入的PLC含有尖锐特性时算法不能保证终止
step3：TetGen使用基本的Delaunay细化方案插入斯坦纳点。使用的算法具有一定的终止性，易于保留清晰的特征。该网格在网格域的大部分区域具有良好的质量(与Shewchuk的Delaunay精化算法相同)。此外，它支持根据(各向同性)网格大小函数生成自适应网格。
三、四面体网格细化代码
void makesegmentendpointsmap();
/*
创建一个从段到其端点的映射。映射保存在数组' segmentendpointslist'中。
这个数组的长度是段数的两倍。每个段分配一个唯一的索引(从0开始)。
*/
REAL set_ridge_vertex_protecting_ball(point);
/*
Calculate the protecting ball for a given ridge vertex.
*/
REAL get_min_angle_at_ridge_vertex(face* seg);
/*
Calculate the minimum face angle at a given ridge vertex.
*/
REAL get_min_diahedral_angle(face* seg);
/*
Calculate the minimum (interior) dihedral angle a given segment.
*/
void create_segment_info_list();
/*
Calculate the minimum dihedral angle at a given segment.
segment_info_list = new double[segmentendpointslist_length * 4];
- [0] min_dihedral_angle (degree) at this segment,
- [1] min_protect_cylinder_radius at this segment (for bookkeeping only),
- [2] min_seg_seg_angle (degree) at its endpoint [0],
- [3] min_seg_seg_angle (degree) at its endpoint [1].
这个函数必须在makesegmentendpointsmap()之后调用。计算唯一段(segmentendpointslist_length)的数量。
*/
void makefacetverticesmap();
/*
创建从facet到其顶点的映射。所有的facet都将被索引(从0开始)。
映射保存在两个全局数组中:'idx2facetlist'和' facetverticeslist'
*/
void create_segment_facet_map();
/*
Create the map from segments to adjacent facet
•创建从段到相邻面的映射。
*/
int ridge_vertices_adjacent(point, point);
/*
检查两个脊顶点是否由输入段连接。
*/
int facet_ridge_vertex_adjacent(face *, point);
/*
检查一个面和一个脊顶点是否与一个输入段相邻。
*/
int segsegadjacent(face *, face *);
/*
检查两个线段，或者一个线段和一个面，或者两个面是否相邻。
*/
int segfacetadjacent(face *checkseg, face *checksh);
/*
检查一个线段和一个面是否相邻。
*/
int facetfacetadjacent(face *, face *);
/*
检查两个面是否相邻。
*/
bool is_sharp_segment(face* seg);
/*
检查给定的线段是否含有尖锐特性。
*/
bool does_seg_contain_acute_vertex(face* seg);
/*
检查给定线段的端点是否为尖角。
*/
bool create_a_shorter_edge(point steinerpt, point nearpt);
/*
我们能否在两个给定顶点之间创建一条边(比minedgelength短)?
*/
void enqueuesubface(memorypool*, face*);
/*
子面或子段队列进行入侵检查。
*/
void enqueuetetrahedron(triface*);
/*
对一个四面体队列进行质量检查。
*/
bool check_encroachment(point pa, point pb, point checkpt);
/*
检查一个给定的点是否侵占线段。
checkpt不应该是假点（dummy point）
*/
bool check_enc_segment(face *chkseg, point *pencpt);
/*
检测给定的段是否被侵占
*/
bool get_steiner_on_segment(face* seg, point encpt, point newpt);
/*
获取斯坦纳点来分割给定的线段。
*/
bool split_segment(face *splitseg, point encpt, REAL *param, int qflag, int, int*);
/*
分割一个给定的段。
// If param != NULL, it contains the circumcenter and its insertion radius //
// of a bad quality tetrahedron. Tghis circumcenter is rejected since it //
// encroaches upon this segment. //
如果param != NULL，它包含一个质量不好的四面体的周心及其插入半径。这个圆周被拒绝，因为它侵犯了这个部分。
*/
void repairencsegs(REAL *param, int qflag, int chkencflag);
/*
修复侵蚀的(子)段。
*/
bool get_subface_ccent(face *chkfac, REAL *ccent);
/*
计算给定子面的直径圆周圆心。
*/
bool check_enc_subface(face *chkfac, point *pencpt, REAL *ccent, REAL *radius);
/*
检查给定子面是否被侵占。
*/
bool check_subface(face *chkfac, REAL *ccent, REAL radius, REAL *param);
/*
给定的子面形状是否糟糕(依据半径-边缘比)
*/
void enqueue_subface(face *bface, point encpt, REAL *ccent, REAL *param);
/*
将形状不佳的子面推入优先队列。
*/
badface* top_subface();
/*
返回队列头的子面。
*/
void dequeue_subface();
/*
从优先级队列中弹出一个形状糟糕的子面。
*/
void parallel_shift(point pa, point pb, point pc, point pt, REAL* ppt);
/*
沿着三角形的法线平行移动一个三角形。
*/
enum locateresult locate_on_surface(point searchpt, face* searchsh);
/*
在一个面中定位一个顶点。
求一个包含给定点的四面体。
从' searchtet'开始搜索，假设从' searchtet'的顶点到查询点' searchpt'有一条线段L，并通过遍历所有与L相交的面来简单地走向' searchpt'。补全时，' searchtet'是一个包含' searchpt'的四面体。
返回值为以下情况之一:
ONVERTEX，搜索点位于searchtet的原点
ONEDGE，搜索点位于' searchtet'的边缘上
ONFACE，搜索点位于' searchtet'的面上
inet，搜索点位于' searchtet '的内部
OUTSIDE，搜索点位于网格外部。
' searchtet'是对于搜索点可见的。
*/
bool split_subface(face *splitfac, point encpt, REAL *ccent, REAL*, int, int, int*);
/*
分割一个子面
param[6], it contains the following data:
// [0],[1],[2] - the location of a rejected circumcent, 被拒绝的环的位置？？？啥东西
// [3] - the samllest edge length ( = insertion radius)
// [4] - ratio-edge ratio (of this subface). If it is zero, it is an encroached subface.
比率-边缘比率(此子面)。如果它是零，它是一个侵蚀子面。
// [5] - no used. 没用为啥要分配
*/
void repairencfacs(REAL *param, int qflag, int chkencflag);
/*
修复被侵占的子面
*/
bool check_tetrahedron(triface *chktet, REAL* param, int& qflag);
/*
检查一个四面体是否需要被拆分
"param[6]" returns the following data:
// [0],[1],[2] - the location of the new point 斯坦点的位置
// [3] - the samllest edge length ( = insertion radius) 插入半径
// [4] - the radius-edge ratio 半径边缘比
// [5] - (optional) edge ratio ？？？？
// "chktet" returns the shortest edge of this tet. 返回这个四面体的最短边
*/
bool checktet4split(triface *chktet, REAL* param, int& qflag);
/*
检查给定的tet是否有坏的形状。（应该是非常重要的）
*/
enum locateresult locate_point_walk(point searchpt, triface*, int chkencflag);
/*
通过直线搜索找到一个点。
*/
bool split_tetrahedron(triface*, REAL*, int, int, insertvertexflags &ivf);
/*
分割一个四面体
// param[6], it contains the following data
// [0],[1],[2] - the location of the new point
// [3] - the samllest edge length ( = insertion radius)
// [4] - radius-edge ratio
// [5] - its volume 四面体体积
// split due to mesh size enforcement. 由于网格尺寸强制而分裂。
*/
void repairbadtets(REAL queratio, int chkencflag);
/*
修复质量不好的四面体
*/
void delaunayrefinement();
/*
通过Delaunay细化网格。
*/
四、确定一个坏质量的四面体
确定坏的四面体可以从四面体体积，四面体顶点张量，四面体的外界球半径以及最短边比值，二面角大小确定。
tetgen中使用了两个版本的函数来判断一个四面体的质量，如下：
bool tetgenmesh::check_tetrahedron(triface *chktet, REAL* param, int &qflag)
bool tetgenmesh::checktet4split(triface *chktet, REAL* param, int& qflag)
前者可能是因为计算的精度不够，或者老旧的生命无法实现当下版本的新功能，因此已经被弃用，后者则是大量的使用了矩阵和向量的运算，但是两者思想相近，前者更加容易理解，因此建议先看 **check_tetrahedro（）**函数。
参数 *chktet是需要检查的四面体指针。
Real *param含有六个属性，分别是，分割该四面体时需要插入的新点的x，y，z坐标，四面体中的最短边的长度（插入半径），半径边缘比，第六个前一个函数代码中声明param[5] = vol；是指四面体的体积，后一个函数代码中声明**param[5] = sqrt(Lmax) / Lmin;**即最长边与最短边的比值
// "param[6]" returns the following data: //
// [0],[1],[2] - the location of the new point //
// [3] - the samllest edge length ( = insertion radius) //
// [4] - the radius-edge ratio //
// [5] - (optional) edge ratio //
五、size function的执行过程
关注-m，与-a开关
执行size function过程中调用的函数流程如下：
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-azlXWvRq-1669347053473)(C:\Users\25442\Desktop\未命名文件.jpg)]
第一步：在调用tetlib库的时候的接口函数中找到使用**-q（网格细化功能的函数，-m开关的执行包含在内）**部分的函数，即delaunayrefinement（）
第二步：delaunayrefinement（）包含两种模式，使用 - r和不适用 -Y 开关，我们要注意的是后面的部分，在这一部分里，会依次判断分割被侵占的段，被侵占面，和质量不好的四面体，并且分割它们，-a和-m开关使用的部分位于第三部分，这里需要依靠rapaiirebadtets（）函数来执行
第三步：在rapaiirebadtets（）函数对于每一个四面体都依靠checktet4split（）函数来检测是否需要分割，当返回结果为true的时候就调用split（）函数来完成分割，将分割的结果返回到repaire（）函数。
第二步：delaunayrefinement（）包含两种模式，使用 - r和不适用 -Y 开关，我们要注意的是后面的部分，在这一部分里，会依次判断分割被侵占的段，被侵占面，和质量不好的四面体，并且分割它们，-a和-m开关使用的部分位于第三部分，这里需要依靠rapaiirebadtets（）函数来执行
第三步：在rapaiirebadtets（）函数对于每一个四面体都依靠checktet4split（）函数来检测是否需要分割，当返回结果为true的时候就调用split（）函数来完成分割，将分割的结果返回到repaire（）函数。