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

摘要：
聚类点云的区域优化；构造包围盒；生成候选平面集；区域朝向面片。

上一节对点云做了聚类处理，在这一节中要开始根据各个聚类区域来生成相应的平面多边形。
提前再叠个甲，后续部分的模型重建，代码基本均来源于

<CGAL/Polygonal_surface_reconstruction.h>
<CGAL/Polygonal_surface_reconstruction/internal/hypothesis.h>

这两个引用文件，它的作者也正是NanLiangLiang老师，我只按我的喜好把它们摘出来做了下封装，打了些注释说明，且印象中只在这个框架下的候选平面那块做了些许改动，所以不要指望它相较源程序有多大优化。

构建假设(候选)平面的整体流程

好，先回main()函数中看一眼，这里创建了一个CGAL::Surface_mesh<Kernel：：Point_3>类型的对象mesh_candidate，它是一个用于表示和处理三维表面网格的模板类。通过调用HypothesisPlaneProcess中的BuildHypothesisPlane()，可构建点云的候选(/假设)面片集合。

HypothesisPlaneProcess obj_hp;
CGAL::Surface_mesh<Kernel::Point_3> mesh_candidate;//多边形网络的候选表面集
obj_hp.BuildHypothesisPlane(cloud_seg, mesh_candidate);//构建点云的面片集合
//obj_hp.Coordination2Original(mesh_candidate, obj_pre._xOffset, obj_pre._yOffset, obj_pre._zOffset);//网络平移至原始坐标
obj_file.SaveColorOBJ(path_output, obj_file.GetFileName(file_path) + "_mesh_candidate", mesh_candidate);
//obj_file.SaveOBJ(path_output, obj_file.GetFileName(file_path) + "_mesh_candidate_nocolor", mesh_candidate);

围绕mesh_candidate，函数BuildHypothesisPlane()基于聚类点云cloud_seg进行面片的生成、处理，整个流程分为聚类点云的区域优化、构建包围盒面片、构建包围盒范围的候选平面、平面相交细分、候选面简化及面片属性计算这6个阶段。

void HypothesisPlaneProcess::BuildHypothesisPlane(const PointCloud & cloud_region, CGAL::Surface_mesh<Kernel::Point_3>& mesh_candidate)
{
	CGAL::Timer timer;//时间计数的实例
	timer.start();

	//0、转格式，cloud_region -> cloud_plane
	std::vector<PNI> point_pni;
	for (int i = 0; i < cloud_region.size(); ++i)
	{
		PNI temp;
		temp.get<0>() = cloud_region.point(i);
		temp.get<1>() = cloud_region.normal(i);
		temp.get<2>() = cloud_region.property_map<int>("region_map").first[i];
		point_pni.push_back(temp);
	}

	//2023.8.21 注意，cloud_plane要全程引用（等同于全局变量），不然容易触发~Point_set_with_planes()从而清掉 planar_segments
	CGAL::internal::Point_set_with_planes<Kernel> cloud_plane(point_pni,//点云对象
		CGAL::Nth_of_tuple_property_map<0, PNI>(),//指定类型为Point_map
		CGAL::Nth_of_tuple_property_map<1, PNI>(),//Normal_map
		CGAL::Nth_of_tuple_property_map<2, PNI>());//Plane_index_map

	if (cloud_plane.planar_segments().size() < 4)
	{
		std::cout << "错误：重建闭合表面模型至少需要4个平面聚类，该点云只有 " + std::to_string(cloud_plane.planar_segments().size()) + " 个聚类" << std::endl;
		return;//检查点集的平面段数量，如果小于4，则返回错误消息
	}
	_cloud_plane = &cloud_plane;


	//1、优化点云的平面聚类（通过合并压低数量）
	//std::cout <<"当前点云的平面段数量："<< _cloud_plane->planar_segments().size() << std::endl;
	PlaneRefine();
	//std::cout << "优化后的点云的平面段数量：" << _cloud_plane->planar_segments().size() << std::endl;


	//2、构建点云的包围盒多边形
	CGAL::Surface_mesh<Kernel::Point_3> mesh_bbox;
	BuildMeshBbox(mesh_bbox);


	//3、构建点云的候选表面多边形（跨幅延伸至包围盒的平面）
	BuildMeshCandidate(mesh_bbox, mesh_candidate);


	//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);

	///输出测试_显示新加的三种属性映射
	//std::cout << mesh_candidate.faces().size() << std::endl;
	//for (auto f : mesh_candidate.faces())
	//{
	//	std::cout << mesh_candidate.property_map<Polygon_mesh::Face_index, std::size_t>("f:num_supporting_points").first[f] << std::endl;
	//	std::cout << mesh_candidate.property_map<Polygon_mesh::Face_index, FT>("f:face_area").first[f] << std::endl;
	//	std::cout << mesh_candidate.property_map<Polygon_mesh::Face_index, FT>("f:covered_area").first[f] << std::endl;
	//	std::cout << std::endl;
	//}


	timer.stop();
	std::cout << "有 " << mesh_candidate.faces().size() << " 个候选表面已生成，运行耗时 " << timer.time() << " 秒" << std::endl;


	///输出测试_候选面的颜色属性//弃用
	//for (auto f : mesh_candidate.faces())
	//	std::cout << mesh_candidate.property_map<Polygon_mesh::Face_index, CGAL::IO::Color>("f:color").first[f] << std::endl;
}

这里演示下基于聚类点云生成的平面基元结果。

聚类点云的区域优化

这一个函数是想通过合并相近、法向近似的聚类来减少区域数量，从而降低计算负担。如何判断两个聚类需要合并？函数内设置了两个参数来简化聚类区域数量：点到区域最大距离的平均值avg_max_dist和区域间的法向量夹角_theta。整个流程主要分以下3步：
（1）计算每个聚类区域拟合平面的平均最大距离。先用最小二乘拟合聚类区域的平面，再记录点到平面的最大距离值max_dist，最后累计所有的max_dist到avg_max_dist，对这个值求平均再除以2后，得到一个衡量点是否应被纳入另一个平面范围内的值。

（2）获取拟合平面的法向量，判断法向量间的夹角是否小于_theta，满足判定的两个聚类朝向近似，可进一步做合并判断。_theta的设置无严格的理论依据，我按经验调整，夹角阈值调的越大则聚类越容易被合并，导致候选面集的面数越少，也就是细节越少。
代码中若对向量内积绝对值的含义仍有疑惑，可以暂停看一下这里的例子。

（3）统计平面s1上是否有足够多的点离s2比较近。平面上的一个点距离另一个平面是否比较近，通过函数calNumber_points_on_plane()进行判断。是否“足够多”则通过num_threshold来判断，该值为s1平面上的点数除以5得到，似乎也是个经验调参值，我未做改动。代码中s1->s2和s2->s1都做了一遍统计，只要有任一个的统计数量能满足num_threshold，则合并两个聚类。

主流程

void HypothesisPlaneProcess::PlaneRefine()
{
	std::vector< CGAL::internal::Planar_segment<Kernel>* >& segments = _cloud_plane->planar_segments();//获取 point_set_ 对象的 planar_segments() 成员变量的引用

	//1、遍历所有平面段，计算每个平面段的平均最大距离
	FT avg_max_dist = 0;
	for (std::size_t i = 0; i < segments.size(); ++i)
	{
		CGAL::internal::Planar_segment<Kernel>* s = segments[i];
		const Kernel::Plane_3* plane = s->fit_supporting_plane();//拟合平面(用户可能提供无效的平面拟合

		FT max_dist = -(std::numeric_limits<FT>::max)();//初始化最大距离为无穷小
		for (std::size_t j = 0; j < s->size(); ++j)
		{
			std::size_t idx = s->at(j);//点索引
			const Kernel::Point_3& p = _cloud_plane->point_map()[idx];
			FT sdist = CGAL::squared_distance(*plane, p);//计算点到平面的距离
			max_dist = (std::max)(max_dist, std::sqrt(sdist));//更新最大距离
		}

		avg_max_dist += max_dist;//累加最大距离
	}
	//计算平均最大距离
	avg_max_dist /= segments.size();
	avg_max_dist /= FT(2.0);


	//2、平面段合并
	//FT theta = static_cast<FT>(CGAL_PI * 10.0 / FT(180.0));//阈值 theta，弧度制//放公有变量处
	bool merged = false;//是否有平面段被合并
	do {
		merged = false;

		//将平面段按照大小进行排序
		//点较少的分段可信度较低，因此应首先合并
		std::sort(segments.begin(), segments.end(), CGAL::internal::SegmentSizeIncreasing<CGAL::internal::Planar_segment<Kernel>>());

		//遍历平面段，plane1为当前，plane2为下一个
		for (std::size_t i = 0; i < segments.size(); ++i)
		{
			CGAL::internal::Planar_segment<Kernel>* s1 = segments[i];
			const Kernel::Plane_3* plane1 = s1->supporting_plane();
			Kernel::Vector_3 n1 = plane1->orthogonal_vector();
			CGAL::internal::normalize<FT, Kernel::Vector_3>(n1);

			FT num_threshold = s1->size() / FT(5.0);//阈值 num_threshold
			for (std::size_t j = i + 1; j < segments.size(); ++j)
			{
				CGAL::internal::Planar_segment<Kernel>* s2 = segments[j];
				const Kernel::Plane_3* plane2 = s2->supporting_plane();
				Kernel::Vector_3 n2 = plane2->orthogonal_vector();
				CGAL::internal::normalize<FT, Kernel::Vector_3>(n2);

				//如果两个平面段的正交向量之间的内积大于 theta，
				if (std::abs(n1 * n2) > std::cos(_theta))
				{
					//且在另一平面上的点数大于 num_threshold，合并两平面
					std::size_t set1on2 = calNumber_points_on_plane(s1, plane2, avg_max_dist);
					std::size_t set2on1 = calNumber_points_on_plane(s2, plane1, avg_max_dist);
					if (set1on2 > num_threshold || set2on1 > num_threshold)
					{
						plane_merge(s1, s2);
						merged = true;
						break;
					}
				}
			}
			if (merged)
				break;
		}
	} while (merged);//merged 若为false，则此时所有平面都不适合被合并，跳出循环

	std::sort(segments.begin(), segments.end(), CGAL::internal::SegmentSizeDecreasing<CGAL::internal::Planar_segment<Kernel>>());//平面段按大小逆序排序

	//3、存储所有的支撑平面->_supporting_planes
	for (std::size_t i = 0; i < segments.size(); ++i)
	{
		CGAL::internal::Planar_segment<Kernel>* s = segments[i];
		const Kernel::Plane_3* plane = s->supporting_plane();
		_supporting_planes.push_back(plane);
	}
}

实现细节_统计共面的点数

函数calNumber_points_on_plane()，计算平面s1上所有点到s2的距离，统计其中值小于avg_max_dist的点数。这些点视为平面s1离平面s2的近点，当这些点数超过一个阈值后，视两个聚类为同一个聚类。

std::size_t HypothesisPlaneProcess::calNumber_points_on_plane(const CGAL::internal::Planar_segment<Kernel>* s, const Kernel::Plane_3 * plane, FT dist_threshold)
{
	CGAL_assertion(const_cast<CGAL::internal::Planar_segment<Kernel>*>(s)->point_set() == _cloud_plane);//检查参数有效性

	//1、初始化变量
	//获取 point_set_ 对象的 point_map() 成员变量的引用
	//可访问并修改 point_set_ 对象的所有点
	const CGAL::internal::Point_set_with_planes<Kernel>::Point_map& points = _cloud_plane->point_map();

	//2、遍历平面段 s 上的所有点
	std::size_t count = 0;
	for (std::size_t i = 0; i < s->size(); ++i)
	{
		std::size_t idx = s->at(i);
		const Kernel::Point_3& p = points[idx];

		//3、计算点 p 到平面 plane_cutting 的距离
		FT sdist = CGAL::squared_distance(*plane, p);
		FT dist = std::sqrt(sdist);
		if (dist < dist_threshold)//如果距离小于阈值，计数+1
			++count;
	}
	return count;
}

实现细节_聚类合并

聚类合并函数plane_merge()，源代码写得很规范，提取、合并两个聚类的点索引，并对新的聚类做一次最小二乘拟合，得到其平面系数。

void HypothesisPlaneProcess::plane_merge(CGAL::internal::Planar_segment<Kernel>* s1, CGAL::internal::Planar_segment<Kernel>* s2)
{
	//检查参数有效性
	CGAL_assertion(const_cast<CGAL::internal::Planar_segment<Kernel>*>(s1)->point_set() == _cloud_plane);
	CGAL_assertion(const_cast<CGAL::internal::Planar_segment<Kernel>*>(s2)->point_set() == _cloud_plane);
	//获取 point_set_ 对象的 planar_segments() 成员变量的引用
	//可访问并修改 point_set_ 对象的所有线段
	std::vector< CGAL::internal::Planar_segment<Kernel>* >& segments = _cloud_plane->planar_segments();

	//1、获取平面 s1 和 s2 的所有点的索引
	std::vector<std::size_t> points_indices;
	points_indices.insert(points_indices.end(), s1->begin(), s1->end());
	points_indices.insert(points_indices.end(), s2->begin(), s2->end());

	//2、创建一个新的平面 s，并设置它的所有点为 s1 + s2
	CGAL::internal::Planar_segment<Kernel>* s = new CGAL::internal::Planar_segment<Kernel>(_cloud_plane);
	s->insert(s->end(), points_indices.begin(), points_indices.end());
	s->fit_supporting_plane();//计算线段的支撑平面
	segments.push_back(s);//将新平面 s 添加到平面段 segments

	//3、从点集内删除平面 s1 和 s2
	typename std::vector< CGAL::internal::Planar_segment<Kernel>* >::iterator pos = std::find(segments.begin(), segments.end(), s1);
	if (pos != segments.end())
	{
		CGAL::internal::Planar_segment<Kernel>* tmp = *pos;
		const Kernel::Plane_3* plane = tmp->supporting_plane();
		segments.erase(pos);
		delete tmp;
		delete plane;
	}
	else
		std::cerr << "Fatal error: should not reach here" << std::endl;

	pos = std::find(segments.begin(), segments.end(), s2);
	if (pos != segments.end())
	{
		CGAL::internal::Planar_segment<Kernel>* tmp = *pos;
		const Kernel::Plane_3* plane = tmp->supporting_plane();
		segments.erase(pos);
		delete tmp;
		delete plane;
	}
	else
		std::cerr << "Fatal error: should not reach here" << std::endl;

}

思路复述

再复述一次合并的具体实现：先是对点云中的各个区域按点数量由小到大排序，再套两层遍历，从中取两个区域s1和s2，分别拟合平面并计算平面归一化后的法向量，再然后，按照向量间的夹角公式，当这两个法向量内积的绝对值大于角度阈值的余弦值时，表示两个区域近似平行。这里的std::abs表示取绝对值，在第一、二象限中余弦值越大表示向量间的夹角越小。
角度值满足条件后，再作判断，统计两个区域中小于avg的点数。若这两个的任一点数大于区域总点数的五分之一时，则对两个区域进行合并，实现上是创建一个新的区域，对原先的两个区域转存点云后擦除。这里的五分之一也是个经验参数。
当发生合并后，将合并的flag值修改为true并跳出上述的两层遍历。最后再套了层do while循环，以这个flag值为判断，表示会一直执行直至所有区域都无法执行合并后停止。

这一节对聚类点云做了个合并优化，点云点数是没变的，但点簇数/区域数量减少了，相应的拟合平面数量也就减少了，能稍稍缓解后续的运算压力。
此外，有个小细节记录下，合并后的聚类可以在空间分布上不连续，毕竟以它们合并前的聚类进行平面拟合，将得到的平面系数表征为平面也是无边界、朝向近似的，不需要额外做连通性判断。平面的边界划分在下一节进行。

构造包围盒

经上一节处理，理论上聚类点云的“分类数量”能有一定程度的减少。当然，到目前为止我们的处理对象仍为点云。到了本小节，开始以此为基础构建面元结构。

如前文所述，聚类点云按最小二乘/PCA等方法按平面方程做拟合，得到的是平面系数，它可视化呈现出来的是一大整块无边界的平面。为了避免后续产生冗余数据，需构建平面的多边形边界。
因此，作者给出了一种边界划分方案：以点云的包围盒边界切分出各个平面的边界。感觉挺合理的，这样就不会在点云外面的很远处生成无用面片。
构建包围盒的算法封装在HypothesisPlaneProcess::BuildMeshBbox()中，首先是构建一个（略大于）点云包围盒范围的多边形mesh_box。先计算出原始点云的轴对称包围盒，取体对角线的一定长度作为偏移量，将包围盒适当地向外扩展。然后以包围盒XYZ的最大最小值创建8个点，添加到mesh_box对象中，再以一定顺序连接这些点，构造出实质为12块面片、表现为6个面下标的立方体。

//1、计算包围盒半径、偏移量
const Kernel::Iso_cuboid_3& box = CGAL::bounding_box(_cloud_plane->point_map().begin(), _cloud_plane->point_map().end());//计算点云的包围盒

FT dx = box.xmax() - box.xmin();
FT dy = box.ymax() - box.ymin();
FT dz = box.zmax() - box.zmin();
FT radius = FT(0.5) * std::sqrt(dx * dx + dy * dy + dz * dz);
FT offset = radius * FT(0.05);


//2、创建一个更大的包围盒，以确保所有点都包含其中
FT xmin = box.xmin() - offset, xmax = box.xmax() + offset;
FT ymin = box.ymin() - offset, ymax = box.ymax() + offset;
FT zmin = box.zmin() - offset, zmax = box.zmax() + offset;


//3、更新 mesh 中的多边形网络为点集 points 的包围盒
mesh.clear();//清空多边形网络
//添加包围盒的8个顶点
Polygon_mesh::Vertex_index v0 = mesh.add_vertex(Kernel::Point_3(xmin, ymin, zmin));  // 0
Polygon_mesh::Vertex_index v1 = mesh.add_vertex(Kernel::Point_3(xmax, ymin, zmin));  // 1
Polygon_mesh::Vertex_index v2 = mesh.add_vertex(Kernel::Point_3(xmax, ymin, zmax));  // 2
Polygon_mesh::Vertex_index v3 = mesh.add_vertex(Kernel::Point_3(xmin, ymin, zmax));  // 3
Polygon_mesh::Vertex_index v4 = mesh.add_vertex(Kernel::Point_3(xmax, ymax, zmax));  // 4
Polygon_mesh::Vertex_index v5 = mesh.add_vertex(Kernel::Point_3(xmax, ymax, zmin));  // 5
Polygon_mesh::Vertex_index v6 = mesh.add_vertex(Kernel::Point_3(xmin, ymax, zmin));  // 6
Polygon_mesh::Vertex_index v7 = mesh.add_vertex(Kernel::Point_3(xmin, ymax, zmax));  // 7

//添加包围盒的6个面
mesh.add_face(v0, v1, v2, v3);
mesh.add_face(v1, v5, v4, v2);
mesh.add_face(v1, v0, v6, v5);
mesh.add_face(v4, v5, v6, v7);
mesh.add_face(v0, v3, v7, v6);
mesh.add_face(v2, v4, v7, v3);

有个要注意的点是，Surface_mesh<Kernel：：Point_3>的add_face()函数，点的添加顺序不宜乱动，因为这里面每个面的顶点都是按顺时针/逆时针排列的。乱改可能会因为合法性的问题无法正常生成面片，大概是因为里面的半边结构HalfEdge对面的合理性做了约束。
有关半边结构的详细理论可看下官网文档^[17]，我当时的理解是，多边形网络是由很多条边构成的嘛，但这些边可能也就只存储了两个点的索引，为了能增加边与边、边与面之间的拓扑关系，CGAL设计了半边结构这种类型，它在边的基础上设置了两条有向线段，两条有向线段的指向方向互相相反，分别以边的两个点为起点/终点，它们对应两个半边结构，称为对偶半边。
而在CGAL的多边形结构中，每个面由一个环形同向的半边链所定义，这些半边按顺序连接起来围成一个闭合环。当在已有闭合环的某一侧另外圈起一个新的“环”时（也就是在已有面的外缘构造新面），若点的顺序不对，就可能使新半边与其对偶半边的方向相同，这就会出现合法性问题，导致这个面无法被创建。

后面还有一部分代码，是对mesh_box这个对象附加属性，分别构建以面、边、点的索引为key值，平面对象为键值的属性，相当于将面、边、点与另外构建的平面做唯一绑定，这个另外构建的平面我称其为支撑平面，由于其相较mesh中的面片结构是没有边界的，能为后续的一些几何运算提供检索便利。

//添加多边形属性
typename Polygon_mesh::template Property_map<Polygon_mesh::Face_index, const Kernel::Plane_3*> face_supporting_planes =
	mesh.template add_property_map<Polygon_mesh::Face_index, const Kernel::Plane_3*>("f:supp_plane").first;//每个面的支撑平面

typename Polygon_mesh::template Property_map<Polygon_mesh::Edge_index, std::set<const Kernel::Plane_3*> > edge_supporting_planes
	= mesh.template add_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 add_property_map<Polygon_mesh::Vertex_index, std::set<const Kernel::Plane_3*> >("v:supp_plane").first;//每个顶点的支撑平面

//为网络中的每个面分配原始平面
for (auto fd : mesh.faces())
{
	Polygon_mesh::Halfedge_index h = mesh.halfedge(fd);
	Polygon_mesh::Vertex_index va = mesh.target(h);        const Kernel::Point_3& pa = mesh.points()[va]; h = mesh.next(h);
	Polygon_mesh::Vertex_index vb = mesh.target(h);        const Kernel::Point_3& pb = mesh.points()[vb]; h = mesh.next(h);
	Polygon_mesh::Vertex_index vc = mesh.target(h);        const Kernel::Point_3& pc = mesh.points()[vc];
	const Kernel::Plane_3* plane = new Kernel::Plane_3(pa, pb, pc);//使用三个顶点的坐标构造一个平面
	_supporting_planes.push_back(plane);
	face_supporting_planes[fd] = plane;
}

//为网络中的每个边分配原始平面
for (auto ed : mesh.edges())
{
	//获取每个边的两个半边
	Polygon_mesh::Halfedge_index h1 = mesh.halfedge(ed);
	Polygon_mesh::Halfedge_index h2 = mesh.opposite(h1);

	//获取半边各自所属的两个平面
	Polygon_mesh::Face_index f1 = mesh.face(h1);
	Polygon_mesh::Face_index f2 = mesh.face(h2);
	CGAL_assertion(f1 != Polygon_mesh::null_face());//bbox网格已关闭
	CGAL_assertion(f2 != Polygon_mesh::null_face());//bbox网格已关闭

	//使用这两个原始平面构造新的平面，存入 edge_supporting_planes 
	const Kernel::Plane_3* plane1 = face_supporting_planes[f1];
	const Kernel::Plane_3* plane2 = face_supporting_planes[f2];
	CGAL_assertion(plane1 && plane2 && plane1 != plane2);

	edge_supporting_planes[ed].insert(plane1);
	edge_supporting_planes[ed].insert(plane2);
	CGAL_assertion(edge_supporting_planes[ed].size() == 2);
}

//为网格中的每个顶点分配原始平面
for (auto vd : mesh.vertices())
{
	CGAL_assertion(vertex_supporting_planes[vd].size() == 0);
	CGAL::Halfedge_around_target_circulator<Polygon_mesh> hbegin(vd, mesh), done(hbegin);
	do {
		Polygon_mesh::Halfedge_index h = *hbegin;//获取每个顶点的所有半边
		Polygon_mesh::Face_index f = mesh.face(h);//获取该半边所属的面

		//基于该原始平面构造新的平面，存入 vertex_supporting_planes
		const Kernel::Plane_3* plane = face_supporting_planes[f];
		vertex_supporting_planes[vd].insert(plane);
		++hbegin;
	} while (hbegin != done);
	CGAL_assertion(vertex_supporting_planes[vd].size() == 3);
}

std::sort(_supporting_planes.begin(), _supporting_planes.end());

CGAL_assertion(mesh.is_valid());//检查多边形网络，确保其是有效的

构建候选表面多边形

有了包围盒多边形作为切割边界，开始基于聚类点云拟合出各个区域的平面。

属性声明

多边形构建的实现函数HypothesisPlaneProcess::BuildMeshCandidate()，迎面撞过来的是一大堆typename别名，他们都是多边形上各个对象与其对应的支撑平面或所属区域的属性。这里可以很明显的看到，候选平面的属性别名用的是add_property_map()方法，而包围盒的属性别名用的都是property_map()方法，这是因为在之前的构造中已创建过包围盒相应的属性，当然改用add_property_map()应该也没问题。

//1、定义 mesh_bbox 、 mesh_candidate 的属性
 Polygon_mesh::Property_map<Polygon_mesh::Edge_index, std::set<const Kernel::Plane_3*> > bbox_edge_supporting_planes
	= bbox_mesh.property_map<Polygon_mesh::Edge_index, std::set<const Kernel::Plane_3*> >("e:supp_plane").first;//每个边的所有支持平面

Polygon_mesh::Property_map<Polygon_mesh::Vertex_index, std::set<const Kernel::Plane_3*> > bbox_vertex_supporting_planes
	= bbox_mesh.property_map<Polygon_mesh::Vertex_index, std::set<const Kernel::Plane_3*> >("v:supp_plane").first;//每个顶点的所有支持平面


mesh_candidate.clear();//清空 mesh_candidate 网格
Polygon_mesh::Property_map<Polygon_mesh::Face_index, const Kernel::Plane_3*> face_supporting_planes
	= mesh_candidate.add_property_map<Polygon_mesh::Face_index, const Kernel::Plane_3*>("f:supp_plane").first;//每个面的支持平面

Polygon_mesh::Property_map<Polygon_mesh::Face_index, CGAL::internal::Planar_segment<Kernel>*> face_supporting_segments
	= mesh_candidate.add_property_map<Polygon_mesh::Face_index, CGAL::internal::Planar_segment<Kernel>*>("f:supp_segment").first;//每个面的支持平面片段

Polygon_mesh::Property_map<Polygon_mesh::Edge_index, std::set<const Kernel::Plane_3*> > edge_supporting_planes
	= mesh_candidate.add_property_map<Polygon_mesh::Edge_index, std::set<const Kernel::Plane_3*> >("e:supp_plane").first;//每个边的所有支持平面

Polygon_mesh::Property_map<Polygon_mesh::Vertex_index, std::set<const Kernel::Plane_3*> > vertex_supporting_planes
	= mesh_candidate.add_property_map<Polygon_mesh::Vertex_index, std::set<const Kernel::Plane_3*> >("v:supp_plane").first;//每个顶点的所有支持平面

这里我额外添加了一个新的属性f:region_map，用于确定面片所属的区域，给以后的面片简化用的，本质上与f:supp_segment近似。

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

计算相交点

接下来，遍历聚类点云的每个区域，将每个区域拟合为切面cutting_plane；

///外部嵌套了一层对区域的遍历
/// for (std::size_t i = 0; i < _cloud_plane->planar_segments().size(); ++i)		
		CGAL::internal::Planar_segment<Kernel>* g = _cloud_plane->planar_segments()[i];
		const Kernel::Plane_3* cutting_plane = _cloud_plane->planar_segments()[i]->supporting_plane();
		std::vector<Kernel::Point_3> intersecting_points;//容器_交叉点
		std::vector< std::set<const Kernel::Plane_3*> > intersecting_points_source_planes;//容器_交叉点的源平面

对这里的每个面，计算其与包围盒的各个交点，也就是确定其边界。具体实现中是遍历包围盒的边，取边的两个端点构成线段，判断该直线与cutting_plane是否存在交点。这里源码有点乱，它分别作了三次交点判断：面与线段st、面分别与s、t两个端点。
我的理解是，它先判断两个端点是否位于切面的两侧，是则构造线段st，求st与cutting_plane的交点；否则，分别判断点s、t是否刚好在cutting_plane上，然后把相交的点存入intersecting_points，并将对应的包围盒支撑平面记录为交叉点的所属平面。

//2、对每个平面，遍历 mesh_bbox 的所有边，判断相交的点并存储
for (auto ed : bbox_mesh.edges())
{
	//获取边的两个顶点
	Polygon_mesh::Vertex_index sd = bbox_mesh.vertex(ed, 0);
	Polygon_mesh::Vertex_index td = bbox_mesh.vertex(ed, 1);
	const Kernel::Point_3& s = bbox_mesh.points()[sd];
	const Kernel::Point_3& t = bbox_mesh.points()[td];
	CGAL::Oriented_side ss = cutting_plane->oriented_side(s);//确定内外侧
	CGAL::Oriented_side st = cutting_plane->oriented_side(t);
	//判断平面 plane_cutting 是否与线段 st 相交
	if ((ss == CGAL::ON_POSITIVE_SIDE && st == CGAL::ON_NEGATIVE_SIDE)
		|| (ss == CGAL::ON_NEGATIVE_SIDE && st == CGAL::ON_POSITIVE_SIDE))
	{
		CGAL::Object obj = CGAL::intersection(*cutting_plane, Kernel::Line_3(s, t));//获取相交点
		if (const Kernel::Point_3* p = CGAL::object_cast<Kernel::Point_3>(&obj))
		{
			intersecting_points.push_back(*p);//相交点添加至 intersecting_points
			std::set<const Kernel::Plane_3*> planes = bbox_edge_supporting_planes[ed];//获取与该相交点相关的所有支持平面
			//平面添加至 planes 、 intersecting_points_source_planes
			planes.insert(cutting_plane);
			CGAL_assertion(planes.size() == 3);
			intersecting_points_source_planes.push_back(planes);
		}
		else
			std::cerr << "错误：平面 plane_cutting 与线段 st 不存在交集" << std::endl;
	}
	else
	{
		//plane_cutting 与线段 s-t 不相交
		//判断 plane_cutting 是否与其中一个顶点相交
		if (ss == CGAL::ON_ORIENTED_BOUNDARY && st != CGAL::ON_ORIENTED_BOUNDARY)
		{
			intersecting_points.push_back(s);
			const std::set<const Kernel::Plane_3*>& planes = bbox_vertex_supporting_planes[sd];
			CGAL_assertion(planes.size() == 3);
			intersecting_points_source_planes.push_back(planes);
		}
		else if (st == CGAL::ON_ORIENTED_BOUNDARY && ss != CGAL::ON_ORIENTED_BOUNDARY)
		{
			intersecting_points.push_back(t);
			const std::set<const Kernel::Plane_3*>& planes = bbox_vertex_supporting_planes[td];
			CGAL_assertion(planes.size() == 3);
			intersecting_points_source_planes.push_back(planes);
		}
	}
}

筛选凸包点

得到切面与点云包围盒的交点后，根据交点坐标构造凸包结构，也就是在切面cutting_plane上确定这些点的连接顺序，以他们为边界生成一块面片。
代码里先做了个合理性判断，size大于3是指面片至少得有3个顶点来构成三角面；然后将这些交点投影至cutting_plane所在的二维坐标系上，调用CGAL库中的convex_hull方法来构成二维凸包；最后，记录凸包点和对应的面片至point_convexHull和ch_source_planes中。

//3、确定 intersecting_points 中的凸包点
if (intersecting_points.size() >= 3)
{
	//将 intersecting_points 的所有点转为二维坐标，存入 pts
	std::list<Kernel::Point_3> pts;
	for (std::size_t i = 0; i < intersecting_points.size(); ++i)
	{
		const Kernel::Point_3& p = intersecting_points[i];
		const Kernel::Point_2& q = cutting_plane->to_2d(p);//返回仿射变换下的投影的图像点，该投影映射到 XY 平面上，并删除 z 坐标
		pts.push_back(Kernel::Point_3(q.x(), q.y(), FT(i))); // z 分量存储点索引
	}
	//计算 pts 的凸包点集，存入 hull
	std::list<Kernel::Point_3> hull;
	CGAL::convex_hull_2(pts.begin(), pts.end(), std::back_inserter(hull), CGAL::Projection_traits_xy_3<Kernel>());
	//将所有凸包点从 intersecting_points 中存入 point_convexHull 和 ch_source_planes 
	std::vector<Kernel::Point_3> point_convexHull;//容器_凸包点
	std::vector< std::set<const Kernel::Plane_3*> > ch_source_planes;//容器_凸包平面
	for (typename std::list<Kernel::Point_3>::iterator it = hull.begin(); it != hull.end(); ++it)
	{
		std::size_t idx = std::size_t(it->z());
		point_convexHull.push_back(intersecting_points[idx]);
		ch_source_planes.push_back(intersecting_points_source_planes[idx]);
	}
	//4、 mesh_candidate 网络创建新面
	//......
}

构造（含边界的）新面

确定出切面边界（凸包）的点坐标后，按顺序连接来构建面片。转存这些点的下标到descriptors后，使用add_face()函数在候选平面集对象mesh_candidate内构造新的面。add_face()的调用示例可参考^[18]，它属于CGAL中的Euler运算^[19]，能在维护原有半边结构的情况下创建多边形网络的新边、面。（注：建议瞅一瞅，里面有些概念图能加深对半边、环、面的理解）
再后面，是对mesh_candidate中刚刚新加入的边赋予支撑面属性。这里作者使用一种围绕面的半边环绕迭代器Halfedge_around_face_circulator来为边ed赋予支撑面属性，即指定这个面片的边是来自哪个切面/拟合的平面系数上。

//4、 mesh_candidate 网络创建新面
if (point_convexHull.size() >= 3)
{
	std::vector<Polygon_mesh::Vertex_index> descriptors;
	for (std::size_t j = 0; j < point_convexHull.size(); ++j)
	{
		Polygon_mesh::Vertex_index vd = mesh_candidate.add_vertex(point_convexHull[j]);
		descriptors.push_back(vd);
		vertex_supporting_planes[vd] = ch_source_planes[j];
		CGAL_assertion(vertex_supporting_planes[vd].size() == 3);
	}
	//将该面上的所有顶点的 supporting_planes 属性设置为 ch_source_planes
	Polygon_mesh::Face_index fd = mesh_candidate.add_face(descriptors);
	face_supporting_segments[fd] = _cloud_plane->planar_segments()[i];
	face_supporting_planes[fd] = cutting_plane;
	//new:为该面添加对应的点云区域下标
	face_region_map[fd] = i;
	//std::cout << "i: " << i << "  region_map_i: " << _cloud_plane->property_map<int>("region_map").second[&i] << std::endl;
	//为该面的每个边缘设置 supporting_planes 属性
	CGAL::Halfedge_around_face_circulator<Polygon_mesh> hbegin(mesh_candidate.halfedge(fd), mesh_candidate), done(hbegin);
	do {
		Polygon_mesh::Halfedge_index hd = *hbegin;
		Polygon_mesh::Edge_index ed = mesh_candidate.edge(hd);
		Polygon_mesh::Vertex_index s_vd = mesh_candidate.source(hd);
		Polygon_mesh::Vertex_index t_vd = mesh_candidate.target(hd);
		const std::set<const Kernel::Plane_3*>& s_planes = vertex_supporting_planes[s_vd];
		const std::set<const Kernel::Plane_3*>& t_planes = vertex_supporting_planes[t_vd];
		std::set<const Kernel::Plane_3*> common_planes;
		std::set_intersection(s_planes.begin(), s_planes.end(), t_planes.begin(), t_planes.end(), std::inserter(common_planes, common_planes.begin()));
		if (common_planes.size() == 2)
		{
			//如果两个顶点的 supporting_planes 属性都包含两个平面，则将这两个平面设置为该边缘的 supporting_planes 属性
			CGAL_assertion(edge_supporting_planes[ed].size() == 0);
			edge_supporting_planes[ed] = common_planes;
			CGAL_assertion(edge_supporting_planes[ed].size() == 2);
		}
		else//如果两个顶点的 supporting_planes 属性不都包含两个平面，则说明发生了拓扑错误
			std::cerr << "topological error 拓扑错误" << std::endl;
		++hbegin;
	} while (hbegin != done);
}

小结

循环每个点云区域并按上述方法构建面片，mesh_candidate内就有若干个跨幅大小为包围盒范围的大块面片。把实际建筑的各个特征区域在脑海里拉长延伸为平面，那此时mesh_candidate中的面片即为这种情况的可视化。
后续的平面相交细分、面片置信度计算以及最优组合求解等步骤，其实都是对这一步生成的面片做细化、剪枝处理。到目前为止已经挺长的了这篇文章，中场休息一下，这些内容就放后面章节再做整理啦。

参考文献&引用

[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