前排瞌睡杂物堆

LKH 求存，探路，躬行，旁通，振兴，延续。
第二阶段，启动！

0%

点云三维重建 PolyFit算法解析与改进（六）

发表于 更新于 分类于
本文字数： 33k 阅读时长 ≈ 1:01

摘要：
平面相交细分；候选面集简化；面片置信度计算；候选平面集。

中场休息结束！！ok，继续分析PolyFit中构造面片的代码。
上一章节，建筑物点云被处理为表示各个区域朝向的面片数据，存储在多边形网络mesh_candidate中。尽管这些提取到的面对象能初步表示建筑物的几何特性，但显然它们之间是孤立的，缺少相交、邻接、连通等拓扑关系，且明显有冗余的部分区域需要被剪除。那么，如上一节末尾所述，接下来便是对网络中的各个面片做细化、去冗处理。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
//4、候选平面相交细分
//2023.11.23 在包围盒平面与细分的平面上，对面的属性增加分区属性
SegMeshCandidate(mesh_candidate);


//5.1、候选面集简化处理
//SimpleMeshCandidate(mesh_candidate);//激进简化，有时重建效果巨差
SimpleMeshCandidate_V2(mesh_candidate);//保守简化

//5.2、附加顶面偏好属性
AddMeshProperty(mesh_candidate);

//6、计算候选面的置信度
//该操作为 mesh_candidate 添加了3个属性映射,用于后续的MIP处理：
//平面的支撑点数	f:num_supporting_points
//平面面积			f:face_area
//覆盖面积			f:covered_area
CGAL::internal::Candidate_confidences<Kernel> conf;
conf.compute(cloud_plane, mesh_candidate);

候选平面相交细分

现在的mesh_candidate中存有若干块表示建筑特征区域朝向的大块面片，若直接以聚类区域边界划分面片，由于点云数据的噪声及区域表面的不平滑特性（可能表征为多个不同方向的平面），该过程往往会切割出一些非预期的平面；且点云覆盖度若较差，容易在模型中产生空洞或不连续的区域，从而降低重建模型的几何精度和视觉质量。而按照这种思路切分出的各区域面，还需要进行面片间的复杂拓扑关系处理后才能整合成模型表面，如相交面的边界融合/平滑、形状规则化[20]、同名相交点的精度控制。
PolyFit算法中的平面细化采用的是暴力穷举的思路。不同于上述那种先取轮廓再不断优化边界的策略，它是先以平面两两相交的方式将mesh_candidate的大块面切成超多细碎面，从而显式表达了聚类区域的所有相交情况；再按照各个面的能量函数值删除、保留面片，最终将隐含的模型表面给筛选出来。这通操作避免了对模型相交边界的优化问题，将模型重建方法从面片结构和面片间的拓扑关系的优化问题，转为从潜在的多边形面集内筛选出组成模型的最优面集。当然缺点也是有的：两两相交会产生大量冗余面，会导致在计算能量项那块出现计算瓶颈。
这种取切面策略的缺陷后面再详述，先聚焦回平面相交细分，它的代码在HypothesisPlaneProcess::SegMeshCandidate()函数内，过程包括计算潜在的三重交点、确定相交面、平面分割这三个步骤。

计算三重交点

按顺序取mesh_candidate中的三个平面，使用CGAL::intersection计算它们三个间是否交于一点，有则将点存入_triplet_intersections。老实说这一步放在这里就计算有点太早了，因为这些点的信息很后面才用到，得到了网络结构重建的ReconstructProcess::GetAdjacency函数中才会调取它们来计算面的邻接情况。

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
//1、预先计算平面三元组（依次选取相邻顺序的三个平面）所有潜在的三重交集点（三平面交于一点的点集）
std::vector<const Kernel::Point_3*> intersecting_points_;//容器_三重交点//8.29_似乎仅做统计，后续用不上
_triplet_intersections.clear();
if (_supporting_planes.size() < 4)
	return;//少于4个平面将无法构成闭合网络，退出本函数

for (std::size_t i = 0; i < _supporting_planes.size(); ++i)
{
	const Kernel::Plane_3* plane1 = _supporting_planes[i];
	for (std::size_t j = i + 1; j < _supporting_planes.size(); ++j)
	{
		const Kernel::Plane_3* plane2 = _supporting_planes[j];
		for (std::size_t k = j + 1; k < _supporting_planes.size(); ++k)
		{
			const Kernel::Plane_3* plane3 = _supporting_planes[k];
			CGAL_assertion(plane1 < plane2 && plane2 < plane3);//断言_三个平面的指针是否按升序排列的

			if (plane1 == plane2 || plane1 == plane3 || plane2 == plane3)
				continue;//任意两个平面相同，则不会有三重交点，跳过此循环

			//计算三个平面的交集，如果交集的对象可以转换为一个Kernel::Point_3类型的指针，则三个平面相交于一点
			CGAL::Object obj = CGAL::intersection(*plane1, *plane2, *plane3);
			if (const Kernel::Point_3* pt = CGAL::object_cast<Kernel::Point_3>(&obj))
			{
				//复制三重交点的坐标，存储点
				Kernel::Point_3* new_point = new Kernel::Point_3(*pt);
				_triplet_intersections[plane1][plane2][plane3] = new_point;
				intersecting_points_.push_back(new_point);
			}
			else
				continue;//未得到交点的两种原因：（1）面是平行的；（2）面相交于同一条线上。

		}
	}
}

确定相交面（及其交点）

这一小节的目的是确定mesh_candidate内某个大块面与哪些大块面存在相交关系。遍历每个面，取出当前面的下标index_face_i及其平面方程plane_face_i。再套一个同样的遍历取候选切面f，通过函数HypothesisPlaneProcess::compute_intersections()计算两个面的交点数量，根据相交点的情况将相交面存入intersecting_faces。
为便于回顾，下述代码省略了一些细枝末节的判断和属性别名。

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
std::vector<Polygon_mesh::Face_index> all_faces(mesh_candidate.faces().begin(), mesh_candidate.faces().end());//容器_所有候选面
for (std::size_t i = 0; i < all_faces.size(); ++i)
{
	//获取第i个面的索引和所在平面
	Polygon_mesh::Face_index index_face_i = all_faces[i];
	const Kernel::Plane_3* plane_face_i = face_supporting_planes[index_face_i];

	//获取相应的属性映射
	typename Polygon_mesh::template Property_map<Polygon_mesh::Face_index, const Kernel::Plane_3*> face_supporting_planes =
		mesh_candidate.template property_map<Polygon_mesh::Face_index, const Kernel::Plane_3*>("f:supp_plane").first;//属性映射_每个面的支撑平面

	//2.1、遍历所有候选面，将与第i个面相交的面存入集合 intersecting_faces
	std::set<Polygon_mesh::Face_index> intersecting_faces;
	for (auto f : mesh_candidate.faces())
	{
		//2.2、判断面 f 是否与平面 plane_cutting 相交
		const Kernel::Plane_3* plane_cutting = face_supporting_planes[index_face_i];//获取给定面作为切割面
		CGAL_assertion(plane_cutting != nullptr);//检查平面是否为空指针

		std::vector<Polygon_mesh::Vertex_index> existing_vts;//平面相交的已有顶点下标
		std::vector<EdgePos> new_vts;//平面相交的新生成的顶点
		compute_intersections(mesh_candidate, f, plane_cutting, existing_vts, new_vts);//计算交点

		bool is_intersect_fPlane = false;
		if (existing_vts.size() == 2)
		{
			if (!halfedge_exists(existing_vts[0], existing_vts[1], mesh_candidate))
				is_intersect_fPlane = true;//2个已有的顶点与plane_cutting相交且两点不在同一条边上
		}
		else if (existing_vts.size() + new_vts.size() == 2)
			is_intersect_fPlane = true;//已有+新生成的共两个点与plane_cutting相交

		if (is_intersect_fPlane)
			intersecting_faces.insert(f);//如果f和第i个面所在的平面相交,将f加入到intersecting_faces

	}

	if (intersecting_faces.empty())
		continue;//对第i个候选面，没有与之相交的候选面，则跳过
	std::vector<Polygon_mesh::Face_index> cutting_faces(intersecting_faces.begin(), intersecting_faces.end());//容器_与第i个候选面相交的其他候选面

	//3、使用与 index_face_i 相交的平面分割 index_face_i
	//......
}

代码中创造了两个变量，存储mesh_candidate已有点下标的existing_vts和存储两平面相交产生的新点的对象new_vts。由于它们的值取决于compute_intersections()，而后面相交关系的if判断也是基于这两变量得到的，不难发现compute_intersections()的交点计算过程才是这小节的重点。来看看compute_intersections函数中对这这两个变量的处理逻辑。

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
void HypothesisPlaneProcess::compute_intersections(const Polygon_mesh & mesh, Polygon_mesh::Face_index face, const Kernel::Plane_3 * plane_cutting, std::vector<Polygon_mesh::Vertex_index>& existing_vts, std::vector<EdgePos>& new_vts)
{
	existing_vts.clear();
	new_vts.clear();

	typename Polygon_mesh::template Property_map<Polygon_mesh::Face_index, const Kernel::Plane_3*> face_supporting_planes =
		mesh.template property_map<Polygon_mesh::Face_index, const Kernel::Plane_3*>("f:supp_plane").first;//属性映射_每个给定面所在的支撑面
	const Kernel::Plane_3* supporting_plane = face_supporting_planes[face];
	if (supporting_plane == plane_cutting)
		return;//给定面与切割平面相同，则不需计算交点，退出函数

	typename Polygon_mesh::template Property_map<Polygon_mesh::Edge_index, std::set<const Kernel::Plane_3*> > edge_supporting_planes
		= mesh.template property_map<Polygon_mesh::Edge_index, std::set<const Kernel::Plane_3*> >("e:supp_plane").first;//属性映射_每条边所在平面的集合


	Polygon_mesh::Halfedge_index cur = mesh.halfedge(face);//给定面的第一个半边的索引
	Polygon_mesh::Halfedge_index end = cur;//记录 cur ，用于结束循环
	do {
		//获取半边和其所在的平面
		Polygon_mesh::Edge_index ed = mesh.edge(cur);
		const std::set<const Kernel::Plane_3*>& supporting_planes = edge_supporting_planes[ed];
		if (supporting_planes.find(plane_cutting) != supporting_planes.end())
			return;//边的平面集合中包含切割平面，则说明当前边在切割平面上，退出函数

		//获取边的两个顶点及坐标
		Polygon_mesh::Vertex_index s_vd = mesh.source(cur);
		Polygon_mesh::Vertex_index t_vd = mesh.target(cur);
		const Kernel::Point_3& s = mesh.points()[s_vd];
		const Kernel::Point_3& t = mesh.points()[t_vd];

		//判断两个顶点相对于切割平面的位置
		CGAL::Oriented_side s_side = plane_cutting->oriented_side(s);
		CGAL::Oriented_side t_side = plane_cutting->oriented_side(t);
		if (t_side == CGAL::ON_ORIENTED_BOUNDARY)//目标顶点在plane_cutting上
		{
			if (s_side == CGAL::ON_ORIENTED_BOUNDARY)//源顶点在plane_cutting上
				return;//边在plane_cutting上，不需计算交点，退出函数
			else
				existing_vts.push_back(t_vd);//目标点为交点，添加到已有顶点集
		}
		else
			if ((s_side == CGAL::ON_POSITIVE_SIDE && t_side == CGAL::ON_NEGATIVE_SIDE) ||
				(s_side == CGAL::ON_NEGATIVE_SIDE && t_side == CGAL::ON_POSITIVE_SIDE))//源顶点和目标顶点在plane_cutting的两侧
			{
				//当前边与plane_cutting在内部相交
				//计算两顶点到plane_cutting的距离的平方
				FT s_sdist = CGAL::squared_distance(*plane_cutting, s);
				FT t_sdist = CGAL::squared_distance(*plane_cutting, t);

				if (s_sdist <= CGAL::snap_squared_distance_threshold<FT>())//源顶点到plane_cutting的距离小于阈值，说明plane_cutting在源顶点处切割
					existing_vts.push_back(s_vd);//源顶点添加到已有顶点集
				else
					if (t_sdist <= CGAL::snap_squared_distance_threshold<FT>())//目标顶点到plane_cutting的距离小于阈值
						existing_vts.push_back(t_vd);//目标顶点添加到已有顶点集
					else
					{
						const Kernel::Plane_3* plane1 = *(supporting_planes.begin());//边所在平面中的第一个平面
						const Kernel::Plane_3* plane2 = *(supporting_planes.rbegin());//边所在平面中的最后一个平面
						const Kernel::Plane_3* plane3 = const_cast<const Kernel::Plane_3*>(plane_cutting);//获取plane_cutting

						//如果切割平面和当前边所在的两个平面不同，则三个平面可构成三重交点
						if (plane3 != plane1 && plane3 != plane2)
						{
							sort_increasing(plane1, plane2, plane3);//三个平面按指针升序的顺序排序，便于查询
							//const Kernel::Point_3* p = query_intersection(plane1, plane2, plane3);//查询三个平面的交点
							//2023.11.5 修改——————
							//计算三个平面的交集，如果交集的对象可以转换为一个Kernel::Point_3类型的指针，则三个平面相交于一点
							CGAL::Object obj = CGAL::intersection(*plane1, *plane2, *plane3);
							Kernel::Point_3* p;
							if (const Kernel::Point_3* pt = CGAL::object_cast<Kernel::Point_3>(&obj))
							{
								//复制三重交点的坐标，存储点
								p = new Kernel::Point_3(*pt);
								_triplet_intersections[plane1][plane2][plane3] = p;
							}
							else
								p = nullptr;
							//————————————
							if (p)//查询到有交点
							{
								if (CGAL::squared_distance(*p, s) <= CGAL::snap_squared_distance_threshold<FT>())//交点和源顶点的距离小于阈值，交点视为源顶点
									existing_vts.push_back(s_vd);
								else if (CGAL::squared_distance(*p, t) <= CGAL::snap_squared_distance_threshold<FT>())//交点和目标顶点的距离小于阈值，视为目标顶点
									existing_vts.push_back(t_vd);//点存入已有顶点集
								else
									new_vts.push_back(EdgePos(ed, p));//将交点与当前边构成EdgePos对象，加入到新生成的顶点集
							}
							else
								std::cerr << "Fatal error: should have intersection" << std::endl;//错误_不存在交点
						}
						else
							std::cerr << "Fatal error: should not have duplicated planes." << std::endl;//错误_出现重复平面
					}
			}
			else
			{
				// Nothing needs to do here, we will test the next edge
			}


		cur = mesh.next(cur);

	} while (cur != end);

}

师弟SY曾与我探讨了这块的内容（感谢.jpg），一并总结下：
先根据当前面在mesh_candidate的下标，循环取该面的半边结构cur，得到该半边的两个顶点s和t。将这两个点与（相交）候选切面的平面方程plane_cutting做位置上的判断，有这么几种情况：
（1）边的两个顶点都在plane_cutting上，则当前面与候选切面刚好相交于一条线上，直接retrun不需要该候选切面来给当前面做切割；
（2）仅有一个点t在plane_cutting上，将这个点记录到existing_vts中。不对点s做这个判断是为了避免在别的半边那里重复记录，不理解的话跳去SegMeshCandidate()中回顾2.2两个面的相交判断。若后续还能记录到另一个点到existing_vts/new_vts，则这两点的连线为切线。
（3）s和t分别在plane_cutting两侧。先分别对两点做了个到候选切面的距离判断，相当于更宽松条件下的（2）情况。若s与t离候选切面较远，那取出生成当前边st的两个平面（应该是当前面和包围盒上的某个面），与plane_cutting一起计算出三重交点p，根据p与s、t的距离判断来将p存入existing_vts或new_vts。
挺绕的，可以说是相当地严谨…不对不对，少了种相切情况的if判断，是后来师弟提出一个问题：为什么不考虑new_vts.size() == 2的情况？当时我少了根筋给了他否定的答复，但现在重读代码，感觉确实是存在候选切面与当前面两个边相交于新点的情况，不知为什么作者没加入这个判断。总之，根据compute_intersections()计算的交点个数，可以确定当前面与候选切面是否会相交，与当前面index_face_i的相交切面将存入到intersecting_faces中。

平面分割

开始对当前面进行切割。以faces_to_be_cut存储待分割的面，按序取一个候选切面进行切割后，原来的面不再存在，它被多个新的面片取代，将它们存入faces_to_be_cut，并在下一轮循环中使用下一个候选切面来切割。

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
//3、使用与 index_face_i 相交的平面分割 index_face_i 
//每次切割后，原来的面不再存在，它被多个块取代，每个块都将被另一个平面切割
std::vector<Polygon_mesh::Face_index> faces_to_be_cut;//存储需要被分割的面
faces_to_be_cut.push_back(index_face_i);//先添加当前的第i个面
while (!intersecting_faces.empty())
{
	//以第一个相交面作为切割面
	Polygon_mesh::Face_index cutting_face = *(intersecting_faces.begin());//切割面的索引
	const Kernel::Plane_3* cutting_plane = face_supporting_planes[cutting_face];//获取切割面所在平面

	std::set<Polygon_mesh::Face_index> new_faces;                //存储切割后产生的新面
	std::set<Polygon_mesh::Face_index> remained_faces;        //存储还需被切割的面
	for (std::size_t j = 0; j < faces_to_be_cut.size(); ++j)
	{
		Polygon_mesh::Face_index current_face = faces_to_be_cut[j];//被切割的面
		std::vector<Polygon_mesh::Face_index> tmp = split_plane(current_face, cutting_plane, mesh_candidate);//调用split_plane函数对这个面进行切割
		new_faces.insert(tmp.begin(), tmp.end());//把这些新的面加入到new_faces中
		if (tmp.empty())
		{
			remained_faces.insert(current_face);//没有发生切割，将该面加入到 remained_faces 待后续进行切割
		}
	}

	//把 new_faces 中的面作为下一轮需要被切割的面
	faces_to_be_cut = std::vector<Polygon_mesh::Face_index>(new_faces.begin(), new_faces.end());

	//把 remained_faces 中的面也加入到需要被切割的面
	faces_to_be_cut.insert(faces_to_be_cut.end(), remained_faces.begin(), remained_faces.end());

	intersecting_faces.erase(cutting_face);//将被切割面从相交的面中移除
}


//4、所有和第i个面相交的面 cutting_faces 都会被第i个面所在的平面 index_face_i 切割
for (std::size_t j = 0; j < cutting_faces.size(); ++j)
	split_plane(cutting_faces[j], plane_face_i, mesh_candidate);

核心的分割过程在HypothesisPlaneProcess::split_plane()函数和HypothesisPlaneProcess::split_edge()函数中，本质是使用CGAL::Euler::split_edge()来在要分割的边上插入一个新顶点，并返回指向新顶点的半边索引；用CGAL::Euler::split_face()在给定面上插入一条新边，最后为新的两个面片赋予所属平面、所属区域等属性。

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
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
std::vector<Polygon_mesh::Face_index> HypothesisPlaneProcess::split_plane(Polygon_mesh::Face_index face, const Kernel::Plane_3 * cutting_plane, Polygon_mesh & mesh)
{
	std::vector<Polygon_mesh::Face_index> new_faces;//容器_切割后产生的新面的索引

	//1、定义属性映射
	typename Polygon_mesh::template Property_map<Polygon_mesh::Face_index, const Kernel::Plane_3*> face_supporting_planes =
		mesh.template property_map<Polygon_mesh::Face_index, const Kernel::Plane_3*>("f:supp_plane").first;//属性映射_每个面的支撑面

	const Kernel::Plane_3* supporting_plane = face_supporting_planes[face];
	if (supporting_plane == cutting_plane)
		return new_faces;//给定面所在的平面和切割平面相同，则不需要切割，返回空向量


	typename Polygon_mesh::template Property_map<Polygon_mesh::Face_index, CGAL::internal::Planar_segment<Kernel>*> face_supporting_segments =
		mesh.template property_map<Polygon_mesh::Face_index, CGAL::internal::Planar_segment<Kernel>*>("f:supp_segment").first;//属性映射_每个面的切割面

	typename Polygon_mesh::template Property_map<Polygon_mesh::Edge_index, std::set<const Kernel::Plane_3*> > edge_supporting_planes =
		mesh.template property_map<Polygon_mesh::Edge_index, std::set<const Kernel::Plane_3*> >("e:supp_plane").first;//属性映射_每条边的支撑面
	CGAL::internal::Planar_segment<Kernel>* supporting_segment = face_supporting_segments[face];

	//2023.11.23new:添加新的属性映射以表示平面所属的区域
	Polygon_mesh::Property_map<Polygon_mesh::Face_index, std::size_t> face_region_map
		= mesh.property_map<Polygon_mesh::Face_index, std::size_t>("f:region_map").first;//每个面的所属区域
	size_t index_region = face_region_map[face];
	//—————————————————————————


	//2、计算交点，分情况分析切割情况
	std::vector<Polygon_mesh::Vertex_index> existing_vts;//容器_与切割面相交的已有顶点
	std::vector<EdgePos> new_vts;//容器_与切割面相交的新生成顶点
	compute_intersections(mesh, face, cutting_plane, existing_vts, new_vts);//计算给定面和切割平面的交点


	if (existing_vts.size() + new_vts.size() != 2)
		return new_faces;//交点数没有两个，则切割失败，返回空向量
	else if (existing_vts.size() == 2)
	{
		if (existing_vts[0] == existing_vts[1])
			return new_faces;//两交点为同一个顶点，返回空向量

		if (halfedge_exists(existing_vts[0], existing_vts[1], mesh))
			return new_faces;//两交点之间已存在一条边，返回空向量
	}

	//确定两交点的半边
	Polygon_mesh::Halfedge_index h0 = Polygon_mesh::null_halfedge();
	Polygon_mesh::Halfedge_index h1 = Polygon_mesh::null_halfedge();
	if (existing_vts.size() == 2)
	{
		//cutting_plane 在两个已有顶点处切割给定面（而不是一条边）
		h0 = mesh.halfedge(existing_vts[0]);//第一个已有顶点所在半边为h0
		h1 = mesh.halfedge(existing_vts[1]);//第二个已有顶点所在半边为h1
	}
	else if (existing_vts.size() == 1)
	{
		//cutting_plane 在一个已有顶点和一条边的内部切割给定面
		h0 = mesh.halfedge(existing_vts[0]);
		h1 = split_edge(mesh, new_vts[0], cutting_plane);//以新生成顶点的边为h1
	}
	else if (new_vts.size() == 2)
	{
		//cutting_plane 在两条边的内部切割给定面
		h0 = split_edge(mesh, new_vts[0], cutting_plane);
		h1 = split_edge(mesh, new_vts[1], cutting_plane);
	}
	CGAL_assertion(h0 != Polygon_mesh::null_halfedge());//断言_h0、h1不是空的半边索引
	CGAL_assertion(h1 != Polygon_mesh::null_halfedge());

	//为了分割平面，边h0和h1必须附着在同一个面上
	if (mesh.face(h0) != face)
	{
		//h0所在的面不是给定面
		Polygon_mesh::Halfedge_index end = h0;//记录h0的初始值，用于结束循环
		do {
			h0 = mesh.opposite(mesh.next(h0));//将h0沿着给定面的边界顺时针移动一步
			if (mesh.face(h0) == face)
				break;
		} while (h0 != end);
	}
	CGAL_assertion(mesh.face(h0) == face);//断言_h0所在的面是给定面

	if (mesh.face(h1) != face)
	{
		Polygon_mesh::Halfedge_index end = h1;
		do {
			h1 = mesh.opposite(mesh.next(h1));
			if (mesh.face(h1) == face)
				break;
		} while (h1 != end);
	}
	CGAL_assertion(mesh.face(h1) == face);

	//以h0和h1为参数，在给定面上插入一条新边，返回新边对应的一个半边索引
	Polygon_mesh::Halfedge_index h = CGAL::Euler::split_face(h0, h1, mesh);
	if (h == Polygon_mesh::null_halfedge() || mesh.face(h) == Polygon_mesh::null_face())
	{
		std::cerr << "Fatal error. could not split face" << std::endl;//错误_边插入失败/新边没有对应的有效面
		return new_faces;
	}


	//3、平面分割/赋值属性映射
	//在边h的属性映射_支撑面中加入给定面和切割平面
	Polygon_mesh::Edge_index e = mesh.edge(h);
	edge_supporting_planes[e].insert(supporting_plane);
	edge_supporting_planes[e].insert(cutting_plane);
	CGAL_assertion(edge_supporting_planes[e].size() == 2);//断言_边h所在的平面集合中只有两个元素

	// Now the two faces
	Polygon_mesh::Face_index f1 = mesh.face(h);//获取新边所在的一个半边对应的面索引
	face_supporting_segments[f1] = supporting_segment;//将给定面所在的平面段赋值给f1所在的平面段
	face_supporting_planes[f1] = supporting_plane;//将给定面所在的平面赋值给f1所在的平面
	face_region_map[f1] = index_region;//给定面对应的点云区域赋值给f1
	new_faces.push_back(f1);//将f1加入到容器_新面中

	Polygon_mesh::Face_index f2 = mesh.face(mesh.opposite(h));//获取新边所在的另一个半边对应的面索引
	face_supporting_segments[f2] = supporting_segment;
	face_supporting_planes[f2] = supporting_plane;
	new_faces.push_back(f2);
	face_region_map[f2] = index_region;

	return new_faces;
}
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
Polygon_mesh::Halfedge_index HypothesisPlaneProcess::split_edge(Polygon_mesh & mesh, const EdgePos & ep, const Kernel::Plane_3 * cutting_plane)
{
	// 这个函数将分割由 'ep' 表示的边
	// - 给新生成的边分配支撑平面
	// - 返回指向新生成顶点的半边索引
	// 函数内部使用了Euler split_edge()操作

	//1、获取mesh的属性映射
	typename Polygon_mesh::template Property_map<Polygon_mesh::Edge_index, std::set<const Kernel::Plane_3*> > edge_supporting_planes
		= mesh.template property_map<Polygon_mesh::Edge_index, std::set<const Kernel::Plane_3*> >("e:supp_plane").first;//属性映射_每个边所在的支撑面

	typename Polygon_mesh::template Property_map<Polygon_mesh::Vertex_index, std::set<const Kernel::Plane_3*> > vertex_supporting_planes
		= mesh.template property_map<Polygon_mesh::Vertex_index, std::set<const Kernel::Plane_3*> >("v:supp_plane").first;//属性映射_每个顶点所在的支撑面

	std::set<const Kernel::Plane_3*> sfs = edge_supporting_planes[ep.edge];//引用_分割边所在的平面//不使用常量引用，因为会在分割后失效
	CGAL_assertion(sfs.size() == 2);//断言_集合仅有两个元素（分割的两个面）


	//2、分割操作
	Polygon_mesh::Halfedge_index h = CGAL::Euler::split_edge(mesh.halfedge(ep.edge), mesh);//在要分割的边上插入一个新顶点，并返回指向新顶点的半边索引
	if (h == Polygon_mesh::null_halfedge())
		return h;//分割失败，返回空的半边索引

	Polygon_mesh::Vertex_index v = mesh.target(h);//新顶点对应的顶点索引
	if (v == Polygon_mesh::null_vertex())
		return Polygon_mesh::null_halfedge();//索引不存在，返回空的半边索引

	typename Polygon_mesh::template Property_map<Polygon_mesh::Vertex_index, Kernel::Point_3>& coords
		= mesh.points();//属性映射_每个顶点坐标
	coords[v] = *ep.pos;//设置新顶点的坐标为EdgePos对象中存储的交点坐标

	//获取新生成的两条边索引和其所在的平面
	Polygon_mesh::Edge_index e1 = mesh.edge(h);
	edge_supporting_planes[e1] = sfs;
	Polygon_mesh::Edge_index e2 = mesh.edge(mesh.next(h));
	edge_supporting_planes[e2] = sfs;

	//分割边所在的平面和切割面加入到新顶点的平面集合中
	vertex_supporting_planes[v] = sfs;
	vertex_supporting_planes[v].insert(cutting_plane);
	CGAL_assertion(vertex_supporting_planes[v].size() == 3);//断言_新顶点的平面集合中有三个元素

	return h;
}

经平面相交细分后，mesh_candidate中的大块面片被切分为若干细碎面片。它们仍称为候选面片集，其中存在某种能表示为建筑表面模型的最优组合，等待着我们去发掘。
候选平面相交细分

候选面集简化

先叠个甲：这一章节的内容并非原PolyFit的内容，是我自己加进去的“改进”内容，目的是去除明显冗余的候选面片集。
从上一小节的尾图可以明显看到，最外层包围盒附近的一圈面片是“被截断的”，这种面片存在某条边不与别的面片相连。由于最终构成的建筑模型是边缘闭合的水密模型，这些被截断的面片显然不属于最终的模型组合，所以多加了这一节的内容对面集进行简化，以减轻后续的计算压力。更学术的说法是对候选面集添加“XX约束”，减少暴力分割过程中产生的明显冗余面，能在[21]、[22]等文献中看到这类优化思路的改进PolyFit算法。
ok来看代码中5.1的HypothesisPlaneProcess::SimpleMeshCandidate_V2函数，整个过程分为区域归类、标记缓冲区面片和边缘面查询三步，算法思路简述如下：

1
2
3
4
①区域归类。对候选平面相交细分产生的每个分割面添加所属区域点云的属性映射，根据属性映射值查找属同一区域的多边形，将区域点云投影至多边形组成的平面。
②标记缓冲区面片。对平面上的每个多边形，判断其内部是否存在投影点并进行标记。提取标记多边形的顶点，对所有使用该类顶点的多边形进行标记。该步骤是为了扩展被标记的候选平面，从而保证最终模型的闭合性。
③边缘面查询。查询缓冲区面片集合中每个顶点被用于构成平面的次数。以次数为1的顶点所构成的平面视为边缘平面，需将其从面集内剔除。
④重复步骤③，直至所有顶点被用于构成平面的次数不为1。输出该多边形网络为简化后的候选平面集合

不过，上述思路是V1版本的。V1版本的代码写得过于激进，现在回看简直一坨。我记得是先统计mesh_candidate中所有顶点所邻接的候选面个数，若有的点只被唯一一个面“所拥有”，则认为该点的相关面为边缘面并予以剔除。方法粗看没什么问题，但在执行这个处理前，我对mesh_candidate做了次初筛，可能破坏了原mesh内部很多能构成封闭结构的面片，导致很容易把候选面删除到只剩一块块分离的小团块。V1先放那了仅做留档参考，V2在剔除面的判断上做了点改动，仍然是一坨，它没实现理想效果。下面来看看代码中的整个流程。

区域归类

先套一层遍历，每轮循环取一块聚类点云。每个聚类区域的点下标在之前被记录至分段分段s，以此提取出一个聚类的点坐标到points_region中。然后，使用supporting_plane()函数来拟合平面，该函数内部是调用CGAL::linear_least_squares_fitting_3的最小二乘拟合方法来计算平面系数。再然后，调用to_2d()方法，将聚类的每个点投影到平面上并存储到points_region_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
std::vector< Kernel::Point_3> vi_covered;//点云覆盖区域及邻近的点下标
for (size_t i_segment = 0; i_segment < _cloud_plane->planar_segments().size(); i_segment++)
{
	//1、获取区域点的覆盖区域
	CGAL::internal::Planar_segment<Kernel>* s = _cloud_plane->planar_segments()[i_segment];//指针_分段区域

	//获取区域点云
	std::vector<Kernel::Point_3> points_region;
	for (std::size_t i = 0; i < s->size(); ++i)
	{
		std::size_t idx = s->at(i);//点索引
		Kernel::Point_3 p = _cloud_plane->point_map()[idx];
		points_region.push_back(p);
	}

	//三维点投影到平面
	const Kernel::Plane_3* plane_fitting = _cloud_plane->planar_segments()[i_segment]->supporting_plane();//拟合平面
	std::vector<Kernel::Point_2> points_region_2;
	for (const Kernel::Point_3& p : points_region)
		points_region_2.push_back(plane_fitting->to_2d(p));


	//2、...
}

标记缓冲区面片

候选的大块面片在两两细分前，我用f:region_map记录了面片所属的聚类，根据下标把所有属于本轮聚类点云的候选面片都找出来，存入到faces_polygon中，相当于把处理对象从整个复杂的mesh_candicate收束到一个平面上。接着，对faces_polygon中的每个面f，将它上面的轮廓点集points_polygon同样投影到拟合面上，将得到的二维点存入points_polygon_2。根据这些点，判断聚类点points_region_2是否落在f内，只要有一个聚类点会投影到f上，那该f需被记录，将它的顶点存入到vi_covered中。

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
//2、遍历多边形上的每个面，将包含投影点及相关的面标记到 faces_covered
std::vector<Polygon_mesh::Face_index> faces_polygon;//存储多边形平面下标
for (const Polygon_mesh::Face_index& f : mesh_candidate.faces())
	if (mesh_candidate.property_map<Polygon_mesh::Face_index, std::size_t>("f:region_map").first[f] == i_segment)
		faces_polygon.push_back(f);


GeometryCal obj_gc;
for (const Polygon_mesh::Face_index& f : faces_polygon)
{
	std::vector<Kernel::Point_3> points_polygon;
	std::vector<Polygon_mesh::Vertex_index> vi_polygon;

	//类模板，创建顶点环绕面的两个索引器
	//vcirc (h, sm)表示从h开始环绕mesh中面的顶点；done(vcirc)用于判断循环是否结束
	CGAL::Vertex_around_face_circulator<Polygon_mesh> vcirc(mesh_candidate.halfedge(f), mesh_candidate), done(vcirc);
	do {
		Polygon_mesh::Vertex_index vi = *vcirc;//顶点索引
		Kernel::Point_3 p = mesh_candidate.point(vi);//顶点坐标
		vi_polygon.push_back(vi);
		points_polygon.push_back(p);

		++vcirc;
	} while (vcirc != done);

	//多边形顶点投影到平面
	std::vector<Kernel::Point_2> points_polygon_2;
	for (const Kernel::Point_3& p : points_polygon)
		points_polygon_2.push_back(plane_fitting->to_2d(p));

	//判断是否有区域点云落在多边形内
	for (Kernel::Point_2 p : points_region_2)
		if (obj_gc.isPointInPolygon(p, points_polygon_2))
		{
			//遍历该面 f 的顶点，存入
			for (const Polygon_mesh::Vertex_index& vi : vi_polygon)
			{
				Kernel::Point_3 p = mesh_candidate.point(vi);//顶点坐标
				if (std::find(vi_covered.begin(), vi_covered.end(), p) == vi_covered.end())
					vi_covered.push_back(p);
			}

			break;//该面 f 的顶点数据已存储，无需再遍历 points_region_2
		}
}

函数GeometryCal::isPointInPolygon()是GPT给的函数，思路参考射线法[23]：从当前点引一条射线穿过多边形边界，记录其穿过的次数，交叉数为偶数时点在外侧，奇数时点在内侧。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
bool GeometryCal::isPointInPolygon(const Kernel::Point_2 & point, const std::vector<Kernel::Point_2>& polygon)
{
	size_t n = polygon.size();//多边形的顶点个数
	int count = 0;//交点个数

	for (int i = 0; i < n; i++)
	{
		//获取边起点和终点
		Kernel::Point_2 p1 = polygon[i];
		Kernel::Point_2 p2 = polygon[(i + 1) % n];
		if (point.y() == p1.y() && point.y() == p2.y() && //点在水平边上
			point.x() >= std::min(p1.x(), p2.x()) && point.x() <= std::max(p1.x(), p2.x()))
			return true;//点在多边形内部

		if (point.y() > std::min(p1.y(), p2.y()) && point.y() <= std::max(p1.y(), p2.y()))
			if (isLineIntersect(point, Kernel::Point_2(point.x() + 1e9, point.y()), p1, p2))
				count++;//点的纵坐标在边的纵坐标范围内+水平线与边相交
	}

	return count % 2 == 1; // 如果交点个数是奇数，点在多边形内部，否则不在
}

此时vi_covered中存储有若干点，它们可构成聚类点云能投影到的候选面片。仍使用半边环绕器逐个取每个f上的点，当该点在vi_covered中有记录就把这个面存入缓冲面集faces_covered中。这部分写得很冗余，从V1那照抄套用属实难绷。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
//3、统计覆盖区域的顶点并保留关联的面，剔除边缘面
std::vector<Polygon_mesh::Face_index> faces_covered;//点云覆盖区域及邻近的面下标
for (const Polygon_mesh::Face_index& f : mesh_candidate.faces())
{
	if (std::find(faces_covered.begin(), faces_covered.end(), f) != faces_covered.end())
		continue;//该面 f 已存入 faces_covered，跳过

	CGAL::Vertex_around_face_circulator<Polygon_mesh> vcirc(mesh_candidate.halfedge(f), mesh_candidate), done(vcirc);
	do {
		Polygon_mesh::Vertex_index vi = *vcirc;//顶点索引
		Kernel::Point_3 p = mesh_candidate.point(vi);//顶点坐标

		if (std::find(vi_covered.begin(), vi_covered.end(), p) != vi_covered.end())
		{
			faces_covered.push_back(f);//faces_covered 中未存在面 f ，将其存入
			break;
		}

		++vcirc;
	} while (vcirc != done);
}

边缘面查询

有了缓冲面集faces_covered，那mesh_candidate中的其他候选面便视为边缘面，将它们下标记录为faces_delete后，使用CGAL::Euler::remove_face()将其移除，并在最后调用collect_garbage()对网络做一次修复。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
//4、简化面
std::vector<Polygon_mesh::Face_index> faces_delete;//面索引_记录需要删除的面
for (const Polygon_mesh::Face_index& f : mesh_candidate.faces())
{
	if (std::find(faces_covered.begin(), faces_covered.end(), f) != faces_covered.end())
		continue;//该面f属于faces_covered，不用删除

	faces_delete.push_back(f);
}

for (std::size_t i = 0; i < faces_delete.size(); ++i)
{
	Polygon_mesh::Face_index f = faces_delete[i];
	Polygon_mesh::Halfedge_index h = mesh_candidate.halfedge(f);
	CGAL::Euler::remove_face(h, mesh_candidate);//根据半边索引删除相应面
}

mesh_candidate.collect_garbage();//重新组织网络，移除被删除的面
候选面集简化

计算候选面的置信度

添加额外的平面属性放下一节详述。那么现在，mesh_candidate中的多边形网络结构如上图所示，原来的多个大块面片经聚类点云切出外轮廓，又经过两两细分变为若干细碎的候选面片。为了从这堆碎片中选出一个最能描述物体几何形状的子集，作者定义了一种称为能量函数的方法来优化这个选择过程，该函数可以由多个能量项构成。
能量函数构成
在作者的设计中，能量函数由数据拟合项、模型复杂性项和点覆盖项这三个能量项构成。数据拟合项评估由点簇所拟合平面的置信度，该值根据候选面内的平面点距及点与面的法向量差值来确定，可反映对应点云的平坦度及信噪比。模型复杂性项定义为模型中锐利边缘的比例。二元函数𝑐𝑜𝑟𝑛𝑒𝑟根据相连面属性决定边e是否为锐边，由于面片均以两两相交分割得到，非孤立边与面片的所有拓扑关系如下图所示（图源[2]），其中图b、c表示边能将相连面组合为一个更大的多边形，𝑐𝑜𝑟𝑛𝑒𝑟值为1；图d到g中的边可表示为模型边界，𝑐𝑜𝑟𝑛𝑒𝑟值为0。点覆盖项指模型中未覆盖点云的面片占比，它根据点簇在拟合面上的面积与拟合面总面积比值来确定。
候选面拓扑关系
能量函数的优化过程属于后一大章的内容先按下不表，先把注意力放在其中的能量项参数上。按公式要求，每个候选面需要增加平面的支撑点数f:num_supporting_points、平面面积f:face_area和覆盖面积f:covered_area这三个新的面片属性。在代码中我是直接调用作者写好的方法，这块的计算对那时的我来说太过耦合，主要是覆盖面积的计算还涉及到Alpha-shape边界提取，这属于另一块CGAL库的算法不太好搬，应该是想了想觉得不会改动里面的计算公式就没把它迁移出来封装。

1
2
3
4
5
6
7
//6、计算候选面的置信度
//该操作为 mesh_candidate 添加了3个属性映射,用于后续的MIP处理：
//平面的支撑点数	f:num_supporting_points
//平面面积			f:face_area
//覆盖面积			f:covered_area
CGAL::internal::Candidate_confidences<Kernel> conf;
conf.compute(cloud_plane, mesh_candidate);

计算置信度源码

候选面属性计算的主流程位于compute_confidences.h的Candidate_confidences::compute()，源码如下：

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
template <typename Kernel>
void Candidate_confidences<Kernel>::compute(const Point_set& point_set, Polygon_mesh& mesh) {
	const unsigned int K = 6;

	const typename Point_set::Point_map& points = point_set.point_map();
	FT avg_spacing = compute_average_spacing<Concurrency_tag>(points, K);

	// The number of supporting points of each face
	typename Polygon_mesh::template Property_map<Face_descriptor, std::size_t> face_num_supporting_points =
		mesh.template add_property_map<Face_descriptor, std::size_t>("f:num_supporting_points").first;

	// The area of each face
	typename Polygon_mesh::template Property_map<Face_descriptor, FT> face_areas =
		mesh.template add_property_map<Face_descriptor, FT>("f:face_area").first;

	// The point covered area of each face
	typename Polygon_mesh::template Property_map<Face_descriptor, FT> face_covered_areas =
		mesh.template add_property_map<Face_descriptor, FT>("f:covered_area").first;

	// The supporting plane of each face
	typename Polygon_mesh::template Property_map<Face_descriptor, const Plane*> face_supporting_planes =
		mesh.template property_map<Face_descriptor, const Plane*>("f:supp_plane").first;

	//距离阈值_判断平面是否太小或接近于0。snap_squared_distance_threshold函数会返回一个距离阈值，判断两个点是否足够接近
	FT degenerate_face_area_threshold = CGAL::snap_squared_distance_threshold<FT>() * CGAL::snap_squared_distance_threshold<FT>();

	for (auto f : mesh.faces()) {
		const Plane* supporting_plane = face_supporting_planes[f];//获取支撑平面=最适合拟合该面的平面
		// Face area
		FT area = face_area(f, mesh);//计算面积
		face_areas[f] = area;

		if (area > degenerate_face_area_threshold) {
			const std::vector<std::size_t>& indices = supporting_points(f, mesh, point_set);//获取每个面上支撑点的索引，即离该面最近点集的点
			face_num_supporting_points[f] = indices.size();

			std::vector<Point> pts;//容器_支撑点坐标
			for (std::size_t i = 0; i < indices.size(); ++i) {
				std::size_t idx = indices[i];
				const Point& p = points[idx];
				pts.push_back(p);
			}

			FT covered_area(0);//double值，存储每个面被支撑点覆盖的面积
			Alpha_shape_mesh<Kernel> alpha_mesh(pts.begin(), pts.end(), *supporting_plane);//构建α网络，包含所有支撑点且不包含空洞的网络
			Polygon_mesh covering_mesh;
			FT radius = avg_spacing * FT(5.0);
			if (alpha_mesh.extract_mesh(radius * radius, covering_mesh)) {
				// We cannot use the area of the 3D faces, because the alpha shape mesh is
				// not perfectly planar
				//根据半径值 radius 提取表面网络，获取该网络每个顶点的坐标
				const typename Polygon_mesh::template Property_map<Vertex_descriptor, Point>& coords = covering_mesh.points();
				for (auto face : covering_mesh.faces()) {
					// We have to use the projected version
					Polygon plg; // the projection of the face onto it supporting plane
					Halfedge_around_face_circulator<Polygon_mesh> cir(covering_mesh.halfedge(face), covering_mesh), done(cir);
					do {
						//通过调取多边形半边，获取顶点的二维坐标，将之加入到多边形中
						Halfedge_descriptor hd = *cir;
						Vertex_descriptor vd = covering_mesh.target(hd);
						const Point& p = coords[vd];
						const Point2& q = supporting_plane->to_2d(p);
						plg.push_back(q);
						++cir;
					} while (cir != done);
					covered_area += std::abs(plg.area());//计算多边形面积，累加到覆盖面积中，使用绝对值避免负数情况
				}
			}

			face_covered_areas[f] = covered_area;//每个面的覆盖面积存储到属性映射face_covered_areas
			if (covered_area > area)
				face_covered_areas[f] = area;//覆盖面积大于原始面积，存在某些超出范围的支撑点，将面积设为原始面积
		}
		else { // For tiny faces, we can simple assign zero supporting points
			face_num_supporting_points[f] = 0;//平面面积小于等于阈值，将每个面上支撑点的数量设置为零
			face_covered_areas[f] = FT(0.0);//将每个面的覆盖面积设置为零
		}
	}
}

添加额外的平面属性

这一小节是自己添加的内容，给候选面多增加一个离顶面距离的属性，试图为能量函数多加一个新能量项做准备。思路很简单，取每个面片的中心点，计算点到包围盒顶面的距离，再将该距离除以包围盒上下面的距离来归一化。实现函数在HypothesisPlaneProcess::AddMeshProperty()。

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
void HypothesisPlaneProcess::AddMeshProperty(Polygon_mesh & mesh_candidate)
{
	typename Polygon_mesh::template Property_map<Polygon_mesh::Face_index, double> face_roof =
		mesh_candidate.template add_property_map<Polygon_mesh::Face_index, double>("f:dist_boxTop").first;//属性映射_面片属顶面

	//候选面顶点集
	std::vector<Kernel::Point_3> vertices(mesh_candidate.number_of_vertices());
	for (int i = 0; i < mesh_candidate.vertices().size(); ++i)
	{
		const Polygon_mesh::Vertex_index v = mesh_candidate.vertices().first[i];
		vertices[i] = mesh_candidate.points()[v];
	}

	//计算包围盒
	const Kernel::Iso_cuboid_3& box = CGAL::bounding_box(vertices.begin(), vertices.end());//计算 vertices 点集的最小包围盒
	FT dx = box.xmax() - box.xmin();//包围盒xyz方向上的长度
	FT dy = box.ymax() - box.ymin();
	FT dz = box.zmax() - box.zmin();

	//计算包围盒顶面
	std::vector<Kernel::Point_3> point_boxTop;//包围盒顶面点
	point_boxTop.push_back(Kernel::Point_3(box.xmin(), box.ymin(), box.zmax()));
	point_boxTop.push_back(Kernel::Point_3(box.xmax(), box.ymin(), box.zmax()));
	point_boxTop.push_back(Kernel::Point_3(box.xmax(), box.ymax(), box.zmax()));
	point_boxTop.push_back(Kernel::Point_3(box.xmin(), box.ymax(), box.zmax()));
	Kernel::Plane_3 plane_boxTop;
	CGAL::linear_least_squares_fitting_3(point_boxTop.begin(), point_boxTop.end(), plane_boxTop, CGAL::Dimension_tag<0>());

	//遍历候选面集的每个面片
	for (const Polygon_mesh::Face_index& f : mesh_candidate.faces())
	{
		double x_center = 0.0f, y_center = 0.0f, z_center = 0.0f;
		int num_p = 0;

		//累计该面上的点坐标
		CGAL::Vertex_around_face_circulator<Polygon_mesh> vcirc(mesh_candidate.halfedge(f), mesh_candidate), done(vcirc);
		do {
			Polygon_mesh::Vertex_index vi = *vcirc;//顶点索引
			Kernel::Point_3 p = mesh_candidate.point(vi);//顶点坐标

			x_center += p.x();
			y_center += p.y();
			z_center += p.z();
			++num_p;

			++vcirc;
		} while (vcirc != done);

		//面片中心点
		Kernel::Point_3 p_center(x_center / num_p, y_center / num_p, z_center / num_p);

		//点到平面的距离 / 包围盒高度
		//除以包围盒高度是想归一化处理，让值处于[0,1]
		double dist_p2plane = std::sqrt(CGAL::squared_distance(p_center, plane_boxTop)) /(box.zmax() - box.zmin());

		//属性映射赋值
		face_roof[f] = dist_p2plane;
	}

}

该属性目的是想提高“屋顶”候选面片在能量函数中的权重，当该项属性越接近0则表明面片越靠近顶部，若最终模型选取这类面的数量越多，那可以期望其顶面细节越丰富。不过单靠这个值是不可靠的，应该还需要加某种约束条件。从以前的测试结果来看并不理性，仍把代码放进来，仅为提供一个添加平面属性的参考。
至此，上一大节的区域朝向面片已被处理为候选面集。下一大节解析能量函数，看他如何从杂乱无章的候选面集挑选出最终模型。
区域朝向面片2候选面集

参考文献&引用

[1] https://3d.bk.tudelft.nl/liangliang/
[2] https://3d.bk.tudelft.nl/liangliang/publications/2017/polyfit/polyfit.html
[3] https://www.youtube.com/watch?v=_0brfDFkIkc
[4] Linfu X ,Han H ,Qing Z , et al. Combined Rule-Based and Hypothesis-Based Method for Building Model Reconstruction from Photogrammetric Point Clouds [J]. Remote Sensing, 2021, 13 (6): 1107-1107.
[5] Nan L , Wonka P .PolyFit: Polygonal Surface Reconstruction from Point Clouds[C]//International Conference on Computer Vision.IEEE, 2017.
[6] https://github.com/Kitware/CMake/releases/download/v3.29.3/cmake-3.29.3-windows-x86_64.zip
[7] https://boostorg.jfrog.io/artifactory/main/release/1.85.0/source/boost_1_85_0.7z
[8] https://github.com/CGAL/cgal/releases/download/v5.6.1/CGAL-5.6.1.zip
[9] https://gitlab.com/libeigen/eigen/-/archive/3.4.0/eigen-3.4.0.zip
[10] https://scipopt.org/download.php?fname=SCIPOptSuite-9.0.0-win64-VS15.exe
[11] https://zenodo.org/records/4390295#.Y0eIodJBxuV
[12] https://www.zhihu.com/question/277599635/answer/2149719454
[13] https://zhuanlan.zhihu.com/p/668272208
[14] https://doc.cgal.org/latest/Polygonal_surface_reconstruction/index.html#Chapter_PolygonalSurfaceReconstruction
[15] https://zhuanlan.zhihu.com/p/90858099
[16] https://zhuanlan.zhihu.com/p/668272208
[17] https://doc.cgal.org/latest/HalfedgeDS/index.html
[18] https://segmentfault.com/q/1010000043819426
[19] https://doc.cgal.org/latest/BGL/group__PkgBGLEulerOperations.html#gaa386d0cdef3b5d6ef43d6b503392dbcd
[20] https://doc.cgal.org/latest/Shape_regularization/index.html
[21] https://doi.org/10.3390/rs13061107
[22] https://doi.org/10.1016/j.isprsjprs.2022.09.017
[23] https://blog.csdn.net/u013279723/article/details/106265948