18.1 点云处理与分析

18.1.1 点云数据简介

点云数据在现代三维应用中无处不在：

工业测量：激光扫描获取零件几何
建筑测绘：建筑物和遗址的三维数字化
自动驾驶：LiDAR传感器的环境感知
医学影像：CT/MRI数据的三维重建
文化遗产：文物和艺术品的数字化保护

点云处理面临的主要挑战包括：

数据量大：扫描数据通常包含数百万甚至上亿个点
噪声干扰：传感器误差导致点位置存在偏差
不均匀采样：扫描密度随距离和角度变化
缺失数据：遮挡区域导致数据不完整
无拓扑信息：需要重建表面连接关系

CGAL的点云处理框架提供了解决这些挑战的完整工具集。

18.1.2 Point_set_3数据结构

CGAL::Point_set_3是CGAL中专门设计用于存储和处理点云数据的数据结构。与简单的点向量不同，Point_set_3提供了以下核心特性：

动态属性系统：支持添加任意类型的属性（法线、颜色、强度等）
高效内存管理：使用属性容器避免内存碎片
惰性删除机制：标记删除点后延迟回收内存
自动参数绑定：与CGAL算法无缝集成

基本使用

以下代码展示了Point_set_3的基本使用方法：

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Point_set_3.h>
 
#include <iostream>
 
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
typedef Kernel::Vector_3 Vector;
 
typedef CGAL::Point_set_3<Point> Point_set;
 
int main(int, char**)
{
  Point_set point_set;
 
  // 添加点
  point_set.insert(Point(0., 0., 0.));
  point_set.insert(Point(1., 0., 0.));
  point_set.insert(Point(0., 1., 0.));
  point_set.insert(Point(1., 1., 0.));
 
  std::cout << "点云包含 " << point_set.size() << " 个点" << std::endl;
 
  // 遍历点云
  for (Point_set::const_iterator it = point_set.begin();
       it != point_set.end(); ++it)
  {
    const Point& p = point_set.point(*it);
    std::cout << "点 " << *it << ": ("
              << p.x() << ", " << p.y() << ", " << p.z() << ")" << std::endl;
  }
 
  return 0;
}

法线属性管理

法线是点云处理中最重要的属性之一，用于描述表面朝向。Point_set_3提供了便捷的法线管理功能：

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Point_set_3.h>
 
#include <iostream>
 
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
typedef Kernel::Vector_3 Vector;
 
typedef CGAL::Point_set_3<Point> Point_set;
 
int main(int, char**)
{
  Point_set point_set;
 
  // 添加点
  point_set.insert(Point(0., 0., 0.));
  point_set.insert(Point(1., 0., 0.));
  point_set.insert(Point(0., 1., 0.));
 
  // 添加法线属性
  point_set.add_normal_map();
 
  // 设置法线值
  Point_set::iterator it = point_set.begin();
  point_set.normal(*(it++)) = Vector(0., 0., 1.);
  point_set.normal(*(it++)) = Vector(0., 0., 1.);
  point_set.normal(*(it++)) = Vector(0., 0., 1.);
 
  // 同时添加点和法线
  point_set.insert(Point(1., 1., 0.), Vector(0., 0., 1.));
 
  // 遍历点和法线
  for (Point_set::const_iterator it = point_set.begin();
       it != point_set.end(); ++it)
  {
    std::cout << "点: " << point_set.point(*it)
              << " 法线: " << point_set.normal(*it) << std::endl;
  }
 
  return 0;
}

自定义属性

Point_set_3的强大之处在于支持任意类型的自定义属性。以下示例展示了如何添加颜色、强度等属性：

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Point_set_3.h>
 
#include <iostream>
 
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
typedef Kernel::Vector_3 Vector;
 
typedef CGAL::Point_set_3<Point> Point_set;
 
// 颜色结构
struct Color
{
  unsigned char r, g, b;
  Color() : r(0), g(0), b(0) {}
  Color(unsigned char _r, unsigned char _g, unsigned char _b)
    : r(_r), g(_g), b(_b) {}
};
 
int main(int, char**)
{
  Point_set point_set;
 
  // 添加自定义属性
  auto [color_map, color_added] = point_set.add_property_map<Color>("color", Color(255, 255, 255));
  auto [intensity_map, intensity_added] = point_set.add_property_map<float>("intensity", 0.0f);
  auto [label_map, label_added] = point_set.add_property_map<int>("label", -1);
 
  // 添加带属性的点
  for (int i = 0; i < 10; ++i)
  {
    Point_set::iterator it = point_set.insert(Point(i * 0.1, 0., 0.));
    color_map[*it] = Color(i * 25, 255 - i * 25, 128);
    intensity_map[*it] = static_cast<float>(i) / 10.0f;
    label_map[*it] = i % 3;  // 0, 1, 2 三类标签
  }
 
  // 查询属性
  std::cout << "属性列表:" << std::endl;
  for (const std::string& prop : point_set.properties())
    std::cout << " - " << prop << std::endl;
 
  // 访问属性
  for (Point_set::const_iterator it = point_set.begin();
       it != point_set.end(); ++it)
  {
    const Color& c = color_map[*it];
    std::cout << "点 " << point_set.point(*it)
              << " RGB(" << (int)c.r << "," << (int)c.g << "," << (int)c.b << ")"
              << " 强度=" << intensity_map[*it]
              << " 标签=" << label_map[*it] << std::endl;
  }
 
  // 删除属性
  point_set.remove_property_map(intensity_map);
 
  return 0;
}

空间索引集成

Point_set_3内部使用KD树实现高效的空间查询。当调用点云处理算法时，CGAL会自动构建和使用这些空间索引结构：

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Point_set_3.h>
#include <CGAL/Point_set_3/IO/XYZ.h>
 
#include <iostream>
#include <fstream>
 
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
 
typedef CGAL::Point_set_3<Point> Point_set;
 
int main(int argc, char** argv)
{
  Point_set point_set;
 
  // 从文件读取点云
  std::ifstream in((argc > 1) ? argv[1] : "data/points.xyz");
  if (!in)
  {
    std::cerr << "无法打开输入文件" << std::endl;
    return 1;
  }
 
  in >> point_set;
  std::cout << "读取了 " << point_set.size() << " 个点" << std::endl;
 
  // 显示点云信息
  std::cout << point_set.info();
 
  // 写入文件
  std::ofstream out("output.xyz");
  out << point_set;
 
  return 0;
}

点集合并

Point_set_3支持高效的点集合并操作：

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Point_set_3.h>
 
#include <iostream>
 
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
 
typedef CGAL::Point_set_3<Point> Point_set;
 
int main(int, char**)
{
  Point_set point_set_1, point_set_2;
 
  // 填充第一个点云
  for (int i = 0; i < 5; ++i)
    point_set_1.insert(Point(i * 1.0, 0., 0.));
 
  // 填充第二个点云
  for (int i = 0; i < 5; ++i)
    point_set_2.insert(Point(i * 1.0, 1., 0.));
 
  std::cout << "合并前: 点云1=" << point_set_1.size()
            << " 点云2=" << point_set_2.size() << std::endl;
 
  // 合并点云
  point_set_1.join(point_set_2);
 
  std::cout << "合并后: 点云1=" << point_set_1.size()
            << " 点云2=" << point_set_2.size() << std::endl;
 
  // 或者使用运算符+=
  Point_set point_set_3;
  point_set_3.insert(Point(0., 2., 0.));
  point_set_1 += point_set_3;
 
  std::cout << "最终点云大小: " << point_set_1.size() << std::endl;
 
  return 0;
}

18.1.3 法线估计

法线估计是点云处理的基础步骤，它为每个点计算表面法线方向。CGAL提供了多种法线估计算法。

PCA局部平面拟合

主成分分析（PCA）是最常用的法线估计方法。它通过拟合局部邻域的最小二乘平面来计算法线：

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Point_set_3.h>
#include <CGAL/pca_estimate_normals.h>
#include <CGAL/mst_orient_normals.h>
 
#include <iostream>
 
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
typedef Kernel::Vector_3 Vector;
 
typedef CGAL::Point_set_3<Point> Point_set;
 
int main(int, char**)
{
  Point_set point_set;
 
  // 生成示例点云（平面）
  for (int i = 0; i < 100; ++i)
    for (int j = 0; j < 100; ++j)
      point_set.insert(Point(i * 0.01, j * 0.01, 0.0));
 
  std::cout << "原始点云: " << point_set.size() << " 个点" << std::endl;
 
  // 添加法线属性
  point_set.add_normal_map();
 
  // 使用PCA估计法线
  // 参数：
  // - Sequential_tag: 串行执行
  // - 12: 每个点使用12个最近邻
  CGAL::pca_estimate_normals<CGAL::Sequential_tag>
    (point_set,
     12,  // 邻域大小
     point_set.parameters().
     neighbor_radius(0.0));  // 0表示不使用半径限制
 
  std::cout << "法线估计完成" << std::endl;
 
  // 显示部分法线
  int count = 0;
  for (Point_set::const_iterator it = point_set.begin();
       it != point_set.end() && count < 5; ++it, ++count)
  {
    std::cout << "点: " << point_set.point(*it)
              << " 法线: " << point_set.normal(*it) << std::endl;
  }
 
  return 0;
}

法线方向一致性

PCA估计的法线方向是随机的（可能指向表面内部或外部）。为了获得一致朝向的法线，需要使用方向传播算法：

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Point_set_3.h>
#include <CGAL/pca_estimate_normals.h>
#include <CGAL/mst_orient_normals.h>
 
#include <iostream>
 
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
typedef Kernel::Vector_3 Vector;
 
typedef CGAL::Point_set_3<Point> Point_set;
 
int main(int, char**)
{
  Point_set point_set;
 
  // 生成球面上的点
  const double radius = 1.0;
  const int n_theta = 50;
  const int n_phi = 50;
 
  for (int i = 0; i < n_theta; ++i)
  {
    double theta = 2.0 * CGAL_PI * i / n_theta;
    for (int j = 0; j < n_phi; ++j)
    {
      double phi = CGAL_PI * j / n_phi;
      double x = radius * sin(phi) * cos(theta);
      double y = radius * sin(phi) * sin(theta);
      double z = radius * cos(phi);
      point_set.insert(Point(x, y, z));
    }
  }
 
  std::cout << "生成球面点云: " << point_set.size() << " 个点" << std::endl;
 
  // 添加法线属性
  point_set.add_normal_map();
 
  // 第一步：估计法线（方向随机）
  CGAL::pca_estimate_normals<CGAL::Sequential_tag>
    (point_set, 12, point_set.parameters());
 
  std::cout << "PCA法线估计完成" << std::endl;
 
  // 第二步：使用最小生成树（MST）定向法线
  // 算法从Z坐标最高的点开始，向外传播法线方向
  Point_set::iterator unoriented_begin =
    CGAL::mst_orient_normals(point_set,
                              12,  // 邻域大小
                              point_set.parameters());
 
  // 移除未能定向的点（如果有）
  point_set.remove(unoriented_begin, point_set.end());
  point_set.collect_garbage();
 
  std::cout << "法线定向完成，剩余 " << point_set.size() << " 个点" << std::endl;
 
  // 验证法线方向（球面法线应指向外侧）
  int count = 0;
  for (Point_set::const_iterator it = point_set.begin();
       it != point_set.end() && count < 5; ++it, ++count)
  {
    Point p = point_set.point(*it);
    Vector n = point_set.normal(*it);
    // 对于球面，法线应与位置向量同向
    Vector expected = Vector(p.x(), p.y(), p.z()) / radius;
    double dot = n * expected;
    std::cout << "点: (" << p.x() << "," << p.y() << "," << p.z() << ")"
              << " 法线: (" << n.x() << "," << n.y() << "," << n.z() << ")"
              << " 一致性: " << (dot > 0 ? "正确" : "错误") << std::endl;
  }
 
  return 0;
}

不同邻域大小的影响

邻域大小的选择对法线估计质量有重要影响：

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Point_set_3.h>
#include <CGAL/pca_estimate_normals.h>
 
#include <iostream>
#include <cmath>
 
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
typedef Kernel::Vector_3 Vector;
 
typedef CGAL::Point_set_3<Point> Point_set;
 
// 计算法线误差（与真实法线的角度差）
double normal_error(const Vector& estimated, const Vector& ground_truth)
{
  double dot = estimated * ground_truth;
  // 处理方向不确定性
  if (dot < 0) dot = -dot;
  return std::acos(std::min(1.0, std::max(-1.0, dot))) * 180.0 / CGAL_PI;
}
 
int main(int, char**)
{
  // 创建带噪声的平面点云
  Point_set point_set;
  const Vector ground_truth_normal(0., 0., 1.);
 
  CGAL::Random random;
  for (int i = 0; i < 1000; ++i)
  {
    double x = random.get_double(0., 1.);
    double y = random.get_double(0., 1.);
    double z = random.get_double(-0.01, 0.01);  // 添加噪声
    point_set.insert(Point(x, y, z));
  }
 
  // 测试不同邻域大小
  std::vector<unsigned int> k_values = {6, 12, 24, 48};
 
  for (unsigned int k : k_values)
  {
    Point_set test_set = point_set;
    test_set.add_normal_map();
 
    CGAL::pca_estimate_normals<CGAL::Sequential_tag>
      (test_set, k, test_set.parameters());
 
    // 计算平均误差
    double total_error = 0.0;
    for (Point_set::const_iterator it = test_set.begin();
         it != test_set.end(); ++it)
    {
      total_error += normal_error(test_set.normal(*it), ground_truth_normal);
    }
    double avg_error = total_error / test_set.size();
 
    std::cout << "k=" << k << ": 平均法线误差 = " << avg_error << "度" << std::endl;
  }
 
  return 0;
}

邻域大小选择指南：

邻域大小	适用场景	优缺点
6-8	高分辨率扫描、锐利特征	细节保留好，但对噪声敏感
12-18	一般用途	平衡噪声鲁棒性和细节保留
24-36	噪声较大的数据	平滑效果好，可能丢失细节
48+	严重噪声数据	高度平滑，特征模糊

Jet拟合估计法线

对于曲面特征明显的点云，可以使用基于Jet拟合的法线估计：

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Point_set_3.h>
#include <CGAL/jet_estimate_normals.h>
 
#include <iostream>
 
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
typedef Kernel::Vector_3 Vector;
 
typedef CGAL::Point_set_3<Point> Point_set;
 
int main(int, char**)
{
  Point_set point_set;
 
  // 生成带曲面的点云
  for (int i = 0; i < 100; ++i)
  {
    double x = i / 100.0;
    for (int j = 0; j < 100; ++j)
    {
      double y = j / 100.0;
      double z = 0.1 * sin(2 * CGAL_PI * x) * sin(2 * CGAL_PI * y);
      point_set.insert(Point(x, y, z));
    }
  }
 
  point_set.add_normal_map();
 
  // 使用Jet拟合估计法线
  // degree_fitting: 拟合多项式的阶数（默认2）
  // degree_monge: Monge形式的阶数（默认2）
  CGAL::jet_estimate_normals<CGAL::Sequential_tag>
    (point_set,
     18,  // 邻域大小
     point_set.parameters().
     degree_fitting(2).
     degree_monge(2));
 
  std::cout << "Jet法线估计完成，共 " << point_set.size() << " 个点" << std::endl;
 
  return 0;
}

18.1.4 移动最小二乘平滑

移动最小二乘（MLS）是一种强大的点云平滑和重采样技术。它通过将点投影到局部拟合的曲面上来实现平滑。

曲面重建原理

MLS方法的核心思想是：对于每个点，在其邻域内拟合一个隐式曲面，然后将该点投影到曲面上。这种方法可以有效去除噪声同时保留几何特征。

投影与平滑

CGAL提供了基于Jet拟合的MLS平滑实现：

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Point_set_3.h>
#include <CGAL/jet_smooth_point_set.h>
#include <CGAL/pca_estimate_normals.h>
 
#include <iostream>
 
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
typedef Kernel::Vector_3 Vector;
 
typedef CGAL::Point_set_3<Point> Point_set;
 
int main(int, char**)
{
  Point_set point_set;
  CGAL::Random random;
 
  // 生成带噪声的曲面点云
  for (int i = 0; i < 100; ++i)
  {
    double x = i / 100.0;
    for (int j = 0; j < 100; ++j)
    {
      double y = j / 100.0;
      // 原始曲面：z = sin(2*pi*x) * sin(2*pi*y)
      double z = sin(2 * CGAL_PI * x) * sin(2 * CGAL_PI * y);
      // 添加噪声
      z += random.get_double(-0.05, 0.05);
      point_set.insert(Point(x, y, z));
    }
  }
 
  std::cout << "原始点云: " << point_set.size() << " 个点（带噪声）" << std::endl;
 
  // 估计法线（平滑前）
  point_set.add_normal_map();
  CGAL::pca_estimate_normals<CGAL::Sequential_tag>
    (point_set, 12, point_set.parameters());
 
  // MLS平滑
  // 参数：
  // - k: 邻域大小
  // - degree_fitting: 拟合多项式阶数
  // - degree_monge: Monge形式阶数
  CGAL::jet_smooth_point_set<CGAL::Sequential_tag>
    (point_set,
     24,  // 邻域大小（平滑需要较大邻域）
     point_set.parameters().
     degree_fitting(2).
     degree_monge(2));
 
  std::cout << "MLS平滑完成" << std::endl;
 
  // 重新估计法线（平滑后）
  CGAL::pca_estimate_normals<CGAL::Sequential_tag>
    (point_set, 12, point_set.parameters());
 
  return 0;
}

参数选择指南

MLS平滑的效果高度依赖于参数选择：

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Point_set_3.h>
#include <CGAL/jet_smooth_point_set.h>
 
#include <iostream>
 
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
 
typedef CGAL::Point_set_3<Point> Point_set;
 
// 计算点云到原始曲面的平均距离（用于评估平滑质量）
double compute_error(const Point_set& point_set)
{
  double total_error = 0.0;
  for (Point_set::const_iterator it = point_set.begin();
       it != point_set.end(); ++it)
  {
    Point p = point_set.point(*it);
    // 原始曲面: z = sin(2*pi*x) * sin(2*pi*y)
    double z_true = sin(2 * CGAL_PI * p.x()) * sin(2 * CGAL_PI * p.y());
    double error = std::abs(p.z() - z_true);
    total_error += error;
  }
  return total_error / point_set.size();
}
 
int main(int, char**)
{
  CGAL::Random random;
 
  // 测试不同参数组合
  std::vector<unsigned int> k_values = {12, 18, 24};
  std::vector<unsigned int> degree_values = {2, 3, 4};
 
  for (unsigned int k : k_values)
  {
    for (unsigned int degree : degree_values)
    {
      // 生成带噪声的点云
      Point_set point_set;
      for (int i = 0; i < 100; ++i)
      {
        double x = i / 100.0;
        for (int j = 0; j < 100; ++j)
        {
          double y = j / 100.0;
          double z = sin(2 * CGAL_PI * x) * sin(2 * CGAL_PI * y);
          z += random.get_double(-0.05, 0.05);
          point_set.insert(Point(x, y, z));
        }
      }
 
      // 应用MLS平滑
      CGAL::jet_smooth_point_set<CGAL::Sequential_tag>
        (point_set, k,
         point_set.parameters().
         degree_fitting(degree).
         degree_monge(degree));
 
      double error = compute_error(point_set);
      std::cout << "k=" << k << ", degree=" << degree
                << ": 平均误差 = " << error << std::endl;
    }
  }
 
  return 0;
}

MLS参数选择建议：

参数	说明	推荐值
k	邻域大小	数据点数的平方根或18-30
degree_fitting	拟合多项式阶数	2（一般），3-4（复杂曲面）
degree_monge	Monge形式阶数	与degree_fitting相同

注意事项：

邻域过小：无法充分平滑噪声
邻域过大：过度平滑，丢失细节
阶数过高：计算开销大，可能过拟合
阶数过低：无法捕捉复杂曲面特征

18.1.5 RANSAC形状检测

RANSAC（随机采样一致性）是检测点云中几何形状的有效方法。CGAL实现了高效的RANSAC算法，可以检测平面、球体、圆柱、圆锥和圆环等形状。

平面、圆柱、球体检测

以下示例展示了如何使用RANSAC检测点云中的基本几何形状：

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Point_set_3.h>
#include <CGAL/Shape_detection/Efficient_RANSAC.h>
#include <CGAL/pca_estimate_normals.h>
 
#include <iostream>
 
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
typedef Kernel::Vector_3 Vector;
 
typedef CGAL::Point_set_3<Point> Point_set;
 
// RANSAC类型定义
typedef CGAL::Shape_detection::Efficient_RANSAC_traits
  <Kernel, Point_set, Point_set::Point_map, Point_set::Vector_map> Traits;
typedef CGAL::Shape_detection::Efficient_RANSAC<Traits>          Efficient_ransac;
 
// 形状类型
typedef CGAL::Shape_detection::Plane<Traits>                     Plane;
typedef CGAL::Shape_detection::Sphere<Traits>                    Sphere;
typedef CGAL::Shape_detection::Cylinder<Traits>                  Cylinder;
typedef CGAL::Shape_detection::Cone<Traits>                      Cone;
typedef CGAL::Shape_detection::Torus<Traits>                     Torus;
 
int main(int, char**)
{
  Point_set point_set;
 
  // 生成测试数据：平面 + 球体
  CGAL::Random random;
 
  // 添加平面点（带少量噪声）
  for (int i = 0; i < 500; ++i)
  {
    double x = random.get_double(-1., 1.);
    double y = random.get_double(-1., 1.);
    double z = random.get_double(-0.01, 0.01);  // 接近z=0平面
    point_set.insert(Point(x, y, z));
  }
 
  // 添加球面点
  double radius = 0.5;
  Point center(0., 0., 1.);
  for (int i = 0; i < 500; ++i)
  {
    // 球坐标随机采样
    double theta = random.get_double(0., 2. * CGAL_PI);
    double phi = random.get_double(0., CGAL_PI);
    double x = center.x() + radius * sin(phi) * cos(theta);
    double y = center.y() + radius * sin(phi) * sin(theta);
    double z = center.z() + radius * cos(phi);
    point_set.insert(Point(x, y, z));
  }
 
  std::cout << "输入点云: " << point_set.size() << " 个点" << std::endl;
 
  // 估计法线（RANSAC需要法线信息）
  point_set.add_normal_map();
  CGAL::pca_estimate_normals<CGAL::Sequential_tag>
    (point_set, 12, point_set.parameters());
 
  // 创建RANSAC检测器
  Efficient_ransac ransac;
  ransac.set_input(point_set,
                   point_set.point_map(),
                   point_set.normal_map());
 
  // 注册要检测的形状类型
  ransac.add_shape_factory<Plane>();
  ransac.add_shape_factory<Sphere>();
  ransac.add_shape_factory<Cylinder>();
 
  // 设置检测参数
  Efficient_ransac::Parameters parameters;
  parameters.probability = 0.05;        // 检测概率（越低越严格）
  parameters.min_points = 100;          // 最小点数阈值
  parameters.epsilon = 0.02;            // 点到形状的最大距离
  parameters.cluster_epsilon = 0.1;     // 聚类距离阈值
  parameters.normal_threshold = 0.9;    // 法线一致性阈值（cos角度）
 
  // 执行检测
  ransac.detect(parameters);
 
  // 输出检测结果
  std::cout << "\n检测到 " << ransac.shapes().size() << " 个形状:" << std::endl;
 
  int shape_index = 0;
  for (const auto& shape : ransac.shapes())
  {
    std::cout << "\n形状 " << (++shape_index) << ":";
 
    // 获取形状信息
    if (Plane* plane = dynamic_cast<Plane*>(shape.get()))
    {
      std::cout << " 类型: 平面" << std::endl;
      std::cout << "  平面方程: " << plane->plane() << std::endl;
      std::cout << "  点数: " << shape->indices_of_assigned_points().size() << std::endl;
    }
    else if (Sphere* sphere = dynamic_cast<Sphere*>(shape.get()))
    {
      std::cout << " 类型: 球体" << std::endl;
      std::cout << "  中心: " << sphere->center() << std::endl;
      std::cout << "  半径: " << sphere->radius() << std::endl;
      std::cout << "  点数: " << shape->indices_of_assigned_points().size() << std::endl;
    }
    else if (Cylinder* cylinder = dynamic_cast<Cylinder*>(shape.get()))
    {
      std::cout << " 类型: 圆柱" << std::endl;
      std::cout << "  轴点: " << cylinder->axis_point() << std::endl;
      std::cout << "  轴方向: " << cylinder->axis_direction() << std::endl;
      std::cout << "  半径: " << cylinder->radius() << std::endl;
      std::cout << "  点数: " << shape->indices_of_assigned_points().size() << std::endl;
    }
  }
 
  std::cout << "\n未分配点数: " << ransac.number_of_unassigned_points() << std::endl;
 
  return 0;
}

参数设置策略

RANSAC算法的性能高度依赖于参数选择。以下是参数调优的详细指南：

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Point_set_3.h>
#include <CGAL/Shape_detection/Efficient_RANSAC.h>
#include <CGAL/pca_estimate_normals.h>
#include <CGAL/compute_average_spacing.h>
 
#include <iostream>
 
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
typedef Kernel::Vector_3 Vector;
 
typedef CGAL::Point_set_3<Point> Point_set;
typedef CGAL::Shape_detection::Efficient_RANSAC_traits
  <Kernel, Point_set, Point_set::Point_map, Point_set::Vector_map> Traits;
typedef CGAL::Shape_detection::Efficient_RANSAC<Traits> Efficient_ransac;
typedef CGAL::Shape_detection::Plane<Traits> Plane;
 
// 自动计算RANSAC参数
Efficient_ransac::Parameters compute_auto_parameters(
  const Point_set& point_set,
  double noise_level = 0.01)  // 噪声水平估计（相对于包围盒对角线）
{
  Efficient_ransac::Parameters params;
 
  // 计算平均点间距
  double avg_spacing = CGAL::compute_average_spacing<CGAL::Sequential_tag>
    (point_set, 6, point_set.parameters());
 
  std::cout << "平均点间距: " << avg_spacing << std::endl;
 
  // 根据点间距设置参数
  params.epsilon = avg_spacing * 2.0;        // 点到形状距离
  params.cluster_epsilon = avg_spacing * 5.0; // 聚类距离
  params.normal_threshold = 0.85;             // 法线阈值（约30度）
  params.min_points = point_set.size() / 50;  // 至少2%的点
  params.probability = 0.01;                  // 高置信度
 
  return params;
}
 
int main(int, char**)
{
  Point_set point_set;
  CGAL::Random random;
 
  // 生成包含多个平面的点云
  // 平面1: z = 0
  for (int i = 0; i < 300; ++i)
    point_set.insert(Point(random.get_double(-1., 1.),
                           random.get_double(-1., 1.),
                           random.get_double(-0.01, 0.01)));
 
  // 平面2: x = 0
  for (int i = 0; i < 300; ++i)
    point_set.insert(Point(random.get_double(-0.01, 0.01),
                           random.get_double(-1., 1.),
                           random.get_double(-1., 1.)));
 
  // 平面3: y = 0
  for (int i = 0; i < 300; ++i)
    point_set.insert(Point(random.get_double(-1., 1.),
                           random.get_double(-0.01, 0.01),
                           random.get_double(-1., 1.)));
 
  point_set.add_normal_map();
  CGAL::pca_estimate_normals<CGAL::Sequential_tag>
    (point_set, 12, point_set.parameters());
 
  // 使用自动参数
  std::cout << "=== 自动参数 ===" << std::endl;
  auto params = compute_auto_parameters(point_set);
 
  Efficient_ransac ransac;
  ransac.set_input(point_set,
                   point_set.point_map(),
                   point_set.normal_map());
  ransac.add_shape_factory<Plane>();
 
  ransac.detect(params);
 
  std::cout << "检测到 " << ransac.shapes().size() << " 个平面" << std::endl;
 
  return 0;
}

RANSAC参数详解：

参数	说明	设置建议
probability	检测失败概率上限	0.01-0.05（值越小越严格）
min_points	形状最小点数	根据数据规模，通常为1-5%
epsilon	点到形状最大距离	2-3倍平均点间距
cluster_epsilon	空间聚类距离	5-10倍平均点间距
normal_threshold	法线一致性阈值	0.8-0.9（cos值，对应25-36度）

结果过滤与后处理

检测完成后，通常需要对结果进行过滤和优化：

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Point_set_3.h>
#include <CGAL/Shape_detection/Efficient_RANSAC.h>
#include <CGAL/pca_estimate_normals.h>
#include <CGAL/structure_point_set.h>
 
#include <iostream>
#include <vector>
 
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
typedef Kernel::Vector_3 Vector;
typedef Kernel::Plane_3 Plane_3;
 
typedef CGAL::Point_set_3<Point> Point_set;
typedef CGAL::Shape_detection::Efficient_RANSAC_traits
  <Kernel, Point_set, Point_set::Point_map, Point_set::Vector_map> Traits;
typedef CGAL::Shape_detection::Efficient_RANSAC<Traits> Efficient_ransac;
typedef CGAL::Shape_detection::Plane<Traits> Plane_shape;
 
// 过滤和优化检测结果
void filter_and_refine_shapes(
  Efficient_ransac& ransac,
  double min_coverage = 0.05,      // 最小覆盖率
  double max_overlap = 0.8)        // 最大重叠率
{
  auto shapes = ransac.shapes();
  std::vector<std::shared_ptr<Efficient_ransac::Shape>> filtered_shapes;
 
  std::size_t total_points = 0;
  for (const auto& shape : shapes)
    total_points += shape->indices_of_assigned_points().size();
 
  std::cout << "原始形状数: " << shapes.size() << std::endl;
 
  for (const auto& shape : shapes)
  {
    std::size_t shape_points = shape->indices_of_assigned_points().size();
    double coverage = double(shape_points) / total_points;
 
    // 过滤小形状
    if (coverage < min_coverage)
    {
      std::cout << "过滤小形状（覆盖率: " << coverage * 100 << "%）" << std::endl;
      continue;
    }
 
    filtered_shapes.push_back(shape);
  }
 
  std::cout << "过滤后形状数: " << filtered_shapes.size() << std::endl;
}
 
// 提取平面并构建结构化点集
void extract_structured_point_set(
  const Point_set& input_points,
  Efficient_ransac& ransac)
{
  // 提取检测到的平面
  std::vector<Plane_3> planes;
  std::vector<int> point_plane_indices(input_points.size(), -1);
 
  int plane_idx = 0;
  for (const auto& shape : ransac.shapes())
  {
    if (Plane_shape* plane = dynamic_cast<Plane_shape*>(shape.get()))
    {
      planes.push_back(plane->plane());
 
      // 标记属于该平面的点
      for (std::size_t idx : shape->indices_of_assigned_points())
        point_plane_indices[idx] = plane_idx;
 
      ++plane_idx;
    }
  }
 
  std::cout << "提取了 " << planes.size() << " 个平面" << std::endl;
 
  // 创建平面索引属性映射
  auto [plane_index_map, added] = const_cast<Point_set&>(input_points)
    .add_property_map<int>("plane_index", -1);
 
  std::size_t idx = 0;
  for (Point_set::const_iterator it = input_points.begin();
       it != input_points.end(); ++it, ++idx)
  {
    plane_index_map[*it] = point_plane_indices[idx];
  }
}
 
int main(int, char**)
{
  Point_set point_set;
  CGAL::Random random;
 
  // 生成建筑结构点云（多个相交平面）
  // 地面
  for (int i = 0; i < 400; ++i)
    point_set.insert(Point(random.get_double(-2., 2.),
                           random.get_double(-2., 2.),
                           random.get_double(-0.01, 0.01)));
 
  // 墙面1
  for (int i = 0; i < 300; ++i)
    point_set.insert(Point(random.get_double(-0.01, 0.01),
                           random.get_double(-2., 2.),
                           random.get_double(0., 2.)));
 
  // 墙面2
  for (int i = 0; i < 300; ++i)
    point_set.insert(Point(random.get_double(-2., 2.),
                           random.get_double(-0.01, 0.01),
                           random.get_double(0., 2.)));
 
  point_set.add_normal_map();
  CGAL::pca_estimate_normals<CGAL::Sequential_tag>
    (point_set, 12, point_set.parameters());
 
  // 检测形状
  Efficient_ransac ransac;
  ransac.set_input(point_set,
                   point_set.point_map(),
                   point_set.normal_map());
  ransac.add_shape_factory<Plane_shape>();
 
  Efficient_ransac::Parameters params;
  params.probability = 0.01;
  params.min_points = 100;
  params.epsilon = 0.05;
  params.cluster_epsilon = 0.2;
  params.normal_threshold = 0.9;
 
  ransac.detect(params);
 
  // 过滤结果
  filter_and_refine_shapes(ransac);
 
  // 提取结构化信息
  extract_structured_point_set(point_set, ransac);
 
  return 0;
}

18.1.6 点云分类基础

点云分类是为每个点分配语义标签的过程，如地面、植被、建筑等。CGAL提供了特征提取和简单分类的工具。

特征提取

特征提取是分类的基础。常用的点云特征包括：

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Point_set_3.h>
#include <CGAL/pca_estimate_normals.h>
#include <CGAL/compute_average_spacing.h>
 
#include <iostream>
#include <vector>
#include <cmath>
 
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
typedef Kernel::Vector_3 Vector;
typedef Kernel::FT FT;
 
typedef CGAL::Point_set_3<Point> Point_set;
 
// 点特征结构
struct PointFeatures
{
  double linearity;      // 线性度（1D特征）
  double planarity;      // 平面度（2D特征）
  double scattering;     // 散射度（3D特征）
  double verticality;    // 垂直度
  double height;         // 高度
  double density;        // 局部密度
};
 
// 基于PCA的特征提取
PointFeatures compute_pca_features(
  const std::vector<Point>& neighbors,
  const Point& query_point)
{
  PointFeatures features = {0, 0, 0, 0, 0, 0};
 
  if (neighbors.size() < 3)
    return features;
 
  // 计算质心
  Point centroid(0, 0, 0);
  for (const Point& p : neighbors)
    centroid = Point(centroid.x() + p.x(),
                     centroid.y() + p.y(),
                     centroid.z() + p.z());
  centroid = Point(centroid.x() / neighbors.size(),
                   centroid.y() / neighbors.size(),
                   centroid.z() / neighbors.size());
 
  // 构建协方差矩阵
  double cov[3][3] = {{0}};
  for (const Point& p : neighbors)
  {
    double dx = p.x() - centroid.x();
    double dy = p.y() - centroid.y();
    double dz = p.z() - centroid.z();
 
    cov[0][0] += dx * dx;
    cov[0][1] += dx * dy;
    cov[0][2] += dx * dz;
    cov[1][1] += dy * dy;
    cov[1][2] += dy * dz;
    cov[2][2] += dz * dz;
  }
 
  for (int i = 0; i < 3; ++i)
    for (int j = 0; j < 3; ++j)
      cov[i][j] /= neighbors.size();
 
  cov[1][0] = cov[0][1];
  cov[2][0] = cov[0][2];
  cov[2][1] = cov[1][2];
 
  // 计算特征值（简化计算，实际应使用特征值分解）
  // 这里使用迹和行列式近似
  double trace = cov[0][0] + cov[1][1] + cov[2][2];
 
  // 简化的特征值估计
  double lambda1 = std::max({cov[0][0], cov[1][1], cov[2][2]});
  double lambda3 = std::min({cov[0][0], cov[1][1], cov[2][2]});
  double lambda2 = trace - lambda1 - lambda3;
 
  // 归一化
  if (trace > 0)
  {
    lambda1 /= trace;
    lambda2 /= trace;
    lambda3 /= trace;
  }
 
  // 计算几何特征
  features.linearity = (lambda1 - lambda2) / lambda1;
  features.planarity = (lambda2 - lambda3) / lambda1;
  features.scattering = lambda3 / lambda1;
 
  // 垂直度（基于主方向）
  Vector vertical(0, 0, 1);
  // 简化：使用z分量近似
  features.verticality = std::abs(centroid.z() - query_point.z());
 
  // 高度
  features.height = query_point.z();
 
  // 密度
  features.density = neighbors.size();
 
  return features;
}
 
int main(int, char**)
{
  Point_set point_set;
  CGAL::Random random;
 
  // 生成混合点云
  // 地面（平面）
  for (int i = 0; i < 200; ++i)
    point_set.insert(Point(random.get_double(-1., 1.),
                           random.get_double(-1., 1.),
                           random.get_double(-0.01, 0.01)));
 
  // 垂直结构（线性）
  for (int i = 0; i < 100; ++i)
    point_set.insert(Point(random.get_double(-0.02, 0.02),
                           random.get_double(-0.02, 0.02),
                           random.get_double(0., 1.)));
 
  // 散乱点
  for (int i = 0; i < 100; ++i)
    point_set.insert(Point(random.get_double(-1., 1.),
                           random.get_double(-1., 1.),
                           random.get_double(0.5, 1.5)));
 
  std::cout << "点云大小: " << point_set.size() << std::endl;
 
  // 为每个点计算特征
  std::vector<PointFeatures> features;
 
  // 这里简化处理，实际应使用KD树查找邻域
  std::cout << "特征提取完成" << std::endl;
 
  return 0;
}

简单分类器实现

基于提取的特征，可以实现简单的规则分类器：

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Point_set_3.h>
 
#include <iostream>
#include <vector>
 
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
 
typedef CGAL::Point_set_3<Point> Point_set;
 
// 语义标签
enum class SemanticLabel
{
  UNKNOWN = 0,
  GROUND = 1,      // 地面
  VEGETATION = 2,  // 植被
  BUILDING = 3,    // 建筑
  POLE = 4         // 杆状物
};
 
// 简单规则分类器
class RuleBasedClassifier
{
public:
  struct Features
  {
    double planarity;
    double linearity;
    double scattering;
    double height;
    double verticality;
  };
 
  SemanticLabel classify(const Features& f) const
  {
    // 地面：低高度 + 高平面度
    if (f.height < 0.2 && f.planarity > 0.7)
      return SemanticLabel::GROUND;
 
    // 杆状物：高线性度 + 高垂直度
    if (f.linearity > 0.8 && f.verticality > 0.9)
      return SemanticLabel::POLE;
 
    // 建筑：中等高度 + 高平面度
    if (f.height > 0.5 && f.planarity > 0.6)
      return SemanticLabel::BUILDING;
 
    // 植被：高散射度
    if (f.scattering > 0.3)
      return SemanticLabel::VEGETATION;
 
    return SemanticLabel::UNKNOWN;
  }
};
 
// 应用分类
void classify_point_set(Point_set& point_set,
                        const std::vector<SemanticLabel>& labels)
{
  // 添加标签属性
  auto [label_map, added] = point_set.add_property_map<int>("label", 0);
 
  std::size_t idx = 0;
  for (Point_set::iterator it = point_set.begin();
       it != point_set.end(); ++it, ++idx)
  {
    label_map[*it] = static_cast<int>(labels[idx]);
  }
 
  // 统计各类别数量
  std::vector<int> counts(5, 0);
  for (SemanticLabel label : labels)
    counts[static_cast<int>(label)]++;
 
  std::cout << "分类统计:" << std::endl;
  std::cout << "  地面: " << counts[1] << std::endl;
  std::cout << "  植被: " << counts[2] << std::endl;
  std::cout << "  建筑: " << counts[3] << std::endl;
  std::cout << "  杆状: " << counts[4] << std::endl;
  std::cout << "  未知: " << counts[0] << std::endl;
}
 
int main(int, char**)
{
  // 示例：使用规则分类器
  RuleBasedClassifier classifier;
 
  // 测试特征
  std::vector<RuleBasedClassifier::Features> test_features = {
    {0.9, 0.1, 0.0, 0.05, 0.1},   // 地面
    {0.1, 0.1, 0.8, 1.2, 0.5},    // 植被
    {0.8, 0.2, 0.0, 2.5, 0.9},    // 建筑墙面
    {0.1, 0.9, 0.0, 1.5, 0.95},   // 杆状
  };
 
  for (const auto& f : test_features)
  {
    SemanticLabel label = classifier.classify(f);
    std::cout << "特征 (平面度=" << f.planarity
              << ", 线性度=" << f.linearity
              << ", 高度=" << f.height << ") -> ";
 
    switch (label)
    {
      case SemanticLabel::GROUND:
        std::cout << "地面" << std::endl;
        break;
      case SemanticLabel::VEGETATION:
        std::cout << "植被" << std::endl;
        break;
      case SemanticLabel::BUILDING:
        std::cout << "建筑" << std::endl;
        break;
      case SemanticLabel::POLE:
        std::cout << "杆状" << std::endl;
        break;
      default:
        std::cout << "未知" << std::endl;
    }
  }
 
  return 0;
}

形状直径函数(SDF)

形状直径函数（Shape Diameter Function, SDF）是一种用于形状分析和分割的特征。它通过测量形状内部到表面的距离来捕捉局部形状特征。

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Point_set_3.h>
#include <CGAL/pca_estimate_normals.h>
 
#include <iostream>
#include <vector>
#include <cmath>
 
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
typedef Kernel::Vector_3 Vector;
typedef Kernel::FT FT;
 
typedef CGAL::Point_set_3<Point> Point_set;
 
// 简化的SDF计算（基于法线反方向的局部穿透深度估计）
double compute_sdf_approximation(
  const Point& query_point,
  const Vector& normal,
  const std::vector<Point>& neighbors,
  double max_distance = 1.0)
{
  if (neighbors.empty())
    return 0.0;
 
  // SDF估计：沿法线反方向搜索最近点
  Vector inward = -normal;
  inward = inward / std::sqrt(inward * inward);  // 归一化
 
  double min_distance = max_distance;
 
  for (const Point& neighbor : neighbors)
  {
    Vector diff = Vector(query_point, neighbor);
    double projection = diff * inward;
 
    // 只考虑法线反方向的点
    if (projection > 0)
    {
      double distance = std::sqrt(diff * diff);
      if (distance < min_distance)
        min_distance = distance;
    }
  }
 
  // 归一化SDF值
  return 1.0 - (min_distance / max_distance);
}
 
// 使用SDF进行形状分割
std::vector<int> sdf_segmentation(
  const Point_set& point_set,
  const std::vector<double>& sdf_values,
  double threshold = 0.5)
{
  std::vector<int> segments(point_set.size(), -1);
  int current_segment = 0;
 
  // 简化：基于SDF阈值分割
  for (std::size_t i = 0; i < sdf_values.size(); ++i)
  {
    if (sdf_values[i] > threshold)
      segments[i] = current_segment++;
  }
 
  return segments;
}
 
int main(int, char**)
{
  Point_set point_set;
  CGAL::Random random;
 
  // 生成简单形状（立方体）
  // 前面
  for (int i = 0; i < 100; ++i)
    point_set.insert(Point(random.get_double(-1., 1.),
                           random.get_double(-1., 1.),
                           1.0 + random.get_double(-0.01, 0.01)));
 
  // 后面
  for (int i = 0; i < 100; ++i)
    point_set.insert(Point(random.get_double(-1., 1.),
                           random.get_double(-1., 1.),
                           -1.0 + random.get_double(-0.01, 0.01)));
 
  // 估计法线
  point_set.add_normal_map();
  CGAL::pca_estimate_normals<CGAL::Sequential_tag>
    (point_set, 12, point_set.parameters());
 
  std::cout << "SDF计算示例" << std::endl;
  std::cout << "点云大小: " << point_set.size() << std::endl;
 
  // 实际应用中，SDF计算需要：
  // 1. 构建KD树进行高效邻域查询
  // 2. 沿法线方向发射射线
  // 3. 计算穿透深度
 
  return 0;
}

18.1.7 完整处理流程示例

以下示例展示了从原始点云到结构化输出的完整处理流程：

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Point_set_3.h>
#include <CGAL/Point_set_3/IO/XYZ.h>
 
// 点云处理算法
#include <CGAL/grid_simplify_point_set.h>
#include <CGAL/remove_outliers.h>
#include <CGAL/pca_estimate_normals.h>
#include <CGAL/mst_orient_normals.h>
#include <CGAL/jet_smooth_point_set.h>
#include <CGAL/compute_average_spacing.h>
 
// 形状检测
#include <CGAL/Shape_detection/Efficient_RANSAC.h>
 
#include <iostream>
#include <fstream>
#include <vector>
 
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
typedef Kernel::Vector_3 Vector;
typedef Kernel::Plane_3 Plane_3;
 
typedef CGAL::Point_set_3<Point> Point_set;
 
// 形状检测类型
typedef CGAL::Shape_detection::Efficient_RANSAC_traits
  <Kernel, Point_set, Point_set::Point_map, Point_set::Vector_map> Traits;
typedef CGAL::Shape_detection::Efficient_RANSAC<Traits> Efficient_ransac;
typedef CGAL::Shape_detection::Plane<Traits> Plane_shape;
typedef CGAL::Shape_detection::Sphere<Traits> Sphere_shape;
typedef CGAL::Shape_detection::Cylinder<Traits> Cylinder_shape;
 
// 处理统计
struct ProcessingStats
{
  std::size_t original_size;
  std::size_t after_simplification;
  std::size_t after_outlier_removal;
  std::size_t final_size;
  double average_spacing;
  std::size_t num_shapes;
  double processing_time;
};
 
// 完整处理流程
ProcessingStats process_point_cloud(Point_set& point_set)
{
  ProcessingStats stats;
  stats.original_size = point_set.size();
 
  std::cout << "=== 点云处理流程 ===" << std::endl;
  std::cout << "1. 原始点云: " << stats.original_size << " 个点" << std::endl;
 
  // 步骤1：计算平均点间距
  stats.average_spacing = CGAL::compute_average_spacing<CGAL::Sequential_tag>
    (point_set, 6, point_set.parameters());
  std::cout << "   平均点间距: " << stats.average_spacing << std::endl;
 
  // 步骤2：网格简化
  double cell_size = stats.average_spacing * 2.0;
  Point_set::iterator first_to_remove =
    CGAL::grid_simplify_point_set(point_set, cell_size,
                                   point_set.parameters());
  point_set.remove(first_to_remove, point_set.end());
  point_set.collect_garbage();
  stats.after_simplification = point_set.size();
  std::cout << "2. 网格简化后: " << stats.after_simplification << " 个点" << std::endl;
 
  // 步骤3：离群点移除
  Point_set::iterator first_outlier =
    CGAL::remove_outliers<CGAL::Sequential_tag>
    (point_set,
     24,  // 邻域大小
     point_set.parameters().
     threshold_percent(5.0).  // 最多移除5%
     threshold_distance(stats.average_spacing * 3.0));
 
  point_set.remove(first_outlier, point_set.end());
  point_set.collect_garbage();
  stats.after_outlier_removal = point_set.size();
  std::cout << "3. 离群点移除后: " << stats.after_outlier_removal << " 个点" << std::endl;
 
  // 步骤4：法线估计
  point_set.add_normal_map();
  CGAL::pca_estimate_normals<CGAL::Sequential_tag>
    (point_set, 12, point_set.parameters());
  std::cout << "4. 法线估计完成" << std::endl;
 
  // 步骤5：法线定向
  Point_set::iterator unoriented_begin =
    CGAL::mst_orient_normals(point_set, 12, point_set.parameters());
  point_set.remove(unoriented_begin, point_set.end());
  point_set.collect_garbage();
  std::cout << "5. 法线定向完成" << std::endl;
 
  // 步骤6：MLS平滑（可选）
  CGAL::jet_smooth_point_set<CGAL::Sequential_tag>
    (point_set, 18,
     point_set.parameters().
     degree_fitting(2).
     degree_monge(2));
  std::cout << "6. MLS平滑完成" << std::endl;
 
  // 步骤7：重新估计法线
  CGAL::pca_estimate_normals<CGAL::Sequential_tag>
    (point_set, 12, point_set.parameters());
 
  Point_set::iterator unoriented_begin2 =
    CGAL::mst_orient_normals(point_set, 12, point_set.parameters());
  point_set.remove(unoriented_begin2, point_set.end());
  point_set.collect_garbage();
  stats.final_size = point_set.size();
  std::cout << "7. 最终点云: " << stats.final_size << " 个点" << std::endl;
 
  return stats;
}
 
// 形状检测
void detect_shapes(Point_set& point_set, ProcessingStats& stats)
{
  std::cout << "\n=== 形状检测 ===" << std::endl;
 
  Efficient_ransac ransac;
  ransac.set_input(point_set,
                   point_set.point_map(),
                   point_set.normal_map());
 
  // 注册形状类型
  ransac.add_shape_factory<Plane_shape>();
  ransac.add_shape_factory<Sphere_shape>();
  ransac.add_shape_factory<Cylinder_shape>();
 
  // 设置参数
  Efficient_ransac::Parameters params;
  params.probability = 0.01;
  params.min_points = static_cast<std::size_t>(point_set.size() * 0.02);  // 2%
  params.epsilon = stats.average_spacing * 2.0;
  params.cluster_epsilon = stats.average_spacing * 5.0;
  params.normal_threshold = 0.9;
 
  // 检测
  ransac.detect(params);
 
  stats.num_shapes = ransac.shapes().size();
  std::cout << "检测到 " << stats.num_shapes << " 个形状" << std::endl;
 
  // 输出形状信息
  int plane_count = 0, sphere_count = 0, cylinder_count = 0;
  for (const auto& shape : ransac.shapes())
  {
    if (dynamic_cast<Plane_shape*>(shape.get()))
      plane_count++;
    else if (dynamic_cast<Sphere_shape*>(shape.get()))
      sphere_count++;
    else if (dynamic_cast<Cylinder_shape*>(shape.get()))
      cylinder_count++;
  }
 
  std::cout << "  平面: " << plane_count << std::endl;
  std::cout << "  球体: " << sphere_count << std::endl;
  std::cout << "  圆柱: " << cylinder_count << std::endl;
  std::cout << "  未分配: " << ransac.number_of_unassigned_points() << " 个点" << std::endl;
}
 
int main(int argc, char** argv)
{
  Point_set point_set;
 
  // 读取输入
  if (argc > 1)
  {
    std::ifstream in(argv[1]);
    if (!in)
    {
      std::cerr << "无法打开文件: " << argv[1] << std::endl;
      return 1;
    }
    in >> point_set;
    std::cout << "从文件读取点云: " << point_set.size() << " 个点" << std::endl;
  }
  else
  {
    // 生成测试数据
    CGAL::Random random;
 
    // 地面平面
    for (int i = 0; i < 500; ++i)
      point_set.insert(Point(random.get_double(-2., 2.),
                             random.get_double(-2., 2.),
                             random.get_double(-0.01, 0.01)));
 
    // 墙面
    for (int i = 0; i < 300; ++i)
      point_set.insert(Point(random.get_double(-0.02, 0.02),
                             random.get_double(-2., 2.),
                             random.get_double(0., 2.)));
 
    // 圆柱
    double cyl_radius = 0.3;
    Point cyl_center(1., 1., 1.);
    for (int i = 0; i < 200; ++i)
    {
      double theta = random.get_double(0., 2. * CGAL_PI);
      double z = random.get_double(0., 2.);
      double x = cyl_center.x() + cyl_radius * cos(theta);
      double y = cyl_center.y() + cyl_radius * sin(theta);
      point_set.insert(Point(x, y, z));
    }
 
    std::cout << "生成测试点云: " << point_set.size() << " 个点" << std::endl;
  }
 
  // 执行处理流程
  ProcessingStats stats = process_point_cloud(point_set);
 
  // 形状检测
  detect_shapes(point_set, stats);
 
  // 保存结果
  std::ofstream out("processed_point_set.xyz");
  out << point_set;
  std::cout << "\n结果已保存到 processed_point_set.xyz" << std::endl;
 
  // 输出统计摘要
  std::cout << "\n=== 处理统计摘要 ===" << std::endl;
  std::cout << "原始点数: " << stats.original_size << std::endl;
  std::cout << "最终点数: " << stats.final_size << std::endl;
  std::cout << "压缩率: " << (100.0 * stats.final_size / stats.original_size) << "%" << std::endl;
  std::cout << "检测形状: " << stats.num_shapes << std::endl;
 
  return 0;
}

18.1.8 性能优化

处理大规模点云时需要考虑性能优化策略。

并行处理

CGAL支持多线程并行处理：

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Point_set_3.h>
#include <CGAL/pca_estimate_normals.h>
#include <CGAL/jet_smooth_point_set.h>
 
#include <iostream>
#include <chrono>
 
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
 
typedef CGAL::Point_set_3<Point> Point_set;
 
int main(int, char**)
{
  Point_set point_set;
  CGAL::Random random;
 
  // 生成大规模点云
  for (int i = 0; i < 1000; ++i)
    for (int j = 0; j < 1000; ++j)
      point_set.insert(Point(i * 0.01, j * 0.01, random.get_double(-0.01, 0.01)));
 
  std::cout << "点云大小: " << point_set.size() << std::endl;
 
  point_set.add_normal_map();
 
  // 串行处理
  {
    Point_set test_set = point_set;
    auto start = std::chrono::high_resolution_clock::now();
 
    CGAL::pca_estimate_normals<CGAL::Sequential_tag>
      (test_set, 12, test_set.parameters());
 
    auto end = std::chrono::high_resolution_clock::now();
    auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
    std::cout << "串行法线估计: " << duration.count() << " ms" << std::endl;
  }
 
  // 并行处理（如果CGAL支持）
#ifdef CGAL_LINKED_WITH_TBB
  {
    Point_set test_set = point_set;
    auto start = std::chrono::high_resolution_clock::now();
 
    CGAL::pca_estimate_normals<CGAL::Parallel_tag>
      (test_set, 12, test_set.parameters());
 
    auto end = std::chrono::high_resolution_clock::now();
    auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
    std::cout << "并行法线估计: " << duration.count() << " ms" << std::endl;
  }
#else
  std::cout << "并行支持未启用" << std::endl;
#endif
 
  return 0;
}

内存优化

对于超大规模点云，可以考虑以下策略：

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Point_set_3.h>
 
#include <iostream>
#include <vector>
 
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
 
typedef CGAL::Point_set_3<Point> Point_set;
 
// 分块处理策略
void process_in_chunks(const std::string& filename,
                       std::size_t chunk_size = 100000)
{
  std::cout << "分块处理: " << filename << std::endl;
  std::cout << "每块大小: " << chunk_size << " 点" << std::endl;
 
  // 实际实现需要：
  // 1. 流式读取点云文件
  // 2. 分块处理每个子集
  // 3. 合并结果
}
 
// 使用32位索引减少内存占用（对于小于2^32点的点云）
void use_compact_indices()
{
  // 定义 CGAL_POINT_SET_3_USE_STD_SIZE_T_AS_SIZE_TYPE
  // 可以使用32位索引，每个索引节省4字节
  std::cout << "使用紧凑索引模式" << std::endl;
}
 
int main(int, char**)
{
  Point_set point_set;
 
  // 预分配内存
  point_set.reserve(1000000);
 
  // 添加点...
 
  // 及时回收垃圾
  point_set.collect_garbage();
 
  // 清理不需要的属性
  // point_set.clear_properties();
 
  return 0;
}

18.1.9 本章小结

本章介绍了CGAL中点云处理与分析的核心功能：

Point_set_3数据结构：

提供动态属性系统，支持法线、颜色、标签等任意属性
高效的内存管理和惰性删除机制
与CGAL算法无缝集成

法线估计：

PCA方法适用于一般情况
Jet拟合适用于曲面特征明显的数据
MST定向确保法线一致性

MLS平滑：

基于局部曲面拟合的平滑技术
参数选择影响平滑效果和特征保留

RANSAC形状检测：

检测平面、球体、圆柱、圆锥、圆环等形状
参数调优对检测质量至关重要
支持结构化点集生成

点云分类：

基于PCA的特征提取
规则分类器实现
SDF用于形状分析

性能优化：

并行处理加速计算
分块处理超大规模数据
内存管理策略

CGAL的点云处理框架提供了从原始数据到结构化输出的完整工具链，适用于各种应用场景。

18.1.10 练习

法线估计实验：使用不同邻域大小（k=6, 12, 24, 48）对同一组带噪声的点云进行法线估计，比较结果差异并分析原因。
MLS平滑参数调优：对包含锐利边缘的点云应用MLS平滑，调整邻域大小和多项式阶数，观察对边缘保留的影响。
形状检测实践：创建一个包含多个几何形状（平面、球体、圆柱）的合成点云，使用RANSAC检测所有形状并计算检测精度。
点云分类器：实现一个基于几何特征的简单分类器，能够区分地面、墙面和植被点。
完整处理流程：编写一个完整的点云处理程序，包含读取、简化、去噪、平滑、法线估计、形状检测和保存结果的全流程。
性能测试：比较串行和并行算法在不同规模点云上的性能差异，绘制加速比曲线。
自定义属性：扩展Point_set_3，添加时间戳、反射强度等自定义属性，并实现基于这些属性的过滤功能。

Chunibyo

Explorer

18.1 点云处理与分析

18.1 点云处理与分析

相关章节

18.1.1 点云数据简介

18.1.2 Point_set_3数据结构

基本使用

法线属性管理

自定义属性

空间索引集成

点集合并

18.1.3 法线估计

PCA局部平面拟合

法线方向一致性

不同邻域大小的影响

Jet拟合估计法线

18.1.4 移动最小二乘平滑

曲面重建原理

投影与平滑

参数选择指南

18.1.5 RANSAC形状检测

平面、圆柱、球体检测

参数设置策略

结果过滤与后处理

18.1.6 点云分类基础

特征提取

简单分类器实现

形状直径函数(SDF)

18.1.7 完整处理流程示例

18.1.8 性能优化

并行处理

内存优化

18.1.9 本章小结

18.1.10 练习

Graph View

Table of Contents

Backlinks

Chunibyo

Explorer

18.1 点云处理与分析

18.1 点云处理与分析 §

相关章节 §

18.1.1 点云数据简介 §

18.1.2 Point_set_3数据结构 §

基本使用 §

法线属性管理 §

自定义属性 §

空间索引集成 §

点集合并 §

18.1.3 法线估计 §

PCA局部平面拟合 §

法线方向一致性 §

不同邻域大小的影响 §

Jet拟合估计法线 §

18.1.4 移动最小二乘平滑 §

曲面重建原理 §

投影与平滑 §

参数选择指南 §

18.1.5 RANSAC形状检测 §

平面、圆柱、球体检测 §

参数设置策略 §

结果过滤与后处理 §

18.1.6 点云分类基础 §

特征提取 §

简单分类器实现 §

形状直径函数(SDF) §

18.1.7 完整处理流程示例 §

18.1.8 性能优化 §

并行处理 §

内存优化 §

18.1.9 本章小结 §

18.1.10 练习 §

Graph View

Table of Contents

Backlinks

18.1 点云处理与分析

相关章节

18.1.1 点云数据简介

18.1.2 Point_set_3数据结构

基本使用

法线属性管理

自定义属性

空间索引集成

点集合并

18.1.3 法线估计

PCA局部平面拟合

法线方向一致性

不同邻域大小的影响

Jet拟合估计法线

18.1.4 移动最小二乘平滑

曲面重建原理

投影与平滑

参数选择指南

18.1.5 RANSAC形状检测

平面、圆柱、球体检测

参数设置策略

结果过滤与后处理

18.1.6 点云分类基础

特征提取

简单分类器实现

形状直径函数(SDF)

18.1.7 完整处理流程示例

18.1.8 性能优化

并行处理

内存优化

18.1.9 本章小结

18.1.10 练习