【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]);
            else
                out.write(reinterpret_cast<const char*>(&c), sizeof(c));
            return out;
        }
    };
} // namespace CGAL

int 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;
    }
    else
        std::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 };
        else
            p_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 format
    CGAL::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;
}