递归最小边界矩形算法

参考文献：Eunju Kwak,Ayman Habib. Automatic representation and reconstruction of DBM from LiDAR data using Recursive Minimum Bounding Rectangle[J].ISPRS Journal of Photogrammetry and Remote Sensing

该算法效果效果如下图：

主要步骤：

 1. 对点云进行预处理，分离出建筑物点云。

2. 计算点云的平均点距，并根据平均点距设置算法的阈值。

3. 提取建筑物的边界点云。直接用PCL::Concavehull提取到单体建筑物的轮廓。如果存在边界断裂，需要进行连接。

4. 计算建筑物边界点云的主要方向。

5. 将边界点云旋转到主方向位置。

6. 找到旋转后点云的x坐标最小值，以该点作为新的起点将边界点云重新排序。

7. 将旋转后的点云order_cloud的端点存入端点列表epList的第一个位置，计算第一级点云的，将其放入MBR列表mbrList，将epList和mbrList返回的索引放入索引队列中。

8. 如果队列不为空，循环执行步骤（9）至步骤（10），否则执行步骤（11）。

9. 从队列中取出一个元素，根据MBR索引到mbrList找到一个MBR（记作curMBR），通过端点索引到端点列表找到点云。根据阈值，计算点云与curMBR的重叠点和非重叠点，将点云划分为多个次级点云，用一个新的端点列表childList存放。

10. 遍历新的端点列表childList。如果点云满足最小建筑物（长和宽不要过小）的条件，计算每个点云的MBR，以及点云端点在curMBR上的投影点，用投影方向校正算法计算得到正确的MBR。将点云放入epList，把修正后的MBR存入mbrList，把返回的索引存入索引队列。

11. 计算规则化轮廓。定义一个空的边列表edgeList，遍历mbrList，取出每个MBR，对每个MBR的四条边进行多边形的拓扑运算，将得到的边存入edgeList。最终得到规则化轮廓，但是，此时轮廓在edgeList中是以边的形式存放，且是无序存放的，存在冗余的数据，即重复的端点。

具体代码如下：

void RecursiveMinimumBoundingRectangle(const PointCloudT::Ptr& in_cloud, std::vector<PointT>& outline, float threshold)
{
	if (!isValidPointCloud(in_cloud, 0, in_cloud->size()))
  {
		return;
  }
	setDelta(10);
	setInputCloud(in_cloud);
	compute();//计算主方向
	PointCloudT::Ptr out_cloud(new PointCloudT);
	rotatePointCloud(in_cloud, out_cloud, this->min_angle);//旋转到主方向
	int size = out_cloud->size(), _strat = 0;
	float minx = out_cloud->points[0].x;
	for (int i = 1; i < size; ++i)
	{
		if (minx > out_cloud->points[i].x)
		{
			minx = out_cloud->points[i].x;
			_strat = i;
		}
	}
	PointCloudT::Ptr order_cloud(new PointCloudT);
	order_cloud->resize(size);
  //以最小值为起点，重新排序
	for (int i = 0; i < size; ++i)
  {
    order_cloud->points[i] = out_cloud->points[(_strat + i) % size];
  }
  //第一级MBR和边界点云信息
	Endpoint first_cloud = { 0, size };
	Rect first_level = getMBR(order_cloud, first_cloud.start, first_cloud.end);
	std::vector<Rect> all_levels_rectangles = { first_level };
	std::queue<RectAndCloud> que;
	que.push({ 0, -1,-1,first_cloud });
	while (!que.empty())
	{
		RectAndCloud front = que.front();
		Rect previous = all_levels_rectangles[front.mbr_index];
		int proj_first = front.proj_first, proj_second = front.proj_secend;
		que.pop();
		std::vector<Endpoint> subsets_cloud = dividePointCloud(order_cloud, previous, front.cloud.start, front.cloud.end, 2 * threshold);
		for (Endpoint ep : subsets_cloud)
		{
			if (ep.end - ep.start > 1) {
				int proj1, proj2;
				std::vector<PointT> projection = getProjectionPoints(order_cloud, previous, ep.start, ep.end, proj_first, proj_second, proj1, proj2);
				Rect m = getMBR(order_cloud, ep.start, ep.end);
				Rect next = fixMBR(order_cloud, projection, ep.start, ep.end);
				all_levels_rectangles.emplace_back(next);
				int index = all_levels_rectangles.size() - 1;
				que.push({ index,proj1, proj2, ep });
			}
		}
	}
	std::vector<Edge> edges;
	for (Rect& m : all_levels_rectangles)
	{
		topologicalOperations(edges, m);//多边形拓扑运算
	}
	outline = edge2Points(edges);//将边存储转为点存储
  //将结果旋转，保证和属于的点云方向一致。
	float theta = pcl::deg2rad(-this->min_angle);
	float cos_angle = cos(theta);
	float sin_angle = sin(theta);
	float z = in_cloud->points[0].z;
	for (int i = 0; i < outline.size(); ++i)
	{
		outline[i] = rotate(outline[i], cos_angle, sin_angle);
		outline[i].z = z;
	}
	return;
}

拟合效果 ：