【CGAL】Region_Growing 检测平面并保存

说明

本篇博客主要介绍CGAL库中使用Region_Growing算法检测平面的算法原理、代码以及最后展示结果。其中，代码部分在CGAL官方库中提供了例子。我在其中做了一些修改，使其可以读取PLY类型的点云文件，并在检测到平面后，为属于同一个平面的点赋予相同的颜色。最后再保存为PLY文件以方便我们查看检测结果。

在CGAL中，Region_Growing算法不仅可以用来检测平面，还可以检测圆、直线、圆锥等基本的几何。除此之外，用户也可以自定义模型并使用算法检测。

环境

Win10/Win11
VS2022
CGAL 5.6.1

上述环境仅为运行此代码时的电脑环境。

一、算法原理

Region_Growing算法应用“贪心”的思想，利用种子点与邻居点的曲率差异来筛选点。在平面检测时，除了利用曲率差异，还会使用当前拟合平面与当前点的距离作为评判标准。具体看算法流程。

算法流程

选取种子点，若不指定种子点，则按索引顺序选取。
搜索种子点的邻居，可以按球形范围搜索，也可以按邻居个数搜索。
包含满足条件的邻居点。
在包含点中重复1-3。
如果区域内点数不再增加，并且还有未分类的点，则在未分类点中重新选取种子点。

参数

search_sphere_radius - 球形邻域搜索半径
max_distance_to_plane - 当前点到拟合平面的最大距离
max_accepted_angle - 当前点与拟合平面的角度阈值
min_region_size - 形成一个平面的最小点数

二、代码展示

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/IO/read_points.h>
#include <CGAL/property_map.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Shape_detection/Region_growing/Region_growing.h>
#include <CGAL/Shape_detection/Region_growing/Point_set.h>
#include <CGAL/Polygonal_surface_reconstruction.h>
#include <CGAL/IO/write_ply_points.h>#include <fstream>
#include <CGAL/Timer.h>
#include <boost/range/irange.hpp>
typedef CGAL::Exact_predicates_inexact_constructions_kernel        Kernel;
typedef Kernel::FT       FT;
typedef Kernel::Point_3         Point;
typedef Kernel::Vector_3 Vector;
// Point with normal, and plane index.
typedef boost::tuple<Point, Vector, int> PNI;
typedef std::vector<PNI> Point_vector;
typedef CGAL::Nth_of_tuple_property_map<0, PNI>        Point_map;
typedef CGAL::Nth_of_tuple_property_map<1, PNI>        Normal_map;
typedef CGAL::Nth_of_tuple_property_map<2, PNI>        Plane_index_map;
using Point_map_region_growing = CGAL::Compose_property_map<CGAL::Random_access_property_map<Point_vector>, Point_map >;
using Normal_map_region_growing = CGAL::Compose_property_map<CGAL::Random_access_property_map<Point_vector>, Normal_map >;
using Region_type = CGAL::Shape_detection::Point_set::Least_squares_plane_fit_region<Kernel, std::size_t, Point_map_region_growing, Normal_map_region_growing>;
using Neighbor_query = CGAL::Shape_detection::Point_set::Sphere_neighbor_query<Kernel, std::size_t, Point_map_region_growing>;
using Region_growing = CGAL::Shape_detection::Region_growing<Neighbor_query, Region_type>;//----------------------------------save_points_with_color_about-------------------------------------------
typedef std::array<unsigned char, 4> Color;
typedef std::tuple<Point, Vector, Color> PNC;
typedef CGAL::Nth_of_tuple_property_map<0, PNC> Save_Point_map;
typedef CGAL::Nth_of_tuple_property_map<1, PNC> Save_Normal_map;
typedef CGAL::Nth_of_tuple_property_map<2, PNC> Save_Color_map;
//----------------------------------save_points_with_color_about-------------------------------------------// Define how a color should be stored
namespace CGAL {template< class F >struct Output_rep< ::Color, F > {const ::Color& c;static const bool is_specialized = true;Output_rep(const ::Color& c) : c(c){ }std::ostream& operator() (std::ostream& out) const{if (IO::is_ascii(out))out << int(c[0]) << " " << int(c[1]) << " " << int(c[2]) << " " << int(c[3]);elseout.write(reinterpret_cast<const char*>(&c), sizeof(c));return out;}};
} // namespace CGALint main(int argc, char* argv[])
{Point_vector points;// Load point set from a file.const std::string input_file = "cube_point_cloud.ply";std::ifstream input_stream(input_file.c_str());if (input_stream.fail()) {std::cerr << "Failed open file \'" << input_file << "\'" << std::endl;return EXIT_FAILURE;}input_stream.close();std::cout << "Loading point cloud: " << input_file << "...";CGAL::Timer t;t.start();if (!CGAL::IO::read_points(input_file.c_str(), std::back_inserter(points),CGAL::parameters::point_map(Point_map()).normal_map(Normal_map()))) {std::cerr << "Error: cannot read file " << input_file << std::endl;return EXIT_FAILURE;}elsestd::cout << " Done. " << points.size() << " points. Time: "<< t.time() << " sec." << std::endl;// Shape detection.// Default parameter values for the data file cube.pwn.const FT          search_sphere_radius = FT(2) / FT(100);const std::size_t k = 12;const FT          max_distance_to_plane = FT(0.5) / FT(100);const FT          max_accepted_angle = FT(30);const std::size_t min_region_size = 500;Point_map_region_growing point_map_rg(CGAL::make_random_access_property_map(points));Normal_map_region_growing normal_map_rg(CGAL::make_random_access_property_map(points));// Create instances of the classes Neighbor_query and Region_type.Neighbor_query neighbor_query(boost::irange<std::size_t>(0, points.size()), CGAL::parameters::sphere_radius(search_sphere_radius).point_map(point_map_rg));//Neighbor_query neighbor_query(//    boost::irange<std::size_t>(0, points.size()), CGAL::parameters::k_neighbors(k).point_map(point_map_rg));Region_type region_type(CGAL::parameters::maximum_distance(max_distance_to_plane).maximum_angle(max_accepted_angle).minimum_region_size(min_region_size).point_map(point_map_rg).normal_map(normal_map_rg));// Create an instance of the region growing class.Region_growing region_growing(boost::irange<std::size_t>(0, points.size()), neighbor_query, region_type);std::cout << "Extracting planes...";std::vector<typename Region_growing::Primitive_and_region> regions;t.reset();region_growing.detect(std::back_inserter(regions));std::cout << " Done. " << regions.size() << " planes extracted. Time: "<< t.time() << " sec." << std::endl;// Stores the plane index of each point as the third element of the tuple.for (std::size_t i = 0; i < points.size(); ++i)// Uses the get function from the property map that accesses the 3rd element of the tuple.points[i].get<2>() = static_cast<int>(get(region_growing.region_map(), i));// 随机生成颜色(要保证颜色种类大于提取的平面个数)std::vector<Color> rand_colors;for (size_t i = 0; i < regions.size() + 1; i++){Color p_color = {static_cast<unsigned char>(rand() % 256),static_cast<unsigned char>(rand() % 256),static_cast<unsigned char>(rand() % 256), 255 };rand_colors.push_back(p_color);}std::vector<PNC> points_with_color;// 为所属平面相同的点赋予相同的颜色for (std::size_t i = 0; i < points.size(); i++){// 获取单个点坐标Point point = points[i].get<0>();// 获得单个点法线Vector normal = points[i].get<1>();// 获取点对应的颜色，所属平面相同的点颜色相同int plane_index = points[i].get<2>();Color p_color;if (plane_index == -1)    // 未分配平面的点为白色p_color = { 255, 255, 255, 255 };elsep_color = rand_colors[plane_index];points_with_color.push_back(std::make_tuple(point, normal, p_color));}std::ofstream f("result.ply", std::ios::binary);CGAL::IO::set_binary_mode(f); // The PLY file will be written in the binary formatCGAL::IO::write_PLY_with_properties(f, points_with_color,CGAL::make_ply_point_writer(Save_Point_map()),CGAL::make_ply_normal_writer(Save_Normal_map()),std::make_tuple(Save_Color_map(),CGAL::IO::PLY_property<unsigned char>("red"),CGAL::IO::PLY_property<unsigned char>("green"),CGAL::IO::PLY_property<unsigned char>("blue"),CGAL::IO::PLY_property<unsigned char>("alpha")));return EXIT_SUCCESS;
}