利用CGAL库函数细分约束三角剖分以及搜寻约束边界内的三角形面(CDT refinement in CGAL)

最近自己要研究一个约束三角剖分的细分功能（对一个服装进行三角剖分、细分并且变形），搜了很久，找到一个比较对口的库——CGAL库。但是国内使用的人很少，资料也很少。所以自己还是处于最基本的摸索阶段。虽然自己也是属于一桶水连一瓢都木有的状态（跟不用说半桶水了orz...），不过还是经过反复试验，实现了自己所需要的功能（中间有一些小问题还尚未解决，写出来希望大牛们一起探讨这么可以解决）。

问题：如下图，提取服装的轮廓之后，对其进行约束三角剖分（轮廓作为约束边）并且细分。



提取的服装轮廓（用关键点的连线多边形表示）如下图：



将以上由关键点组成的polygon作为边界约束，进行三角剖分以及细化。为了方便，首先加入一些头文件并且定义一些CGAL数据（其实我也不是很懂这些参数的意义，反正看到CGAL用户手册上总是定义一堆数据，下面的这些定义都是差不多，我直接copy过来的。你们也照着抄吧，哈哈哈）。原版CGAL实例中定义CDT::Point Point，但是因为我的应用中还用到openCV，里面也有一个cv::Point，为了防止冲突我就改成了typedef CDT::Point
CdtPoint.

然后对以上的这些关键点组成的封闭多边形进行约束三角剖分(CDT: Constrained Delaunay Triangulation)，所以先在主函数定义一个CDT cdt; 然后将这些点将入CDT中，并且每对相邻点是一个边约束。

主要用到的函数：

<span style="font-size:14px;">Vertex_handle insert(const Point& p, Face_handle start = Face_handle());
Constrained_triangulation_2<Gt,Tds,Itag>::insert_constraint(Vertex_handle  vaa, Vertex_handle vbb);</span>

我直接写了一个InsertPolygon()的函数：

<span style="font-size:14px;">void InsertPolygonInCDT(CDT& cdt, vector<CdtPoint> polygon)
{
int s = polygon.size();
for (int i=0; i<s; i++)
{
Vertex_handle v_pt = cdt.insert(polygon[i]);
Vertex_handle v_pt2 = cdt.insert(polygon[(i+1)%s]);//依次插入两个相邻的关键点

cdt.insert_constraint(v_pt, v_pt2);//将相邻点作为一个边约束插入CDT
}
}</span>

接着，对cdt进行细分，用到的函数为：

<span style="font-size:14px;">CGAL::refine_Delaunay_mesh_2(CDT cdt, Criteria(0.125, min_len));</span>

其中Criteria指示的是细分度，第一个参数默认为0.125（具体含义看手册好像是和B有关的细分参数，其余我也不是很清楚，索性都直接写0.125好了），第二个参数表示细分中最大的边长长度，用户可以自己指定。

至此，主函数前面为：

<span style="font-size:14px;">int main()
{
vector<CdtPoint> polygon;
polygon.push_back(CdtPoint(...,...));
......
CDT cdt;
InsertPolygonInCDT(cdt, polygon);
CGAL::refine_Delaunay_mesh_2(cdt,Criteria(0.125, min_len));
}</span>

显示细分三角剖分，可以用遍历CDT的方法（如何访问CDT的边以及顶点数据，具体函数可参考我的前一篇博文），这里我自定义了两个函数，用来CDT图的描绘：

<span style="font-size:14px;">void DrawEdge(Mat& img, Edge e, cv::Scalar color)
{
Face_handle fh = e.first;
int e_id = e.second;

Vertex_handle vh1 = fh->vertex(fh->ccw(e_id));
Vertex_handle vh2 = fh->vertex(fh->cw(e_id));
DrawEdge(img, vh1, vh2, color);
}

void DrawEdgesInCDT(cv::Mat& img, Scalar color, CDT& cdt)
{
CDT::Finite_edges_iterator eit;
for (eit=cdt.finite_edges_begin(); eit!=cdt.finite_edges_end(); eit++)
{
DrawEdge(img, *eit, color);
}
}</span>

最后显示如下：



从结果来看，CGAL::refine_Delaunay_mesh_2()可以得到很好的细分结果，但是它的结果是一个“凸包”，不能正确地将边界约束当作真正的三角剖分内部剖分约束（或者是我水平太low，还没有找到相关的函数）。而我真正需要的是以下的这种结果：



这部分功能实现只能自己想了。于是我用类似于“漫水法”的思想，找到包含在约束边界内的全部三角面/边/点。为了说明自己的实现思想，我们先来考虑另一个问题：如何在一个多边形边界内找到内部的所有像素点？

基本思路就是：从多边形的内部找到一个像素点，然后不断搜索其周围的点进行扩大，如果遇到了边界，就停止生长。最后由一个点开始，将多边形中的所有像素点都占满。

具体步骤：

1. 将多边形所有的边界所占的像素点值设为1（可以直接用openCV中的画线函数进行设置），图像其余点的像素值都为0。

2. 然后，在多边形内部找到一个初始的像素点，将其设为1，并作为第一个种子点，压入栈。

3. 每次在栈中取出一个种子点，将其所在的像素位置值设为1，并且搜寻这个种子点的（上、下、左、右）4个相邻像素点，如果相邻像素点为0，则它就成为了一个新的种子点，压入栈；如果为1（表示已经生长的像素点，或者是遇到了边界像素点），舍弃。

4. 重复步骤3，直至栈中的种子点被访问结束。

类似于以上的思想，搜索约束边界内的所有三角面/边，步骤如下：

1. 在细分的三角剖分图（CDT Graph）中找到一个内部的三角面

2. 将这个三角作为种子三角面压入栈，并任意选取其中一条边作为起始种子边压入栈（共有两个栈）；另外用一个容器存储已经搜寻的三角面。

3. 从栈中分别取出种子三角面和种子边（两者是对应的，该边是该三角面的一条边），用该边去搜索该三角面的另外两条边：如果其余两条边不是约束边界边（遇到了边界），则搜索边的相邻面，并且判断该两边和对应的相邻三角面是否已经被访问过（即是否存在于容器中），如果没有被访问过，则将边和相邻面作为种子压入栈；否则停止增加。

4.重复以上步骤3，直至栈中的种子边和种子三角面为空。

关键代码：

在步骤1. 用CDT::is_constrained(e)来判断边是不是约束边，如果三条边都是非约束边，说明是一个内部三角面，如下

<span style="font-size:14px;"><span style="font-size:14px;"><pre name="code" class="cpp">for (e_iter=cdt.finite_edges_begin(); e_iter!=cdt.finite_edges_end(); e_iter++)
{
Face_handle fhd = e_iter->first;
int e_idx = e_iter->second;

int con_count = 0;
for (int i=0; i<3; i++)
{
Edge e(fhd, i);
if(cdt.is_constrained(e))//
con_count++;
}

if(con_count==0)//三条边都是非约束边
{
...//第2步的主要代码
}
<span style="white-space:pre">	</span>break;
}</span></span>

在步骤2. 关键代码，edges_in和edges_con是std::set<Edge>类型，分别表示搜索到的所有内部边和边界边。这两个容器是作为参数传入的，所以没有写creation function

<span style="font-size:14px;"><span style="white-space:pre">			</span>//从获取的第一个种子边开始</span><span style="font-size:14px; font-family: Arial, Helvetica, sans-serif;">			</span>

<span style="font-size:14px;"><span style="white-space:pre">			edges_in.insert(*e_iter);</span>
std::set<Face_handle> fset;//用来判断存储已经访问过的三角面
std::set<Vertex_handle> vset; <span style="font-family: Arial, Helvetica, sans-serif;">//用来存储顶点</span>

std::list<Face_handle> seed_faces;//种子三角面的栈
std::list<Vertex_handle> seed_vertices;//种子边的栈，用两个顶点vertex_handle来替代

//压入第一个seed face
Face_handle face_start = e_iter->first;
int e_idx = e_iter->second;
seed_faces.push_back(face_start);

//第一个三角面已被访问
fset.insert(face_start);
vset.insert(face_start->vertex(0));
vset.insert(face_start->vertex(1));
vset.insert(face_start->vertex(2));

//压入第一个seed edge的两个顶点
Vertex_handle vh1 = face_start->vertex(face_start->ccw(e_idx));
Vertex_handle vh2 = face_start->vertex(face_start->cw(e_idx));
seed_vertices.push_back(vh1);
seed_vertices.push_back(vh2);</span>

在步骤3中，具体用到的CGAL库函数，主要是：种子边e，用e.first找到它所在的三角面，用e.second指示这条边在这个三角面的index，然后用这个index可以找到两边的index，从而得到另两条边的信息。如下：

<span style="font-size:14px;">		//从种子三角面开始，不断搜索周围的面
while(!seed_faces.empty())
{
count++;
//取出第一个种子面
Face_handle fh = seed_faces.front();
seed_faces.pop_front();

//取出第一个种子边的两个顶点
Vertex_handle vh1 = seed_vertices.front();
seed_vertices.pop_front();

Vertex_handle vh2 = seed_vertices.front();
seed_vertices.pop_front();

//得到两个顶点在该种子面的index,从而计算得到第三个顶点的index
int vh1_idx = fh->index(vh1);
int vh2_idx = fh->index(vh2);
int vh3_idx = (fh->ccw(vh1_idx)==vh2_idx)? fh->cw(vh1_idx):fh->ccw(vh1_idx);

//得到第三个顶点
Vertex_handle vh3 = fh->vertex(vh3_idx);

if(vset.find(vh3)!=vset.end());
else
vset.insert(vh3);

//除了取出的种子边之外，如果另两条边未被访问，则存储
Edge e1(fh, vh1_idx);
Edge e2(fh, vh2_idx);
if(edges_in.find(e1)!=edges_in.end());
else
{
edges_in.insert(e1);
}

if(edges_in.find(e2)!=edges_in.end())
edge_repeat++;
else
{
edge_count++;
edges_in.insert(e2);
}

//访问另两边，若是constrained, 表示达到了边界, 不再访问相邻三角面
if(!cdt.is_constrained(e1))
{
Face_handle nf1 = fh->neighbor(vh1_idx);//e1的相邻面

if(fset.find(nf1)!=fset.end());
else//未访问,压入种子面和对应的种子边
{
fset.insert(nf1);
seed_faces.push_back(nf1);
seed_vertices.push_back(vh2);
seed_vertices.push_back(vh3);
}
}
else
{
edges_con.insert(e1);
}

if(!cdt.is_constrained(e2))
{
Face_handle nf2 = fh->neighbor(vh2_idx);

if(fset.find(nf2)!=fset.end());
else
{
fset.insert(nf2);
seed_faces.push_back(nf2);
seed_vertices.push_back(vh1);
seed_vertices.push_back(vh3);
}
}
else
{
edges_con.insert(e2);
}
}
</span>

最后得到的结果如下: 分别用绿色和红色表示了约束边界的内部边和边界边



不过以上的代码中有一些问题：用set储存Face_handle和Vertex_handle类型的三角面和边，都没有什么问题。但是在存储Edge类型的时候，因为Edge的定义是一个二元组<Face_handle, int>, 第一个是它所在的一个面，第二个是它所在该面的index，以上步骤3中在判断该边是否已经存在的时候（edges_in.find(e1)!=edges_in.end()）其实是一个无效函数，这是因为每次在存储一条边的时候，它对应的Face_handle,
int组合和之前是不同的，也就是说，同一条边它会被相邻三角面经过两次储存（尽管它实际上指示的是同一条边）。因此，以上代码最后得到的结果edges_in内部边是有重复性的。（比较笨的方法就是将得到所有边进行两两对比判断，去除重复的数据。但是我觉得这个方法是后续的补救方法，我想要的是在内部实现的时候就能够解决这个问题。）这个哪位大牛来帮忙解决一下？

另外，还有一个问题就是内部三角形的搜索问题。以上这个实例中，因为除了服装内部，约束边界外正好没有一个三边都是非约束边的三角面，所以运行正确。但是如果衣服的手臂是张开的（空隙比较大），那么，在约束边界外部可能也存在这样符合条件的“内部三角面”，这样得到是边界外自成独立的一块三角面集合。这个问题我也没有想到很好的解决方法，请各位赐教！