如何在C++中对带孔洞的2D点云进行网格划分？

一、从Alpha Shapes提取内外边界多边形 #

CGAL的Alpha Shapes没有直接返回多边形列表的接口，但可以通过遍历边界边手动组装闭合多边形，步骤如下：

• 初始化Alpha Shapes对象，可通过Alpha_shape_2::find_optimal_alpha()自动计算合适的alpha值，或根据点云密度手动调整
• 遍历所有有限边，筛选出边界边（调用is_boundary()方法判断）
• 从任意未访问的边界边开始，沿着相邻边界边遍历，直到回到起点，组装成闭合多边形
• 通过多边形面积符号区分外边界与孔洞：面积为正对应逆时针外边界，面积为负对应顺时针孔洞（CGAL约束三角剖分要求孔洞为顺时针）

关键代码片段：

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Alpha_shape_2.h>
#include <vector>
#include <unordered_set>

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_2 Point;
typedef CGAL::Alpha_shape_2<K> Alpha_shape_2;
typedef Alpha_shape_2::Alpha_shape_edge Alpha_shape_edge;

std::vector<std::vector<Point>> extract_boundary_polygons(Alpha_shape_2& alpha_shape) {
    std::vector<std::vector<Point>> polygons;
    std::unordered_set<Alpha_shape_edge> visited_edges;

    for (auto it = alpha_shape.finite_edges_begin(); it != alpha_shape.finite_edges_end(); ++it) {
        if (alpha_shape.is_boundary(*it) && visited_edges.find(*it) == visited_edges.end()) {
            std::vector<Point> polygon;
            Alpha_shape_edge current_edge = *it;
            do {
                Point p1 = alpha_shape.point(current_edge.first, current_edge.second);
                Point p2 = alpha_shape.point(current_edge.first, current_edge.third);
                if (polygon.empty() || polygon.back() != p1) polygon.push_back(p1);
                polygon.push_back(p2);
                visited_edges.insert(current_edge);
                current_edge = alpha_shape.next_boundary_edge(current_edge);
            } while (current_edge != *it);
            if (!polygon.empty()) polygon.pop_back(); // 移除闭合重复点
            polygons.push_back(polygon);
        }
    }
    return polygons;
}

二、用约束三角剖分生成带孔洞的网格 #

拿到边界多边形后，使用CGAL约束Delaunay三角剖分生成符合要求的网格：

• 初始化带面信息的约束三角剖分对象，用于标记内部/外部区域
• 添加外边界的所有边作为约束（确保逆时针）
• 添加每个孔洞的所有边作为约束（确保顺时针）
• 标记内部区域，遍历筛选出内部三角形

关键代码片段：

#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Triangulation_face_base_with_info_2.h>
#include <list>

struct FaceInfo2 {
    int nesting_level = -1;
    bool in_domain() { return nesting_level % 2 == 1; }
};

typedef CGAL::Triangulation_vertex_base_2<K> Vb;
typedef CGAL::Triangulation_face_base_with_info_2<FaceInfo2, K> Fb;
typedef CGAL::Constrained_triangulation_face_base_2<K, Fb> CFb;
typedef CGAL::Triangulation_data_structure_2<Vb, CFb> Tds;
typedef CGAL::Constrained_Delaunay_triangulation_2<K, Tds> CDT;
typedef CDT::Face_handle Face_handle;

void mark_domains(CDT& cdt, Face_handle start, int index, std::list<CDT::Edge>& border) {
    if (start->info().nesting_level != -1) return;
    std::list<Face_handle> queue;
    queue.push_back(start);
    start->info().nesting_level = index;
    while (!queue.empty()) {
        Face_handle fh = queue.front();
        queue.pop_front();
        for (int i = 0; i < 3; ++i) {
            CDT::Edge e(fh, i);
            Face_handle n = fh->neighbor(i);
            if (n->info().nesting_level == -1) {
                cdt.is_constrained(e) ? border.push_back(e) : (n->info().nesting_level = fh->info().nesting_level, queue.push_back(n));
            }
        }
    }
}

void mark_domains(CDT& cdt) {
    for (auto fh = cdt.finite_faces_begin(); fh != cdt.finite_faces_end(); ++fh) fh->info().nesting_level = -1;
    std::list<CDT::Edge> border;
    mark_domains(cdt, cdt.infinite_face()->neighbor(0), 0, border);
    while (!border.empty()) {
        CDT::Edge e = border.front();
        border.pop_front();
        Face_handle n = e.first->neighbor(e.second);
        if (n->info().nesting_level == -1) mark_domains(cdt, n, e.first->info().nesting_level + 1, border);
    }
}

std::vector<std::vector<Point>> generate_triangle_mesh(std::vector<std::vector<Point>>& polygons) {
    CDT cdt;
    // 添加外边界（假设第一个多边形为外边界）
    auto& outer_poly = polygons[0];
    for (size_t i = 0; i < outer_poly.size(); ++i) {
        size_t j = (i + 1) % outer_poly.size();
        cdt.insert_constraint(outer_poly[i], outer_poly[j]);
    }
    // 添加孔洞
    for (size_t p = 1; p < polygons.size(); ++p) {
        auto& hole_poly = polygons[p];
        for (size_t i = 0; i < hole_poly.size(); ++i) {
            size_t j = (i + 1) % hole_poly.size();
            cdt.insert_constraint(hole_poly[i], hole_poly[j]);
        }
    }
    // 标记内部区域并提取三角形
    mark_domains(cdt);
    std::vector<std::vector<Point>> triangles;
    for (auto fh = cdt.finite_faces_begin(); fh != cdt.finite_faces_end(); ++fh) {
        if (fh->info().in_domain()) {
            triangles.push_back({fh->vertex(0)->point(), fh->vertex(1)->point(), fh->vertex(2)->point()});
        }
    }
    return triangles;
}

注意事项 #

• Alpha值需适配点云密度：太小会生成细碎边界，太大可能合并孔洞
• 必须保证多边形方向正确：外边界逆时针，孔洞顺时针，否则剖分无法识别孔洞
• 若点云密度不均，可先对稀疏区域补点，再进行后续处理

内容的提问来源于stack exchange，提问作者Hugal31