Mzh's Studio.

三维表型--网格孔洞识别填补(一)————孔洞识别

几何, 三维重建, 植物
字数统计: 2.4k阅读时长: 11 min
2019/04/22 Share

半边排序-填补孔洞

记录如何得到半边的逻辑

根据delauney算法，若一条边最多被两个面共用，若一条边只被共用一次那么即为半边，只有半边可以形成孔洞，需要定义边类型，和半边类型

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
class myedge
{
public:
	int p1id;
	int p2id;
	int edge_times = 0;// the use times of this edge;
	int edge_valid = 0;// edge_valid = 0， edge is valid, otherwise  invalid;
	myedge();
	myedge(int p1, int p2)
	{
		if (p1 == p2) {
			cout << "Wrong edge!" << endl;
			edge_valid = 1;
		}
		else if (p1 > p2) {
			p1id = p2;
			p2id = p1;
		}
		else
		{
			p1id = p1;
			p2id = p2;
		}
	}
	bool operator == (const myedge &x)
	{
		if ((this->p1id == x.p1id) && (this->p2id) == x.p2id)
			return true;
		else
			return false;
	}
};

先三角化方法得到初步的网格数据vertics，然后转化为边类型
此处需要定义默认构造函数和重载==，因为边的两个点的顺序不确定，因此重载==时即考虑边相等时，两个情况都得考虑

注意：这里面半边数据结构集中的每个点最多使用两次，若要形成封闭的缺口，那么一定每个点使用两次

通过收入所有三角网格的边，并以端点的形式保存，v1<v2,然后计数每条边出现的次数，出现一次的为半边，其余的都为非半边删除

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
bool isSameEdge(myedge e1, myedge e2)
{
	if (e1.p1id == e2.p1id&&e1.p2id == e2.p2id)
		return true;
	else
		return false;
}

bool isInfront(myedge e1, myedge e2)
{
	if (e1.p1id == e2.p1id)
	{
		if (e1.p2id <= e2.p2id)
			return true;
		else
			return false;
	}
	else if (e1.p1id < e2.p1id)
		return true;
	else
		return false;
}

// get half edges，and save into vector
void getHalfEdge(PolygonMesh &mesh1, std::vector<myedge> &halfedges_out)
{
	// get halfedge;
	std::vector<myedge> halfedges;
	for (int i = 0; i < mesh1.polygons.size(); i++)
	{
		int p1 = mesh1.polygons.at(i).vertices[0];
		int p2 = mesh1.polygons.at(i).vertices[1];
		int p3 = mesh1.polygons.at(i).vertices[2];
		myedge e1(p1, p2);
		myedge e2(p2, p3);
		myedge e3(p3, p1);
		if (e1.edge_valid == 0)
			halfedges.push_back(e1);
		if (e2.edge_valid == 0)
			halfedges.push_back(e2);
		if (e3.edge_valid == 0)
			halfedges.push_back(e3);
	}
	std::vector<myedge>::iterator begin_pos = halfedges.begin();
	std::vector<myedge>::iterator end_pos = halfedges.end();
	std::sort(halfedges.begin(), halfedges.end(), isInfront);
	// above process is to get all edges(includeing halfedges and whole edges) and take them in order
	// if the edgd showed twice, it is not the halfedges

	std::vector<myedge> halfedges_copy;
	halfedges_copy = halfedges;

	std::vector<myedge>::iterator new_end = std::unique(halfedges.begin(), halfedges.end(), isSameEdge);
	halfedges.erase(new_end, end_pos);

	// Counts
	for (int i = 0; i < halfedges.size() - 1; i++)
	{
		auto begin_temp =
			find(halfedges_copy.begin(), halfedges_copy.end(), halfedges.at(i));
		auto end_temp =
			find(halfedges_copy.begin(), halfedges_copy.end(), halfedges.at(i + 1));
		halfedges.at(i).edge_times = end_temp - begin_temp;
	}

	// get halfedges
	for (int i = 0; i < halfedges.size() - 1; i++)
	{
		if (halfedges.at(i).edge_times == 1)
			halfedges_out.push_back(halfedges.at(i));
	}
	cout << halfedges_out.size() << endl;
}

检索孔洞逻辑

步骤：

1、通过邻接表保存半边

2、通过递归的方式剔除半边中的端点，即不能形成闭环的点，即只被一条边共用的端点；

3、识别孔洞

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
int recursion_times = 0;
void DeleteEndPoint(Lgraph gp, int searchsize)
{
	cout << "Recursion times is: " << recursion_times << endl;
	std::vector<int> deleteID;
	// 由于find的得到的是第几个节点，因此得到n还需+1，因此链表下标从0开始
	for (int i1 = 0; i1 < searchsize; i1++)
	{
		int len = getlen(gp->G[i1]);
		if (len==1)
		{
			deleteID.push_back(i1);
		}
	}
	if (deleteID.size() != 0)
	{
		for (int i = 0; i < deleteID.size(); i++)
		{
			for (int j = 5485; j < searchsize; j++)
			{
				int len = getlen(gp->G[j]);//当邻接表表长为0的时候说明是空表，
				if(len==0)
					continue;
				auto temp = delete_Vertex(deleteID.at(i), gp->G[j].FirstEdge);
				gp->G[j].FirstEdge = temp;
				gp->G[deleteID.at(i)].FirstEdge = NULL;
				/*if (gp->G[j].FirstEdge->Next == NULL)
					gp->G[j].FirstEdge = NULL;*/
			}
		}
		recursion_times++;
		DeleteEndPoint(gp, searchsize);
		
	}
	else
		return;

}

void getHolesmap(std::vector<myedge> &halfedges, Lgraph &ph, std::vector<std::vector<int>> &holemap2)
{
    // save sa adjacent list
    std::vector<int>pts;
	for (int i = 0; i < halfedges.size(); i++)
	{
		pts.push_back(halfedges.at(i).p1id);
		pts.push_back(halfedges.at(i).p2id);
	}
	sort(pts.begin(), pts.end());
	pts.erase(unique(pts.begin(), pts.end()), pts.end());
	
	int searchSize = FindMaxValue(pts);

	Lgraph graph;
	graph = CreatGraph(pts.size());
	graph->Ne = halfedges.size();
	if (graph->Ne != 0)
	{
		for (int i = 0; i < graph->Ne; i++)
		{
			InsertEdge(graph, halfedges.at(i));
		}
	}

    // delete end point using recursion
	DeleteEndPoint(graph, searchSize);
	ph = graph;

	std::vector<int> lens;
	int len4 = 0;
	int len6 = 0;
	for (int i = 0; i < searchSize; i++)
	{
		int lentep = getlen(ph->G[i]);
		if (lentep != 0)
		{
			lens.push_back(lentep);
		}
		else
			continue;;
		if (lentep == 4)
		{
			len4++;
		}
		if (lentep == 6)
		{
			len6++;
		}
	}
	cout << "Number of points have 4 adjacent points: " << len4 << endl;
	cout << "Number of points have 6 adjacent points: " << len6 << endl;
}

在得到孔洞的半边数据后，由于有些点被共用多次，如下图所示

hole_situation

  1. 情况1表示孔洞所有点只被共用2次；
  2. 情况2表示孔洞红色点被共用4次，其余点被共用2次；
  3. 情况3表示孔洞红色点被共用6次，其余点被共用2次；
  4. 情况4表示多种情况混合，即图中红色点分别被共用6,4次（也可能存在4（4,4），（6,6）等情况的发生），其余点被共用2次；

因此在采用深度优先搜索算法时需要进行分情况说明，伪码描述为：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
vector<vector<node>> holes;
DFS_hole(vertex v,vector<node> out)//此处仅考虑len=2的情况
{
    int len=getlen(V); //获取节点后面的长度，长度为2,4,6
    if(len==2)//若这个节点长度是2，说明该节点被4条边共用，属于情况1
    {
        if(v.visit==0)//visit=0说明该节点还未访问过，添加到输出vector out中
        {
            vector.push_back(v);
            v.visit++;
            v=v.nest;
            DFS_hole(v,out);
        }
        //visit=1说明该节点已经被访问，且又回到该节点，那说明该节点是起点且又回到起点，则out里一定有一个闭环，则holes中添加该闭环，并退出梯柜
        else if(v.visit==1)
        {
            holes.push_back(out);
            return;
        }
        //如果访问≥1次，在len=2时不存在这种情况，所以直接return；
        //不存在的原因时候len=4或len=6时会在生成一个闭环时终止递归，因此对于所有len=2的点最多只可能变量2次
        else
            return；
    }
}

上文仅考虑len=2时的情况，即情况1的情况，然而实际情况中是存在2,3,4三种情况，因此在执行上述深度优先搜索算法提取空洞前，先考虑将这些复杂情况的孔洞情况删除，即从邻接表中删除这些节点并保存至另一复杂孔洞的领接表中，伪码描述为：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
//形参：检索起始节点，输入图
void dfs_complex_holes(vnode v, Lgraph listin)
{
    //因为进入这个递归过程中的点，都是复杂情况下的闭环的点，因此标记这些点为复杂闭环圈；
    //标记完后递归
    adjacentlist w;
    v.isComplex=1;
    for(w=v.firstedge;w;w=w->next)
    {
        if(w.next.isComplex==0)
        {
            dfs_complex_holes(listin->g[w.vertex],listin);
        }
    }
}
void get_complex_holes(Lgraph listin，int searchtimes，Lgraph listout)
{
    
    for（int i=0;i<searchtimes;i++)
    {
      int len=getlen(listin->v[i]);
        if(len==0)
            continue;
        else
        {
            //将遍历点都从复杂点开始，即被共用次数为4,6等，除0外
            if(len!=2)
                dfs_complex_holes(listin->g[i],listin);
        }
    }
     
}

项目代码，可运行Ver1.0：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
//标记含有复杂点的闭环圈，即圈内存在被共用次数＞2的点
void dfs_complex_holes(Vnode v, Lgraph listin)
{
	//因为进入这个递归过程中的点，都是复杂情况下的闭环的点，因此标记这些点为复杂闭环圈；
	//标记完后递归
	PtrToAdjVnode w;
	if (v.isComplex==0)
	{
		listin->G[v.id].isComplex = 1;// v.isComplex = 1;
	}
	else
	{
		return;
	}
	for (w = v.FirstEdge; w!=NULL; )
	{
		auto temp = listin->G[w->AdjV];
		if ((listin->G[w->AdjV].isComplex == 0))
		{
			dfs_complex_holes(listin->G[w->AdjV], listin);
			if (w->Next == NULL)
			{
				break;
			}
			else
			{
				w = w->Next;
			}
		}
		else
		{
			if (w->Next == NULL)
			{
				break;
			}
			else
			{
				w = w->Next;
				continue;
			}
		}
	}
}
//获取复杂点
void get_complex_holes(Lgraph listin, int searchtimes)
{

	for(int i = 0; i < searchtimes; i++)
	{
		int len = getlen(listin->G[i]);
		if (len == 0)
			continue;
		else
		{
			//将遍历点都从复杂点开始，即被共用次数为4,6等，除0外
			if (len != 2)
				dfs_complex_holes(listin->G[i], listin);
		}
	}
}
//深度优先搜索检索孔洞
void DFS_hole(Vnode v, Lgraph listin, std::vector<Vnode> &out,
	std::vector<std::vector<Vnode>> &holes)//此处仅考虑len=2的情况
{
	//如果该节点是复杂点，则不进入搜索闭环的环节；
	if (v.isComplex==1)
	{
		return;
	}
	int len = getlen(v); //获取节点后面的长度，长度为2,4,6
	if (len == 2)//若这个节点长度是2，说明该节点被4条边共用，属于情况1
	{
		PtrToAdjVnode w;
		if (v.Visited==0)
		{
			listin->G[v.id].Visited = 1;
			out.push_back(v);
		}
		for (w = v.FirstEdge; w != NULL; )
		{
			if (listin->G[w->AdjV].Visited == 0)
			{
				DFS_hole(listin->G[w->AdjV], listin, out, holes);
			}
			if (w->Next == NULL)
			{
				break;
			}
			else
			{
				w = w->Next;
			}
		}
	}
}
//获取正常孔洞
void get_holes(Lgraph ph, std::vector<std::vector<Vnode>> &holes)
{
	for (int i = 0; i < MaxVertexuUm; i++)
	{
		std::vector<Vnode> temp;
		DFS_hole(ph->G[i], ph, temp, holes);
		if (temp.size()!=0)
		{
			holes.push_back(temp);
		}
	}
}

识别效果如下：对于正常孔洞已经能较好的识别，但是存在将最外边界的点即形成封闭的叶片轮廓的点也识别为孔洞；同时还有一些其他的面片也易被识别为孔洞，后续需要根据几何性质进行情况分析。

hole3

下一章主要为填补孔洞

豪斯多夫距离Hausdorff distance of convex polygons这部分没看完，先空着，与评估方法效果有关

原文作者：MZH

原文链接：https://singlemzh.github.io/2019/04/22/Note1/

发表日期：April 22nd 2019, 10:12:33 am

更新日期：October 29th 2020, 9:37:55 am

版权声明：本文采用知识共享署名-非商业性使用 4.0 国际许可协议进行许可

Total : 16
2021
2020
2019
2018
Polygon, Mesh, Processing 点云处理 股票 Deep learning, AI 几何, 三维重建, 植物 三维重建, 植物 生活，思考 论文写作, 汇报注意 几何, 三维重建
MeshProcess HandlePoints 股市 DeepLearning 3DReconstruction Others Life Writing 理财