VS+ QT（一）使用CGAL实现点云分类

一、依赖的添加

1、附加包含文件目录

$(CGAL)\include
$(CGAL)\auxiliary\gmp\include

$(CGAL)代表CGAL的文件夹地址
例如我的是：($(SolutionDir)为解决方案文件夹目录）
$(SolutionDir)…\dependences\cgal\include;
$(SolutionDir)…\dependences\cgal\auxiliary\gmp\include;

2、附加包含库文件目录

$(CGAL)\auxiliary\gmp\lib;

同上

3、附加包含库文件

libmpfr-4.lib;libgmp-10.lib;

4、头文件

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Classification.h> //分类
#include <CGAL/bounding_box.h> //生成包含点云的矩形（平行于坐标轴）
#include <CGAL/tags.h>

#include <CGAL/IO/read_points.h> //读写ply文件
#include <CGAL/IO/write_ply_points.h>

#include <CGAL/Real_timer.h>//用于计算运行时间

#include <lasreader.hpp>
#include <laswriter.hpp> //LASTools读取文件

#include <cstdlib>
#include <fstream>
#include <iostream>
#include <string>

二、代码

1、类名定义

typedef CGAL::Parallel_if_available_tag Concurrency_tag;

typedef CGAL::Simple_cartesian<double> Kernel;
typedef Kernel::Iso_cuboid_3 Iso_cuboid_3;
typedef Kernel::Point_3 Point;
typedef std::vector<Point> Point_range;
typedef CGAL::Identity_property_map<Point> Pmap; //将 property map 映射到自身的属性映射（通过引用）

namespace Classification = CGAL::Classification;
typedef Classification::Sum_of_weighted_features_classifier Classifier;
typedef Classification::Planimetric_grid<Kernel, Point_range, Pmap>             Planimetric_grid;
typedef Classification::Point_set_neighborhood<Kernel, Point_range, Pmap>       Neighborhood;
typedef Classification::Local_eigen_analysis                                    Local_eigen_analysis;
typedef Classification::Label_handle                                            Label_handle;
typedef Classification::Feature_handle                                          Feature_handle;
typedef Classification::Label_set                                               Label_set;
typedef Classification::Feature_set                                             Feature_set;
typedef Classification::Feature::Distance_to_plane<Point_range, Pmap>           Distance_to_plane;
typedef Classification::Feature::Elevation<Kernel, Point_range, Pmap>           Elevation;
typedef Classification::Feature::Vertical_dispersion<Kernel, Point_range, Pmap> Dispersion;

重新定义类名，减少后续的代码量

2、点云读取

//点云读取
LASreadOpener lasreadopener;

cout << "read File " << Path << endl;
lasreadopener.set_file_name(Path.c_str());
LASreader* lasreader = lasreadopener.open();
//打印头文件信息
// 初始化 header
LASheader lasHeader = lasreader->header;

//声明存储点云XYZ变量
vector<Point> pts;
//转换点云到容器
while (lasreader->read_point())
{
    LASpoint& p = lasreader->point;
    pts.push_back(Point(p.get_x(), p.get_y(), p.get_z()));
}

读取las点云并使用其xyz三个位置坐标赋值给CGAL的点云数组

2、点云预处理

float grid_resolution = 0.34f; //平面网格分辨率
unsigned int number_of_neighbors = 6; //邻域核心数

cout << "Computing " << pts.size() << " Point useful structures" << endl;

Iso_cuboid_3 bbox = CGAL::bounding_box(pts.begin(), pts.end());//计算点的轴对齐边界框

Planimetric_grid grid(pts, Pmap(), bbox, grid_resolution); //基于输入范围构造平面格网
Neighborhood neighborhood(pts, Pmap()); //计算输入点集的空间搜索结构，并允许以索引集的形式访问点的局部邻域。
                                        //它允许用户基于固定范围邻域或固定 K 个相邻域生成模型
// k_neighbor_query 构造 K 邻域查询对象
//预计算并存储使用局部邻域集的点集的所有点的协方差矩阵的特征向量和特征值
Local_eigen_analysis eigen = Local_eigen_analysis::create_from_point_set(pts, Pmap(), neighborhood.k_neighbor_query(number_of_neighbors));

3、添加特征类型计算

    float radius_neighbors = 1.7f;//局部领域的半径范围
    float radius_dtm = 15.0f; //用于数字地形建模的半径(应大于最大建筑物的宽度和长度)

    cout << "Computing features" << endl;
    Feature_set features;

    features.begin_parallel_additions();//初始化结构以并行计算feature
    //具有非平面性水平有助于识别输入的嘈杂部分，例如植被。此功能计算点到局部拟合平面的距离。
    Feature_handle distance_to_plane = features.add<Distance_to_plane>(pts, Pmap(), eigen);//基于与拟合平面的局部距离的feature
    Feature_handle dispersion = features.add<Dispersion>(pts, Pmap(), grid, radius_neighbors);//基于点的局部分散的feature
    Feature_handle elevation = features.add<Elevation>(pts, Pmap(), grid, radius_dtm);//基于本地高度的feature
    features.end_parallel_additions(); //等待并行特征计算结束，然后清除专用数据结构。

4、添加标签名

Label_set labels;
//初始化标签名
Label_handle ground = labels.add("ground");//名
Label_handle vegetation = labels.add("vegetation", CGAL::IO::Color(0, 255, 0));//名，颜色
Label_handle roof = labels.add("roof", CGAL::IO::Color(255, 0, 0), 6);//名，颜色，标准索引（这里是ASPRS建筑索引）

5、初始化分类变量

 Classifier classifier(labels, features);

6、设置权重

//设置权重和影响程度
cout << "Setting weights" << endl;
classifier.set_weight(distance_to_plane, 6.75e-2f);
classifier.set_weight(dispersion, 5.45e-1f);
classifier.set_weight(elevation, 1.47e1f);

cout << "Setting effects" << endl;
classifier.set_effect(ground, distance_to_plane, Classifier::NEUTRAL);
classifier.set_effect(ground, dispersion, Classifier::NEUTRAL);
classifier.set_effect(ground, elevation, Classifier::PENALIZING);

classifier.set_effect(vegetation, distance_to_plane, Classifier::FAVORING);
classifier.set_effect(vegetation, dispersion, Classifier::FAVORING);
classifier.set_effect(vegetation, elevation, Classifier::NEUTRAL);

classifier.set_effect(roof, distance_to_plane, Classifier::NEUTRAL);
classifier.set_effect(roof, dispersion, Classifier::NEUTRAL);
classifier.set_effect(roof, elevation, Classifier::FAVORING);

7、计算并分类

    // 开始计算
    cout << "Classifying" << endl;
    vector<int> label_indices(pts.size(), -1);//标签索引数组

    CGAL::Real_timer t;

    t.start();//计算类型1
    //运行分类算法，而不进行任何正则化。
    //项目之间没有关系，分类能量只是逐项最小化。这种方法很快，但产生的结果是次优的。
    Classification::classify<CGAL::Parallel_if_available_tag>(pts, labels, classifier, label_indices);
    t.stop();
    cout << "Raw classification performed in " << t.time() << " second(s)" << endl;
    t.reset();

    t.start();//计算类型2
    //运行具有局部平滑的分类算法。
    //计算的分类能量在用户定义的项目局部邻域上进行平滑处理。这种方法是效率和更高质量的结果之间的折中。
    Classification::classify_with_local_smoothing<CGAL::Parallel_if_available_tag>
        (pts, Pmap(), labels, classifier, neighborhood.sphere_neighbor_query(radius_neighbors), label_indices);
    t.stop();
    cout << "Classification with local smoothing performed in " << t.time() << " second(s)" << endl;
    t.reset();

    t.start();//计算类型3
    //运行基于图形剪切的全局正则化的分类算法。
    //计算出的分类能量通过 alpha 展开算法进行全局正则化。此方法很慢，但为用户提供了高质量的结果。
    //为了加快计算速度，可以将输入域细分为较小的子集，以便应用几个较小的图形切割而不是一个大的图形切割。
    //这些较小的图形切割的计算可以并行完成。增加子集的数量可以缩短计算时间，但也会降低结果的质量
    Classification::classify_with_graphcut<CGAL::Parallel_if_available_tag>
        (pts, Pmap(), labels, classifier, neighborhood.k_neighbor_query(12), 0.2f, 4, label_indices);
    t.stop();

    cout << "Classification with graphcut performed in " << t.time() << " second(s)" << endl;

三种计算类型任选一种，1到3，精度逐次提高，但计算复杂度也相应提高，导致所需内存也是逐步增加，时间亦是逐步增加

8、标签转换的颜色信息

vector<unsigned char> red, green, blue;
red.reserve(pts.size());
green.reserve(pts.size());
blue.reserve(pts.size());
for (std::size_t i = 0; i < pts.size(); ++i)
{
    Label_handle label = labels[std::size_t(label_indices[i])];
    unsigned r = 0, g = 0, b = 0;
    if (label == ground)
    {
        r = 245; g = 180; b = 0;
    }
    else if (label == vegetation)
    {
        r = 0; g = 255; b = 27;
    }
    else if (label == roof)
    {
        r = 255; g = 0; b = 170;
    }
    red.push_back(r);
    green.push_back(g);
    blue.push_back(b);
}

9、点云保存

    cout << "write File 1.las" << endl;
    //点云写入
    LASwriteOpener lasWriterOpener;
    lasWriterOpener.set_file_name("2.las");
    LASwriter* lasWriter = lasWriterOpener.open(&lasHeader);    // 打开 laswriter

    lasreadopener.set_file_name(Path.c_str());
    lasreader = lasreadopener.open();
    int i = 0;
    // 初始化 point
    while (lasreader->read_point() && i < red.size())
    {
        LASpoint& p = lasreader->point;
        p.set_x(pts[i].x());
        p.set_y(pts[i].y());
        p.set_z(pts[i].z());
        p.set_R(red.at(i));
        p.set_G(green.at(i));
        p.set_B(blue.at(i));
        //写入
        lasWriter->write_point(&p);
        lasWriter->update_inventory(&p);
        ++i;
    }
     关闭
     update the boundary
    lasHeader.set_bounding_box(minX, minY, minZ, maxX, maxY, maxZ);
     update the header
    lasWriter->update_header(&lasHeader, true);

    lasWriter->close();
    delete lasWriter;
    lasWriter = nullptr;

这里只有改变颜色信息，所以直接读取后修改其颜色信息就行