前排瞌睡杂物堆

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

0%

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

发表于 更新于 分类于
本文字数： 23k 阅读时长 ≈ 42 分钟

摘要：
线性整数规划；提取邻接关系；能量项设置及累计；构建方程约束；重建模型。

历经点云聚类、多边形网络重建等处理，手中这花花绿绿的候选面片数据还剩一个主要的处理步骤：线性整数规划。以我从运筹学/最优化理论那粗浅借来的理解，这个过程将在某些硬约束条件下，以最小化能量函数中能量项的加权和为目标，通过求解器的解算来得出无边界、闭合流型的三维模型。

线性规划之胡言乱语

还是更具体地拆解下这个过程中的思路吧。为描述面片选取过程，以集合中的面要素和边要素为变量，构建二元的混合整数线性规划问题，感觉就是把多边形网络中的面和线想象成方程中的变量x、y，设置多个包含不同x和y的方程构成方程组？线性规划就像是求解出其中最想要的（多个）x和y值。
三角锥
怎么形容呢，以上面的一个三角锥几何体为例，这是很久之前我测试ply文件读写功能时绘制的，它由4个面和6条边（/12条半边）构成，相当于这个几何体包含了4个x、6个y。这些变量是已知的全部条件，把它们全部丢给求解器后，为了让它返回一个完整闭合的几何体而不是分散的面/线，就需要设置称为“条件约束”的额外要求。一个较容易理解的约束是要求每条线y必须与两个面x相连，即连续闭合的几何体中每个面都是两两相接的。在当前的这个三角锥中很容易能靠这个约束得到唯一解，毕竟唯一符合该要求的x、y组合只有一种情况。倘若把这个几何体沿中间腰部横切一刀，得到一个小三角锥和梯形体，满足这个约束条件的就有两组解，这时候就要引入“能量函数”或者应该说是“目标函数”的概念，它应该是通过给每个x/y赋予不同的属性值，比方说我期望得到的几何体还是三棱锥，那在求解前我会对所有面x做一次三角面/四角面检测，给“三角面”赋一个较高的属性值，“四角面”赋一个较低的属性值，当解算器筛选出了两组（或更多）解后，进行能量函数求解的过程，将每组xy的属性值累加起来，取值最大的那组解为输出结果。
那么这些便是我对线性规划概念的理解了，不是很严谨肯定是有问题的，欢迎斧正及补充。主要还是自己怠惰了，没认真找过这方面的系统性学习资料，仅看到过的两份基础知识入门[24]、[25]也是一头雾水，以后有兴趣了再翻翻看吧（笑）。

PolyFit的线性规划术语

胡扯了这么久，也该来看看论文中关于这部分的阐述。文中作者设置了一个硬约束来保证模型的闭合性，见下图的闭合模型硬约束公式，第一条的等式左侧变量是所有的边ei连接的面数量，讲人话来表述是它要求所有被选取的边只能连接两个面片，或所连接面片能构成更大的多边形（合并后所接的面数为0），从而使选取的边均位于闭合多边形上而不是外凸面。但如前所述，仅在该约束下候选平面集内存在多种符合条件的组合，为求得最优结果，PolyFit算法基于面片支撑点数量、面片面积和点覆盖面积这三类面片属性构建优化问题的能量函数，并通过能量最小化过程得到最优解。其中的最小化过程可以理解为，每个面片根据上述参数赋予一个值，在筛选出满足约束的所有组合后，取累计值最小的组合为规划问题的解。
硬约束与能量方程
因此，如何定义合理的能量函数也是该方法的关键。PolyFit算法的能量函数由这三个能量项构成：
能量函数项
数据拟合项评估由点簇所拟合平面的置信度，该值根据候选面内的平面点距及点与面的法向量差值来确定，可反映对应点云的平坦度及信噪比。
模型复杂性项定义为模型中锐利边缘的比例。二元函数𝑐𝑜𝑟𝑛𝑒𝑟根据相连面属性决定边e是否为锐边，由于面片均以两两相交分割得到，非孤立边与面片的所有拓扑关系如下图所示，图b、c表示边能将相连面组合为一个更大的多边形，𝑐𝑜𝑟𝑛𝑒𝑟值为1；图d到g中的边可表示为模型边界，𝑐𝑜𝑟𝑛𝑒𝑟值为0。
点覆盖项指模型中未覆盖点云的面片占比，它根据点簇在拟合面上的面积与拟合面总面积比值来确定。
候选面拓扑关系

面片最优重建的整体流程

终于，理论部分基本总结完毕，进入代码解析环节。有关面片最优重建的源码封装在ReconstructProcess文件内，代码的原始版本来自Polygonal_surface_reconstruction/polyfit_example_model_complexity_control.cpp，网址[14]往下翻的3.3章节便能看到，作者列出了好3种能构成不同表面复杂度的指导参数，这些参数是上面提及的能量项系数，我这里设置的值与论文建议值稍微有点不同：数据拟合项wt_fitting = 0.43，模型复杂性项wt_complexity = 0.30，点覆盖项wt_coverage = 0.27，论文建议值是 0.46，0.27，0.27，相比起来相当于把拟合项的权重匀了点给复杂性项，可能是当时的我更看重模型的细节吧，印象中复杂性项的权重越高那期望的模型锐边数目/局部细节会更多。
主程序中对该部分的调用代码如下：

1
2
3
4
5
ReconstructProcess obj_recon;
CGAL::Surface_mesh<Kernel::Point_3> mesh_model;//模型
//obj_recon.UseExtraFactor(0.05f);
//obj_recon.Reconstruct(mesh_candidate, mesh_model, obj_hp.GetTriplet(), 0.43f, 0.27f, 0.25f);//筛选构成模型的面片集合
obj_recon.Reconstruct(mesh_candidate, mesh_model, obj_hp.GetTriplet(), 0.43f, 0.27f, 0.30f);//默认的推荐参数

基于上一大节中的候选面集mesh_candidate，函数Reconstruct()将对其进行邻接信息提取、设置约束条件、构建约束方程和求解方程4步，将得到的最优面集组合存入到model中。

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
bool ReconstructProcess::Reconstruct(Polygon_mesh & mesh_candidate, Polygon_mesh & model,
	const Planes_intersections & map_planes_intersection,
	double wt_fitting, double wt_coverage, double wt_complexity)
{
	CGAL::Timer timer;//时间计数的实例
	timer.start();

	//0、检验数据有效性
	if (mesh_candidate.num_faces() < 4)
	{
		std::cout << "错误：至少需要4个候选表面才能重建出紧密模型，当前的表面集仅有 " + std::to_string(mesh_candidate.num_faces()) + " 个面";
		return false;
	}
	//参数初始化
	InitParameter(wt_fitting, wt_coverage, wt_complexity);


	//1、计算表面的邻接信息：mesh_candidate的三重交点及其边的所在面
	GetAdjacency(mesh_candidate, map_planes_intersection);


	//2、设置约束条件
	if (_isUseRoofObj)
	{
		DefineConstraintFactors_v2(mesh_candidate);//各能量项归一化处理
		DefineConstraintFactor_TopDis(mesh_candidate);//添加顶面偏好项
	}
	else
		DefineConstraintFactors(mesh_candidate);//原PolyFit方法
	///输出测试_能量项累计值
	//CheckObjectiveValue(mesh_candidate);


	//3、构建约束方程
	BuildConstraint(mesh_candidate);


	//4、求解线性约束整数方程
	SolveConstraint(mesh_candidate, model);

	timer.stop();
	std::cout << "重建模型的面数为 " << model.faces().size() << " 个，运行耗时 " << timer.time() << " 秒" << std::endl;

	return true;
}

邻接信息提取

第一步是把候选面集的半边分类为交叉边和边缘边，便于确定所选边、面的关系。
该小节的代码实现在函数ReconstructProcess::GetAdjacency()，其中2部分的内容没太弄懂其中逻辑，大概是先取候选面集内的每个半边，依据前一节的三重交点集合map_planes_intersection快速定位出两个端点的坐标，以两个端点坐标为键、半边下标为值存入到哈希表face_pool。但是，我要说但是，代码里写了很多断言CGAL_assertion()，看似是做了些判断筛选，但这些都是debug用的，除非前面有些步骤导致数据出现缺失不然都不会触发，真正做了判断的仅两个：去除边界半边（只与一个面相接）、调整端点顺序，况且已知mesh中点下标的情况下可以直接拿取对应的点坐标，感觉根本不需要传入map_planes_intersection来获取。不过倒也体现了CGAL的严密性，如果前面搞了些删属性的骚操作到这一步就会跑不通。

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
//1、定义、声明
typename Polygon_mesh::template Property_map<Polygon_mesh::Vertex_index, std::set<const Kernel::Plane_3*> > vertex_supporting_planes
	= mesh_candidate.template property_map<Polygon_mesh::Vertex_index, std::set<const Kernel::Plane_3*> >("v:supp_plane").first;//属性映射_每个点的支撑面

// 边由它的两个端点表示
typedef typename std::unordered_map<const Kernel::Point_3*, std::set<Polygon_mesh::Halfedge_index> >	Edge_map;//无序映射_点指针对应半边索引集合，表示边与它的半边
typedef typename std::unordered_map<const Kernel::Point_3*, Edge_map >									Face_pool;//无序映射_点指针对应Edge_map，表示一组邻接关系
Face_pool face_pool;


//2、向 face_pool 中存入非边界边的所在面及边的两个顶点坐标，每个顶点为三重交点
for (auto h : mesh_candidate.halfedges())
{
	Polygon_mesh::Face_index f = mesh_candidate.face(h);//获取每个半边所在的面索引
	if (f == Polygon_mesh::null_face())
		continue;//半边没有对应的面，则为边界半边，跳过

	//获取半边的目标顶点+源顶点、其所在的平面
	Polygon_mesh::Vertex_index sd = mesh_candidate.source(h);
	Polygon_mesh::Vertex_index td = mesh_candidate.target(h);
	const std::set<const Kernel::Plane_3*>& set_s = vertex_supporting_planes[sd];
	const std::set<const Kernel::Plane_3*>& set_t = vertex_supporting_planes[td];
	CGAL_assertion(set_s.size() == 3);//断言_点所在的平面集合有三个元素，即顶点是由三个平面相交而成的
	CGAL_assertion(set_t.size() == 3);

	std::vector<const Kernel::Plane_3*> s_planes(set_s.begin(), set_s.end());//将平面集合转存至容器中，方便调取
	CGAL_assertion(s_planes[0] < s_planes[1]);//断言_平面按序排列
	CGAL_assertion(s_planes[1] < s_planes[2]);
	const Kernel::Point_3* s = map_planes_intersection[s_planes[0]][s_planes[1]][s_planes[2]];//获取交点坐标

	std::vector<const Kernel::Plane_3*> t_planes(set_t.begin(), set_t.end());
	CGAL_assertion(t_planes[0] < t_planes[1]);
	CGAL_assertion(t_planes[1] < t_planes[2]);
	const Kernel::Point_3* t = map_planes_intersection[t_planes[0]][t_planes[1]][t_planes[2]];

	if (s > t)
		std::swap(s, t);//源顶点坐标大于目标顶点，则它们的顺序不一致，交换两者的值
	face_pool[s][t].insert(mesh_candidate.halfedge(f));//源顶点和目标顶点坐标为键，半边所在的面索引为值，插入到 face_pool
}

第3部分则把face_pool换了个结构体Intersection转存到std::vector _adjacency。

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
//结构体_平面相交半边的点结构，继承自容器_多边形网络半边索引
struct Intersection : public std::vector<Polygon_mesh::Halfedge_index>
{
	const Kernel::Point_3* s;//指针_相交半边的起点坐标
	const Kernel::Point_3* t;//指针_相交半边的终点坐标
};


//3、使用 _adjacency 转存 face_pool 中的信息
for (Face_pool::const_iterator it = face_pool.begin(); it != face_pool.end(); ++it)
{
	const Kernel::Point_3* s = it->first;//源顶点坐标
	const Edge_map& tmp = it->second;
	typename Edge_map::const_iterator cur = tmp.begin();
	for (; cur != tmp.end(); ++cur)
	{
		const Kernel::Point_3* t = cur->first;//目标顶点坐标
		const std::set<Polygon_mesh::Halfedge_index>& faces = cur->second;//集合_半边索引
		Intersection fan;
		fan.s = s;
		fan.t = t;
		fan.insert(fan.end(), faces.begin(), faces.end());
		_adjacency.push_back(fan);
	}
}

设置约束条件

先看DefineConstraintFactors()中的内容，这里面封装的是原PolyFit的约束方案。第1部分是常规的变量声明，提一嘴题外话，更老练的变量声明应该放在头文件中类的私有变量那，函数用到的时候再给它定义/初始化/分配内存，印象中PolyFit的源码也是这么写的，当初大概是为了理清算法里面的过程，把这些变量的声明定义绑到一起塞到各个函数内做封装，现在回看有点吃力不讨好。

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
//1、参数定义/起别名
//mesh_candidate属性映射
typename Polygon_mesh::template Property_map<Polygon_mesh::Face_index, std::size_t> face_num_supporting_points =
	mesh_candidate.template property_map<Polygon_mesh::Face_index, std::size_t>("f:num_supporting_points").first;//属性映射_面的支撑点数量

typename Polygon_mesh::template Property_map<Polygon_mesh::Face_index, FT> face_areas =
	mesh_candidate.template property_map<Polygon_mesh::Face_index, FT>("f:face_area").first;//属性映射_面的面积

typename Polygon_mesh::template Property_map<Polygon_mesh::Face_index, FT> face_covered_areas =
	mesh_candidate.template property_map<Polygon_mesh::Face_index, FT>("f:covered_area").first;//属性映射_面的点覆盖面积

typename Polygon_mesh::template Property_map<Polygon_mesh::Face_index, std::size_t> face_indices =
	mesh_candidate.template add_property_map<Polygon_mesh::Face_index, std::size_t>("f:index").first;//属性映射_面索引

//支撑点总数
double total_points = 0.0;
std::size_t idx = 0;
for (auto f : mesh_candidate.faces())
{
	total_points += face_num_supporting_points[f];
	face_indices[f] = idx;//给每个面附上新的索引值
	++idx;
}

std::size_t num_faces = mesh_candidate.number_of_faces();//候选表面的面数量
std::size_t num_edges(0);//交叉边数量，初始值为0

第2部分，计算一些中间变量，主要是为交叉边索引edge_usage_status和锐利边索引edge_sharp_status取值，以此设置线性规划求解器MIP_Solver的变量对象variables。这里作者设置约束的变量总个数为候选面数+两倍交叉边数，每个变量为二进制类型只可取值0或1。这段代码对应图2硬约束的第二个式子，其中，每个变量取二进制好理解，1和0表示在提取最优面集的过程中，是否保留该候选面/边，比较摸不着头脑的是它的变量总数是num_faces + num_edges * 2 而不是 num_faces + num_edges，问了下ai猜测原因如下：
（1）每个候选面对应一个二进制变量（数量为num_faces），用于决定该面是否被选中。
（2）每个边需要两个变量进行约束，第一个为保留状态变量（数量为num_edges），决定该边是否保留；第二个为锐/交叉边状态变量（数量为num_edges），用于判断该边是否为交叉边。

1
2
3
private:
	std::unordered_map<const Intersection*, std::size_t> edge_usage_status;//无序映射_交叉边索引，表示该交叉边被保留或移除
	std::unordered_map<const Intersection*, std::size_t> edge_sharp_status;//无序映射_锐利边索引，表示其是否为锐利边
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
//2、计算约束条件的值
// edge_usage_status 赋值
for (std::size_t i = 0; i < _adjacency.size(); ++i)
{
	//遍历所有含邻接关系的边，并将包含4个半边的交叉边索引添加到 edge_usage_status
	//该边由4个候选面相交而成
	const Intersection& fan = _adjacency[i];
	if (fan.size() == 4)
	{
		//交叉边索引 为键，候选面数量+交叉边数量 为值
		std::size_t var_idx = num_faces + num_edges;
		edge_usage_status[&fan] = var_idx;
		++num_edges;
	}
}

// variables 定义、赋值
std::size_t total_variables = num_faces + num_edges + num_edges;//总变量数量 = 候选面数量 + 2 * 交叉边数量
variables = solver.create_variables(total_variables);//容器_变量个数为 total_variables
for (std::size_t i = 0; i < total_variables; ++i)
{
	//为每个变量设置类型为 BINARY二进制，只能取0或1
	MIP_Solver::Variable* v = variables[i];
	v->set_variable_type(MIP_Solver::Variable::BINARY);
}

第3部分，设置线性目标函数对象objective，它的类型（目标）为线性最小化，这意味着在整个混合整数规划（MIP）问题中，求解器将会寻找一组变量取值，使得目标函数的值尽可能低。这块值得注意的点是对权重系数的再次计算，在我摘抄的代码中拟合项的权重不变，余下两项乘上点云的总点数，其中点覆盖项再除以包围盒表面积（area(M)），复杂度项除以非边界的半边数量（|E|）。它们俩的除数处理对应能量函数项公式中各自的分子项；乘以总点数的理由按作者在GitHub上的答疑[26]，是为了将权重系数从[0,1.0f]的范围映射到[0,total_points]，以避免浮点数的精度损失。
（注：我不太确定当时摘抄是否改动了权重系数的内容，应该是没有的没那么扎实的理论储备我不敢乱调，且github上对形类似式的权重系数代码，作者的答疑是“所有三个权重都简单地乘以总点数”，或许拟合项的权重系数_wt_fitting在别的地方做了“乘以总点数”的处理）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
//3、添加MIP_Solver的 objective（各个能量项？）
//候选面集的顶点坐标及其包围盒参数
std::vector<Kernel::Point_3> vertices(mesh_candidate.number_of_vertices());
idx = 0;
for (auto v : mesh_candidate.vertices())
{
	vertices[idx] = mesh_candidate.points()[v];
	++idx;
}
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();
FT box_area = FT(2.0) * (dx * dy + dy * dz + dz * dx);//包围盒表面积

//选择更好的尺度
double coeff_data_fitting = _wt_fitting;//拟合项的权重系数
double coeff_coverage = total_points * _wt_coverage / box_area;//点覆盖项的权重系数
double coeff_complexity = total_points * _wt_complexity / double(_adjacency.size());//复杂度项的权重系数

///输出测试
//std::cout << coeff_data_fitting << "  " << coeff_coverage << "  " << coeff_complexity << std::endl;

MIP_Solver::Linear_objective * objective = solver.create_objective(MIP_Solver::Linear_objective::MINIMIZE);//指针_线性目标函数对象，并设置类型为最小化

余下的一小部分内容，是对线性目标函数objective中的每个变量variables累计能量项值，然后设置约束条件MIP_Solver::Linear_constraint* c。代码中每个variables累计能量项的过程为：找到非边界半边集合_adjacency中的所有锐边，属于锐边的变量累加一个复杂度项的值coeff_complexity；遍历候选面集mesh_candidate中的每个面片，按面下标找到对应的variables并累计数据拟合项的值：
负的复杂度项权重*该面片的整成点数量。支撑点数量的计算过程可跳转到compute_confidences.h的Candidate_confidences::compute查阅（201-202行），负号意味着支持点越多，拟合能量越低，从而鼓励目标函数选取支持点数更多的面；
最后一个点覆盖项，先计算每个面的未覆盖面积（面片的总面积减去已覆盖面积），并以coeff_coverage为权重将该值累计进目标函数的MIP变量。这一项的目的是惩罚未被覆盖的区域，从而引导最终模型更好地覆盖点云数据。

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
// edge_sharp_status 赋值
std::size_t num_sharp_edges = 0;
for (std::size_t i = 0; i < _adjacency.size(); ++i)
{
	const Intersection& fan = _adjacency[i];
	if (fan.size() == 4)
	{
		std::size_t var_idx = num_faces + num_edges + num_sharp_edges;//尖锐边数量可能小于交叉边，部分variables全程0值？
		edge_sharp_status[&fan] = var_idx;

		//累计复杂度项
		objective->add_coefficient(variables[var_idx], coeff_complexity);//为每个MIP变量添加系数
																		 //系数 = 复杂度项的权重，表示复杂度项是最小化目标
		++num_sharp_edges;
	}
}
CGAL_assertion(num_edges == num_sharp_edges);

for (auto f : mesh_candidate.faces())
{
	std::size_t var_idx = face_indices[f];

	//累计数据拟合项
	std::size_t num = face_num_supporting_points[f];//面的支撑点数量
	objective->add_coefficient(variables[var_idx], -coeff_data_fitting * num);//为每个MIP变量添加系数
																			  //系数=拟合能量项的权重*（-支撑点数），表示拟合能量项是最小化目标
	//累计点覆盖项
	double uncovered_area = (face_areas[f] - face_covered_areas[f]);//面的未覆盖面积
	objective->add_coefficient(variables[var_idx], coeff_coverage * uncovered_area);
}

最后为了满足闭合网格要求（内部边恰好被两个面共享），添加约束条件MIP_Solver::Linear_constraint* c确保“与边关联的面数必须为2或0”。在面变量部分，对fan中每个候选面，将对应的面变量以系数1添加到约束中，这样约束项累加了fan内所有面的变量；边变量部分，若为内部交叉边（fan.size() == 4，fan由四个面构成），取出对应的边变量并添加系数-2。由此构成的约束表达式如下所示：
闭合硬约束
由于边界边不属于封闭模型的一部分，若该变量为边界边，则不对其累计额外的约束值。（注：这里add_coefficient的对象为线性约束变量c，上面的add_coefficient对象为线性目标函数objective，注意累计值的区分）

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
//添加约束：与边关联的面数必须为2或0
std::size_t var_edge_used_idx = 0;//记录当前交叉边对应的变量索引
for (std::size_t i = 0; i < _adjacency.size(); ++i)
{
	MIP_Solver::Linear_constraint* c = solver.create_constraint(0.0, 0.0);//线性约束_约束上界和下界 = 0，表示无约束条件约束(无取值范围限制
	const Intersection& fan = _adjacency[i];
	for (std::size_t j = 0; j < fan.size(); ++j)
	{
		//对面添加约束
		Polygon_mesh::Face_index f = mesh_candidate.face(fan[j]);
		std::size_t var_idx = face_indices[f];
		c->add_coefficient(variables[var_idx], 1.0);//表示约束中包含变量 variables[var_idx]
	}

	//交叉边由四个候选面相交而成
	if (fan.size() == 4)
	{
		//对交叉边添加约束
		std::size_t var_idx = num_faces + var_edge_used_idx;
		c->add_coefficient(variables[var_idx], -2.0);//表示约束中包含这个变量，并且它与其他变量之和必须等于零？
		++var_edge_used_idx;
	}
	else
	{
		//交叉边为边界边，不添加任何系数，因为生成的模型不为开放的表面模型
	}
}

设置额外的能量项

首要说明的是，代码中的DefineConstraintFactors_v2()与DefineConstraintFactor_TopDis()现在看来均是有问题的。DefineConstraintFactors_v2()相较于DefineConstraintFactors()少了尺度处理部分，本意是想直接按归一化的权重来累计objective中的能量值，但佘略的步骤有两个隐患：一是代码与能量项公式不符，每个项的系数均少了分母项；二是如[26]所说会有潜在的精度损失。而另一个函数DefineConstraintFactor_TopDis()是我对修改能量公式的一次失败尝试，从最终结果上看，顶部表面凹凸不平坑坑洼洼而不是更丰富、规整的细化结构，或许是权重分配的不合理/未设置额外的条件约束。
不管怎么说，还是把这一坨保留下来存档了，属于是相当简陋了（笑）。代码中先取出指向MIP求解器的线性对象指针*objective，为mesh_candidate中的所有面对象累计值_wt_roof * dist_p2boxTop，_wt_roof是预设的权重系数，dist_p2boxTop是该面片到包围盒顶部的归一化距离值，详见上一大节的“添加额外的平面属性”，理论上会让求解器更倾向于选取靠近顶部的面（这类variables累计的值相对而言更小）。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
void ReconstructProcess::DefineConstraintFactor_TopDis(const Polygon_mesh & mesh_candidate)
{
	typename Polygon_mesh::template Property_map<Polygon_mesh::Face_index, double> face_roof =
		mesh_candidate.template property_map<Polygon_mesh::Face_index, double>("f:dist_boxTop").first;//属性映射_面片属顶面

	typename Polygon_mesh::template Property_map<Polygon_mesh::Face_index, std::size_t> face_indices =
		mesh_candidate.template property_map<Polygon_mesh::Face_index, std::size_t>("f:index").first;//属性映射_面索引


	//double total_points = 0.0;//支撑点总数
	//std::size_t num_faces = mesh_candidate.number_of_faces();//候选表面的面数量

	MIP_Solver::Linear_objective * objective = solver.objective();

	//对每个能量函数，累计顶面偏好值
	for (auto f : mesh_candidate.faces())
	{
		std::size_t var_idx = face_indices[f];

		double dist_p2boxTop = face_roof[f];
		objective->add_coefficient(variables[var_idx], _wt_roof * dist_p2boxTop);//离包围盒顶面越远，值越大
	}
}

这部分默认是不使用的，可将_isUseRoofObj修改为true来启用。

构建约束方程

这一部分的内容那会翻译错名字了，应该是“补充约束”更贴合些，其核心目的是为候选面片mesh_candidate中那些与锐边相关的对象添加约束，从而保证在重构过程中如果某个边被认定为锐边，那么其相邻的面必须满足一定的选择条件，避免出现不连续或不合理的表面。
代码封装在ReconstructProcess::BuildConstraint()，代码解析大多借助GPT。松弛因子M用于构造带有条件逻辑的线性约束，使得约束可以用线性表达式表示。然后是按照不等式“X[var_edge_usage_idx] - X[var_edge_sharp_idx] >= 0”构建约束，它表明如果一个边被标记为锐边（X[var_edge_usage_idx] = 1），则它一定要被选中（X[var_edge_sharp_idx] = 1），以保持约束成立。若边没有被使用，则无论其是否为锐边，这个约束都会自动满足（可能此时约束的值为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
void ReconstructProcess::BuildConstraint(const Polygon_mesh & mesh_candidate)
{
	//为尖锐的边缘添加约束。提出这个限制的解释可以在这里找到:
	// https://user-images.githubusercontent.com/15526536/30185644-12085a9c-942b-11e7-831d-290dd2a4d50c.png


	// mesh_candidate 属性映射
	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;//属性映射_面的支撑面

	typename Polygon_mesh::template Property_map<Polygon_mesh::Face_index, std::size_t> face_indices =
		mesh_candidate.template property_map<Polygon_mesh::Face_index, std::size_t>("f:index").first;//属性映射_面索引


	double M = 1.0;//线性规划中的松弛变量
	for (std::size_t i = 0; i < _adjacency.size(); ++i)
	{
		const Intersection& fan = _adjacency[i];//获取容器_邻接关系中的交叉边
		if (fan.size() != 4)
			continue;

		//构建约束：X[var_edge_usage_idx] >= X[var_edge_sharp_idx]。表示如果一个交叉边是锐利的，那么它必须被选择
		MIP_Solver::Linear_constraint* c = solver.create_constraint(0.0);//线性约束对象_下界，表示约束是大于等于0的不等式约束
		std::size_t var_edge_usage_idx = edge_usage_status[&fan];
		c->add_coefficient(variables[var_edge_usage_idx], 1.0);//对使用状态变量添加系数，系数值为1.0
		std::size_t var_edge_sharp_idx = edge_sharp_status[&fan];
		c->add_coefficient(variables[var_edge_sharp_idx], -1.0);//对锐利状态变量添加系数，系数值为-1.0

		//额外约束
		//......
	}

	///输出测试
	//std::cout << "最后的各个变量项" << std::endl;
	//CheckObjectiveValue(mesh_candidate);

}

内层的双重循环是为每对邻接关系的相关面添加额外约束，取当前邻接fan关联的所有面两两组合，外层循环的为f1内层循环为f2，判断组合的这两个面是否共面，若对应的支撑面face_supporting_planes不同则不共面，则两个面在几何上存在不连续的可能，需进一步设置约束。该约束的表达式见注释部分，其中，为满足锐边的约束（X[var_edge_sharp_idx] = 1的情况），要求该边相关的四个面至少有两个面被选中。代码的写法应该是在该约束c中，所选的边、面变量的累加值不小于1 - 3M。

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
//额外约束
for (std::size_t j = 0; j < fan.size(); ++j)
{
	//获取相邻次序的两个交叉边及其支撑面
	Polygon_mesh::Face_index f1 = mesh_candidate.face(fan[j]);
	const Kernel::Plane_3* plane1 = face_supporting_planes[f1];
	std::size_t fid1 = face_indices[f1];
	for (std::size_t k = j + 1; k < fan.size(); ++k)
	{
		Polygon_mesh::Face_index f2 = mesh_candidate.face(fan[k]);
		const Kernel::Plane_3* plane2 = face_supporting_planes[f2];
		std::size_t fid2 = face_indices[f2];

		if (plane1 != plane2)
		{
			//两个支撑面不共面的情况下，构成约束：X[var_edge_sharp_idx] + M * (3 - (X[fid1] + X[fid2] + X[var_edge_usage_idx])) >= 1，
			//等价于：X[var_edge_sharp_idx] - M * X[fid1] - M * X[fid2] - M * X[var_edge_usage_idx] >= 1 - 3M
			//表示如果一个交叉边是锐利的，那么它所在的四个面中至少有两个被选择，否则就会导致表面不连续
			c = solver.create_constraint(1.0 - 3.0 * M);//线性约束对象_下界，表示约束是大于等于这个值的不等式约束
			c->add_coefficient(variables[var_edge_sharp_idx], 1.0);//线性约束对象_对锐利状态变量添加一个系数，设为1.0
			c->add_coefficient(variables[fid1], -M);//线性约束对象_设置第一面的变量系数为-M
			c->add_coefficient(variables[fid2], -M);//线性约束对象_设置第二面的变量系数为-M
			c->add_coefficient(variables[var_edge_usage_idx], -M);
		}
	}
}

综上，BuildConstraint()额外设置了两条约束，第一条约束：若一条边为锐边，则该边必须被激活使用；第二条约束：锐边的邻近关系中必须有两个面被选择，来确保网络的表面整体连续。

约束方程求解

第四部分，在列好线性规划问题的变量对象及约束条件后，调用MIP求解器解出候选网络的最优组合，进而生成三维模型。实现代码封装在ReconstructProcess::SolveConstraint()中，首先，先调用求解器的solve()方法，只有在求解器成功找到一个满足所有约束的解时，才对候选网格进行处理。

1
2
3
4
5
6
7
if (solver.solve())//求解线性规划问题成功
{
	//后续的处理基于求解成功后的解
	//...
}
else
	std::cerr << "solving the binary program failed" << std::endl;

求解完毕后，将所有变量（面、边的选择、尖锐边状态）的值用solver.solution()取出并存入到解向量X中，这些解应该是0或1的二元解。然后，
遍历候选面集的下标并从X中取相应的值，若其值四舍五入后为0，说明这个面在解中没有被选中，将其记录到删除列表to_delete。接着调用CGAL::Euler::remove_face()方法，根据to_delete中面的半边索引移除mesh_candidate所有未被选中的面。

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
//1、根据solve中的解，删除 mesh_candidate 中相应的边
typename Polygon_mesh::template Property_map<Polygon_mesh::Face_index, std::size_t> face_indices =
	mesh_candidate.template property_map<Polygon_mesh::Face_index, std::size_t>("f:index").first;//属性映射_面索引

const std::vector<double>& X = solver.solution();//线性规划问题的解

std::vector<Polygon_mesh::Face_index> to_delete;//容器_面索引_记录需要删除的面
std::size_t f_idx(0);
for (auto f : mesh_candidate.faces())
{
	///输出测试2.9
	//std::cout << X[f_idx] << "  ";

	//解的值四舍五入后等于0，说明这个面没有被选择，将该面的索引存入 to_delete
	if (static_cast<int>(std::round(X[f_idx])) == 0)
		to_delete.push_back(f);
	++f_idx;
}

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

还有一小部分内容是更新面集的锐边信息，对mesh_candidate中的边增加一个属性”e:sharp_edges”，然后标记MIP解中满足锐边条件的边，设该属性值为true，方便后续处理中对这些边做特殊处理（提取边界特征、细化锐角等）。本程序后续未有效利用该属性，可视情况注释掉。

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
//2、设置 mesh_candidate 中的锐利边属性
typename Polygon_mesh::template Property_map<Polygon_mesh::Edge_index, bool> edge_is_sharp =
	mesh_candidate.template add_property_map<Polygon_mesh::Edge_index, bool>("e:sharp_edges").first;//属性映射_是否锐利边
for (auto e : mesh_candidate.edges())
	edge_is_sharp[e] = false;//该属性映射默认没有锐利边

for (std::size_t i = 0; i < _adjacency.size(); ++i)
{
	const Intersection& fan = _adjacency[i];
	if (fan.size() != 4)
		continue;//跳过边界边

	//获取交叉边的锐利状态变量索引，若该值为1（该边是锐利的）
	std::size_t idx_sharp_var = edge_sharp_status[&fan];
	if (static_cast<int>(X[idx_sharp_var]) == 1)
	{
		for (std::size_t j = 0; j < fan.size(); ++j)
		{
			//获取每个半边的索引和相应的面索引，若面索引非空（该面未被删除）
			Polygon_mesh::Halfedge_index h = fan[j];
			Polygon_mesh::Face_index f = mesh_candidate.face(h);
			if (f != Polygon_mesh::null_face())
			{
				//获取该面，若果解的fid元素四舍五入后值为1
				std::size_t fid = face_indices[f];
				if (static_cast<int>(std::round(X[fid])) == 1)
				{
					Polygon_mesh::Edge_index e = mesh_candidate.edge(h);
					edge_is_sharp[e] = true;//该面被选择，设置该边的 属性映射_锐利 值为真
					break;
				}
			}
		}
	}
}

此时的mesh_candidate便是最终的网络模型，调用CGAL::copy_face_graph将其复制到model中，它反映了求解器给出的最佳选择结果。
候选面集2闭合模型
OK兄弟们，PolyFit的整个流程总算走完一遍了，传入的点云数据历经区域聚类、平面拟合、两两相交细分、线性规划求解后，终于生成了一整块闭合且贴近原始表面的多边型网络，真是可喜可贺可口可乐~~不过现在它这个花花绿绿的样子还有待简化合并下，下一节再整理整理模型的面片合并和模型导出方法。

参考文献&引用

[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
[24] https://www.gurobi.com/resources/mixed-integer-programming-mip-a-primer-on-the-basics/
[25] https://docs.mosek.com/modeling-cookbook/mio.html
[26] https://github.com/LiangliangNan/PolyFit/issues/30