Pybind11封装CGAL Delaunay三角化扩展段错误排查求助

我正在开发基于Pybind11的Python C++扩展，调用CGAL实现3D Delaunay三角化，为每个顶点附加索引信息并实现顶点增删改操作。运行Python脚本导出顶点时频繁出现段错误，仅偶尔正常；日志显示vertex_handle->info()获取到异常数值（如678559744），疑似数组越界。请求排查原因并提供解决方案，相关代码及报错如下：

核心C++代码 #

#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Triangulation_3.h>
#include <CGAL/Triangulation_vertex_base_with_info_3.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/convex_hull_3_to_face_graph.h>
#include <CGAL/Surface_mesh/IO/PLY.h>

#include <CGAL/compute_average_spacing.h>

#include <vector>
#include <unordered_map>

#include "utils/vec_math.h"

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Triangulation_vertex_base_with_info_3<unsigned int, K> Vb;
typedef CGAL::Triangulation_data_structure_3<Vb> Tds;
typedef CGAL::Delaunay_triangulation_3<K, Tds> Triangulation;

typedef K::FT FT;
typedef CGAL::Parallel_if_available_tag Concurrency_tag;

typedef Triangulation::Cell_handle Cell_handle;
typedef Triangulation::Vertex_handle Vertex_handle;
typedef Triangulation::Locate_type Locate_type;
typedef Triangulation::Point Point;

typedef CGAL::Surface_mesh<Point> Surface_mesh;

class TetrahedraBuilder {
public:
    TetrahedraBuilder();
    ~TetrahedraBuilder();

    void add_points(size_t num_points, float3 *points);
    void delete_points(size_t num_points, unsigned int *index);
    void move_points(size_t num_points, unsigned int *index, float3 *point);

    std::vector<float3> export_points_idx();

private:
    std::unordered_map<size_t, Vertex_handle> idx_to_vertex;
    Triangulation delaunay_tetra;
    unsigned int cur_size = 0;
};

/* 实现部分如下 */
TetrahedraBuilder::TetrahedraBuilder() {
    idx_to_vertex == std::unordered_map<size_t, Vertex_handle>();
    delaunay_tetra = Triangulation();
    cur_size = 0;
}

TetrahedraBuilder::~TetrahedraBuilder() {
}

void TetrahedraBuilder::add_points(size_t num_points, float3 *points) {
    /* 按顺序添加点 */
    for (int i = 0; i < num_points; i++) {
        auto p = Point(points[i].x, points[i].y, points[i].z);
        auto vertex_handle = delaunay_tetra.insert(p);
        vertex_handle->info() = cur_size;
        idx_to_vertex[cur_size] = vertex_handle;
        cur_size++;
    }
}

void TetrahedraBuilder::delete_points(size_t num_points, unsigned int *index) {
    /* 删除部分点，保留剩余点的连续索引 */
    // TODO：可优化
    std::unordered_set<unsigned int> uset(index, index + num_points);
    size_t new_size = 0;
    for  (int i = 0; i < cur_size; i++) {
        if (uset.find(i) == uset.end()) {
            auto vertex_handle = idx_to_vertex[i];
            vertex_handle->info() = new_size;
            idx_to_vertex[new_size] = vertex_handle;
            new_size++;
        }
    }
    cur_size = new_size;
}

void TetrahedraBuilder::move_points(size_t num_points, unsigned int *index, float3 *point) {
    for (int i = 0; i < num_points; i++) {
        auto vertex_handle = idx_to_vertex[index[i]];
        const auto p = Point(point[i].x, point[i].y, point[i].z);
        auto rt_vhandle = delaunay_tetra.move_if_no_collision(vertex_handle, p);
        if (rt_vhandle != vertex_handle) {
            std::cout << "目标位置已存在点！" << std::endl;
        }
    }
}
std::vector<float3> TetrahedraBuilder::export_points_idx() {
    std::vector<float3> result(cur_size);
    float *result_float = reinterpret_cast<float*>(result.data());
    std::cout<<"hello in here 1!" << std::endl;

    // size_t i = 0;
    for (auto vertex_handle : delaunay_tetra.all_vertex_handles()) {
        
        const auto point = vertex_handle->point();
        // result_float[vertex_handle->info()] = make_float3(point.x(), point.y(), point.z());
        const unsigned int i = vertex_handle->info();
        std::cout<<"hello in index " << i << std::endl;
        result_float[i * 3 + 0] = point.x();
        result_float[i * 3 + 1] = point.y();
        result_float[i * 3 + 2] = point.z();
    }
    return result;
}

Pybind11封装代码 #

struct PyTetrahedraBuilder {
public:
    PyTetrahedraBuilder() {
        builder = std::make_unique<TetrahedraBuilder>();
    }

    ~PyTetrahedraBuilder() {
        builder.reset();
    }
    void add_points(const torch::Tensor &points) {
        CHECK_CONTIGUOUS(points);
        CHECK_CPU(points);
        CHECK_FLOAT(points);
        TORCH_CHECK(points.dim() == 2 && points.size(1) == 3, "points must have shape [num_points, 3]");    // [num_points, 3]
        builder->add_points(
            points.numel() / 3,
            reinterpret_cast<float3 *>(points.data_ptr()));
    }
    torch::Tensor export_points_idx() {
        // std::cout<<"hello 1!" << std::endl;
        std::vector<float3> points = builder->export_points_idx();
        // std::cout<<"hello 2!" << std::endl;
        if (points.size() >= (size_t)std::numeric_limits<int>::max) {
            throw Exception("Too many points!");
        }
        // for (int i = 0; i < points.size(); i++) {
        //     std::cout<< i <<": " << points[i].x << " " << points[i].y << " " << points[i].z << std::endl;
        // }
        auto points_out = torch::empty({(long)points.size(), 3}, torch::dtype(torch::kFloat32).device(torch::kCPU));
        memcpy(
            points_out.data_ptr(),
            reinterpret_cast<void *>(points.data()),
            points.size() * sizeof(float3));
        return points_out;
    }
PYBIND11_MODULE(dynamic_tetra_cpp_extension, m) {
    py::class_<PyTetrahedraBuilder>(m, "TetrahedraBuilder")
        .def(py::init<>())
        .def("add_points", &PyTetrahedraBuilder::add_points)
        .def("export_points_idx", &PyTetrahedraBuilder::export_points_idx);
}

Python测试脚本 #

import numpy as np
import torch 
import trimesh


from utils.extension import TetrahedraTracer, TetrahedraBuilder

tetra_build = TetrahedraBuilder()

points = torch.rand(100, 3)
tetra_build.add_points(points)

vertices = tetra_build.export_points_idx()
print(vertices)

报错输出 #

hello in here 1!
hello in index 678559744
[1]    1369488 segmentation fault (core dumped)  python test_cpp_extension.py

问题原因分析 #

1. CGAL幽灵顶点未过滤：delaunay_tetra.all_vertex_handles()会遍历CGAL三角化中所有顶点，包括内部用于维护数据结构的幽灵顶点（ghost vertices）。这些顶点从未被初始化过info字段，因此读取到的是内存中的随机垃圾值，直接用来索引数组会导致越界。

2. 删除点逻辑无效：delete_points方法仅修改了idx_to_vertex哈希表和顶点的info值，但没有实际调用CGAL的remove方法从三角化中删除顶点。这导致三角化中仍然存在被标记为"删除"的顶点，遍历这些顶点时info值可能已被篡改或超出cur_size范围。

3. 构造函数语法错误：TetrahedraBuilder构造函数中使用了idx_to_vertex == std::unordered_map<size_t, Vertex_handle>();（双等号），这是一个无效的比较操作，并未初始化哈希表，可能导致后续哈希表操作出现未定义行为。

4. 数组越界无检查：export_points_idx中result的大小为cur_size，但遍历的顶点数量可能远大于这个值，且未对vertex_handle->info()的数值做边界校验，一旦索引超出cur_size-1就会触发段错误。

解决方案 #

1. 过滤无效顶点，遍历有效顶点集合

替换all_vertex_handles()为自己维护的idx_to_vertex哈希表，避免遍历幽灵顶点和无效顶点：

std::vector<float3> TetrahedraBuilder::export_points_idx() {
    std::vector<float3> result(cur_size);
    float *result_float = reinterpret_cast<float*>(result.data());
    std::cout<<"hello in here 1!" << std::endl;

    // 遍历自己维护的有效顶点映射
    for (const auto& pair : idx_to_vertex) {
        size_t idx = pair.first;
        Vertex_handle vertex_handle = pair.second;
        const auto point = vertex_handle->point();
        result_float[idx * 3 + 0] = point.x();
        result_float[idx * 3 + 1] = point.y();
        result_float[idx * 3 + 2] = point.z();
    }
    return result;
}

2. 正确实现顶点删除逻辑

调用CGAL的remove方法实际删除顶点，并同步清理哈希表：

void TetrahedraBuilder::delete_points(size_t num_points, unsigned int *index) {
    std::unordered_set<unsigned int> uset(index, index + num_points);
    // 先从CGAL三角化中删除目标顶点
    for (size_t i = 0; i < num_points; ++i) {
        unsigned int idx = index[i];
        auto it = idx_to_vertex.find(idx);
        if (it != idx_to_vertex.end()) {
            delaunay_tetra.remove(it->second);
            idx_to_vertex.erase(it);
        }
    }
    // 重新分配连续索引并更新映射
    size_t new_size = 0;
    std::unordered_map<size_t, Vertex_handle> new_idx_map;
    for (const auto& pair : idx_to_vertex) {
        pair.second->info() = new_size;
        new_idx_map[new_size] = pair.second;
        new_size++;
    }
    idx_to_vertex.swap(new_idx_map);
    cur_size = new_size;
}

注意：CGAL的remove操作会触发三角化重构，需确保操作的原子性。

3. 修复构造函数语法错误

将双等号改为单等号，或使用初始化列表：

// 方式1：初始化列表
TetrahedraBuilder::TetrahedraBuilder() 
    : idx_to_vertex(), delaunay_tetra(), cur_size(0) {
}

// 方式2：赋值初始化
TetrahedraBuilder::TetrahedraBuilder() {
    idx_to_vertex = std::unordered_map<size_t, Vertex_handle>();
    delaunay_tetra = Triangulation();
    cur_size = 0;
}

4. 添加边界校验

在export_points_idx中增加索引范围检查，避免越界：

const unsigned int i = vertex_handle->info();
if (i >= cur_size) {
    std::cerr << "Invalid vertex info index: " << i << ", cur_size: " << cur_size << std::endl;
    continue; // 或抛出异常终止程序
}
std::cout<<"hello in index " << i << std::endl;
result_float[i * 3 + 0] = point.x();
// ...

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