CGAL项目GIS教程：从点云到地形模型与分类

CGAL项目GIS教程：从点云到地形模型与分类

cgal The public CGAL repository, see the README below 项目地址: https://gitcode.com/gh_mirrors/cg/cgal

前言

在GIS（地理信息系统）应用中，处理密集点云数据是一项常见任务。本文将基于CGAL库，详细介绍如何从原始点云数据出发，构建数字表面模型(DSM)、数字地形模型(DTM)，并进行地形分类处理。本教程适合有一定CGAL基础，希望了解GIS相关处理的开发者。

核心概念

在开始前，我们需要明确几个关键术语：

• TIN（不规则三角网）：基于水平面投影连接的2D三角网格结构，但保留3D点坐标信息
• DSM（数字表面模型）：包含建筑物、植被等所有地表特征的表面模型
• DTM（数字地形模型）：仅包含裸露地面的地形模型
• DEM（数字高程模型）：包含DSM和DTM的通用术语

1. 构建不规则三角网(TIN)

CGAL提供了多种三角剖分数据结构和算法。我们可以通过结合2D Delaunay三角剖分与投影特性来构建TIN：

using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
using Projection_traits = CGAL::Projection_traits_xy_3<Kernel>;
using Point_2 = Kernel::Point_2;
using Point_3 = Kernel::Point_3;
using Segment_3 = Kernel::Segment_3;

using TIN = CGAL::Delaunay_triangulation_2<Projection_traits>;

2. 创建数字表面模型(DSM)

从点云数据创建DSM的过程非常直接：

std::ifstream ifs("data.xyz");
CGAL::Point_set_3<Point_3> points;
ifs >> points;

TIN tin;
tin.insert(points.points().begin(), points.points().end());

由于CGAL的Delaunay三角剖分符合FaceGraph概念，我们可以轻松将其转换为网格结构并保存：

std::ofstream dsm_ofs("dsm.ply");
CGAL::Surface_mesh<Point_3> dsm;
CGAL::copy_face_graph(tin, dsm);
CGAL::IO::write_PLY(dsm_ofs, dsm);

3. 生成数字地形模型(DTM)

3.1 扩展TIN数据结构

为了支持后续处理，我们需要扩展TIN数据结构，存储额外信息：

using Vbi = std::size_t;
using Fbi = bool;
using TIN_with_info = CGAL::Triangulation_2<Projection_traits, 
                                           CGAL::Triangulation_vertex_base_with_info_2<Vbi, Projection_traits>,
                                           CGAL::Triangulation_face_base_with_info_2<Fbi, Projection_traits>>;

3.2 识别连通组件

通过洪水填充算法识别连通组件：

double height_threshold = 0.5; // 建筑物最小高度阈值
std::vector<std::size_t> component(tin.number_of_faces(), -1);
std::size_t current_component = 0;

for (auto face_handle : tin.finite_face_handles()) {
    if (component[face_handle->info()] != -1)
        continue;
    
    std::queue<TIN_with_info::Face_handle> queue;
    queue.push(face_handle);
    
    while (!queue.empty()) {
        auto current_face = queue.front();
        queue.pop();
        
        if (face_height(current_face) > height_threshold)
            continue;
            
        component[current_face->info()] = current_component;
        
        for (int i = 0; i < 3; ++i) {
            auto neighbor = current_face->neighbor(i);
            if (component[neighbor->info()] == -1)
                queue.push(neighbor);
        }
    }
    ++current_component;
}

3.3 过滤小组件

移除小于建筑物最大尺寸的组件：

double perimeter_threshold = 50; // 建筑物最大周长阈值
std::vector<std::size_t> component_size(current_component, 0);
for (auto face_handle : tin.finite_face_handles()) {
    if (component[face_handle->info()] != -1)
        component_size[component[face_handle->info()]]++;
}

std::size_t largest_component = *std::max_element(component_size.begin(), 
                                                component_size.end());

3.4 孔洞填充与网格优化

// 过滤过大面片
CGAL::Surface_mesh<Point_3> dtm;
std::unordered_map<TIN_with_info::Vertex_handle, 
                  CGAL::Surface_mesh<Point_3>::Vertex_index> vertex_map;

for (auto face_handle : tin_with_info.finite_face_handles()) {
    if (face_height(face_handle) < height_threshold) {
        std::vector<CGAL::Surface_mesh<Point_3>::Vertex_index> vertices;
        for (int i = 0; i < 3; ++i) {
            auto vh = face_handle->vertex(i);
            if (vertex_map.find(vh) == vertex_map.end())
                vertex_map[vh] = dtm.add_vertex(vh->point());
            vertices.push_back(vertex_map[vh]);
        }
        dtm.add_face(vertices);
    }
}

// 孔洞填充
unsigned int nb_holes = 0;
for (auto hole : CGAL::holes(dtm)) {
    if (hole != largest_hole) {
        CGAL::triangulate_and_refine_hole(dtm, hole, std::back_inserter(dtm));
        ++nb_holes;
    }
}

// 各向同性网格重划分
CGAL::Polygon_mesh_processing::isotropic_remeshing(
    faces(dtm), 0.02, dtm,
    CGAL::parameters::number_of_iterations(3));

4. 栅格化处理

将TIN转换为栅格高度图：

double spacing = 0.5; // 栅格间距
CGAL::Bbox_2 bbox = CGAL::bbox_2(points.points().begin(), points.points().end());

std::size_t width = std::size_t((bbox.xmax() - bbox.xmin()) / spacing);
std::size_t height = std::size_t((bbox.ymax() - bbox.ymin()) / spacing);

std::vector<float> raster(width * height, -std::numeric_limits<float>::infinity());

for (std::size_t i = 0; i < width; ++i) {
    for (std::size_t j = 0; j < height; ++j) {
        Point_2 query(bbox.xmin() + i * spacing, 
                     bbox.ymin() + j * spacing);
        auto face_handle = tin.locate(query);
        if (tin.is_infinite(face_handle))
            continue;
            
        auto bary = CGAL::Barycentric_coordinates::compute_coordinates(
            face_handle->vertex(0)->point(), 
            face_handle->vertex(1)->point(),
            face_handle->vertex(2)->point(),
            Point_3(query.x(), query.y(), 0));
            
        if (bary[0] >= 0 && bary[1] >= 0 && bary[2] >= 0) {
            float z = bary[0] * face_handle->vertex(0)->point().z() +
                      bary[1] * face_handle->vertex(1)->point().z() +
                      bary[2] * face_handle->vertex(2)->point().z();
            raster[i + j * width] = z;
        }
    }
}

5. 等高线提取

5.1 构建等高线图

using Graph = boost::adjacency_list<boost::listS, boost::vecS, 
                                  boost::undirectedS,
                                  Point_3>;
using vertex_descriptor = boost::graph_traits<Graph>::vertex_descriptor;
using edge_descriptor = boost::graph_traits<Graph>::edge_descriptor;

Graph graph;
std::map<Point_3, vertex_descriptor> point_to_vertex;

double min_height = ...; // 计算最小高度
double max_height = ...; // 计算最大高度
int num_isolines = 50;
double height_increment = (max_height - min_height) / num_isolines;

for (int i = 0; i <= num_isolines; ++i) {
    double height = min_height + i * height_increment;
    
    for (auto face_handle : tin.finite_face_handles()) {
        if (face_has_isoline(face_handle, height)) {
            Segment_3 segment = extract_isoline(face_handle, height);
            
            vertex_descriptor v1, v2;
            auto it = point_to_vertex.find(segment.source());
            if (it == point_to_vertex.end()) {
                v1 = boost::add_vertex(segment.source(), graph);
                point_to_vertex[segment.source()] = v1;
            } else {
                v1 = it->second;
            }
            
            // 类似处理v2...
            boost::add_edge(v1, v2, graph);
        }
    }
}

5.2 简化和优化等高线

std::vector<std::vector<Point_3>> contours;
Contour_visitor<Point_3> visitor(contours);
CGAL::split_graph_into_polylines(graph, visitor);

// 使用约束三角剖分进行简化
using CDT = CGAL::Constrained_triangulation_plus_2<
    CGAL::Constrained_Delaunay_triangulation_2<Kernel>>;

double tolerance = 4.0 * average_spacing;
for (auto& polyline : contours) {
    CDT cdt;
    cdt.insert_constraint(polyline.begin(), polyline.end());
    
    CGAL::Polyline_simplification_2::simplify(
        cdt, CGAL::Polyline_simplification_2::Squared_distance_cost(),
        CGAL::Polyline_simplification_2::Stop_above_cost_threshold(tolerance * tolerance));
    
    // 更新简化后的折线...
}

6. 点云分类

CGAL提供了分类功能，可以将点云分割为用户定义的标签集：

using Point_set = CGAL::Point_set_3<Point_3>;
Point_set points;

// 加载点云和训练数据
std::ifstream ifs("classified.ply");
ifs >> points;

auto label_map = points.property_map<int>("label").first;
auto training_map = points.property_map<int>("training").first;

// 设置分类器
CGAL::Classification::ETHZ::Random_forest_classifier classifier(
    CGAL::Classification::Label_set({"ground", "vegetation", "roof"}));

// 训练分类器
classifier.train(points.range(), training_map, label_map);

// 执行分类
CGAL::Classification::classify_with_graphcut<CGAL::Parallel_if_available_tag>(
    points.range(), classifier, 0.2f, 10);

总结

本教程详细介绍了使用CGAL处理GIS数据的完整流程，从原始点云到最终的地形分类。通过结合CGAL的多种算法和数据结构，我们可以构建强大的GIS处理管线。实际应用中，可能需要根据具体数据特点调整参数和算法细节。

CGAL的强大之处在于其模块化设计和算法的高效实现，使得开发者可以灵活组合不同功能来满足特定需求。对于GIS应用开发者来说，掌握这些核心技术将大大提升处理地理空间数据的能力。

cgal The public CGAL repository, see the README below 项目地址: https://gitcode.com/gh_mirrors/cg/cgal