pybind11 绑定CGAL几何算法（numpy数据交换）

[TOC]

一 前言

对于CGAL，前段时间也用过相关的Python绑定，详情见文章：【CGAL+Python】安装CGAL的Python绑定。如果我们想要调用[Polygon Mesh Processing](CGAL 5.5.1 - Polygon Mesh Processing: User Manual)中的算法，在不就地读取网格文件的前提下，需要在Python代码中构建CGAL的多边形网格对象，显得非常的不方便和蹩脚。正好，我在之前使用过libigl库的Python绑定，其数据通过numpy进行交换，即输入和输出都是numpy数组。于是，在对pybind11进行了一段时间的学习后，我开始尝试通过pybind11对CGAL的相关函数进行绑定生成Python调用接口，更重要的是使用numpy数组进行数据交换。

二 numpy数据交换

2.1 pybind11对numpy的支持

在pybind11/numpy.h头文件中，提供了对Numpy array的支持。我们可以通过py::array_t<T>来实现一个Numpy array。

py::array_t支持一些基于Numpy的API：

  • .dtype()返回数组元素的类型。
  • .strides()返回数组strides的指针。
  • .reshape({i, j, ...})返回指定shape的数组视图。.resize({})也可以。
  • .index_at(i, j, ...)获取数组指定索引的元素。

为了更高效地访问大型数组，pybind11提供了不带检查的代理类unchecked<N>和mutable_unchecked<N>，其中N为数组所需的维数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
m.def("sum_3d", [](py::array_t<double> x) {
    auto r = x.unchecked<3>(); // x must have ndim = 3; can be non-writeable
    double sum = 0;
    for (py::ssize_t i = 0; i < r.shape(0); i++)
        for (py::ssize_t j = 0; j < r.shape(1); j++)
            for (py::ssize_t k = 0; k < r.shape(2); k++)
                sum += r(i, j, k);
    return sum;
});
m.def("increment_3d", [](py::array_t<double> x) {
    auto r = x.mutable_unchecked<3>(); // Will throw if ndim != 3 or flags.writeable is false
    for (py::ssize_t i = 0; i < r.shape(0); i++)
        for (py::ssize_t j = 0; j < r.shape(1); j++)
            for (py::ssize_t k = 0; k < r.shape(2); k++)
                r(i, j, k) += 1.0;
}, py::arg().noconvert());

这两个代理类的区别就是：当只用从一个array_t<T>对象中读取数据时，我们使用unchecked<N>对它进行代理访问；当我们需要修改一个array_t<T>对象中的数据时，我们使用mutable_unchecked<N>对它进行代理访问。

同时Numpy array支持缓冲协议（buffer protocol），因此我们也可以通过.request()对其的缓冲区信息进行访问。

1
2
3
4
5
6
7
8
struct buffer_info {
    void *ptr;	// Pointer to the underlying storage
    py::ssize_t itemsize;	// Size of individual items in bytes
    std::string format;
    py::ssize_t ndim;		// Number of dimensions
    std::vector<py::ssize_t> shape; // Shape of the tensor (1 entry per dimension)
    std::vector<py::ssize_t> strides;	//Number of bytes between adjacent entries
};

其他详细信息参见文档：NumPy - pybind11 documentation。

2.2 Numpy VF(py::array_t)与CGAL mesh(Surface Mesh)之间的转换

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
// 将Numpy VF转换成CGAL mesh
CGAL_Mesh convert_mesh_from_Numpy_to_CGAL(py::array_t<double>& V, py::array_t<int>& F)
{
    CGAL_Mesh CGAL_mesh;
    // 获取V, F的信息
    py::buffer_info buf_v = V.request();
    py::buffer_info buf_f = F.request();
    if (buf_v.shape[1] != 3)
        throw std::runtime_error("vertex must be 3 dims ");
    if (buf_f.shape[1] != 3)
        throw std::runtime_error("face must be 3 dims ");

    auto v = V.unchecked<2>();
    auto f = F.unchecked<2>();

    // clear and reserve the mesh
    CGAL_mesh.clear();
    int vn = buf_v.shape[0];    // 顶点个数
    int fn = buf_f.shape[0];    // 面个数
    int e = 0;
    CGAL_mesh.reserve(vn, 2 * fn, e);

    //copy the vertices
    double x, y, z;
    for (py::ssize_t i = 0; i < v.shape(0); i++)
    {
        Point p;
        x = v(i, 0);
        y = v(i, 1);
        z = v(i, 2);
        p = Point(x, y, z);
        CGAL_mesh.add_vertex(p);
    }

    // copy the faces
    std::vector <int> vertices;
    for (py::ssize_t i = 0; i < f.shape(0); i++)
    {
        vertices.resize(3);
        vertices[0] = f(i, 0);
        vertices[1] = f(i, 1);
        vertices[2] = f(i, 2);
        CGAL_mesh.add_face(CGAL_Mesh::Vertex_index(vertices[0]),
            CGAL_Mesh::Vertex_index(vertices[1]),
            CGAL_Mesh::Vertex_index(vertices[2]));
    }

    return CGAL_mesh;
}

// 将CGAL mesh转换成Numpy VF
std::pair<py::array_t<double>, py::array_t<int>> convert_mesh_from_CGAL_to_Numpy(CGAL_Mesh& CGAL_mesh)
{
    //申请内存
    py::array_t<double> V = py::array_t<double>(CGAL_mesh.number_of_vertices() * 3);
    py::array_t<int> F = py::array_t<int>(CGAL_mesh.number_of_faces() * 3);
    std::pair<py::array_t<double>, py::array_t<int>> VF(V,F);

    V.resize({ int(CGAL_mesh.number_of_vertices()), 3 });
    F.resize({ int(CGAL_mesh.number_of_faces()), 3 });

    auto v = V.mutable_unchecked<2>();  // mutable_unchecked: can be writeable
    auto f = F.mutable_unchecked<2>();

    int i = 0;
    for (CGAL_Mesh::Vertex_index vd : vertices(CGAL_mesh))
    {
        v(i, 0) = CGAL_mesh.point(vd).x();
        v(i, 1) = CGAL_mesh.point(vd).y();
        v(i, 2) = CGAL_mesh.point(vd).z();
        i++;
    }

    i = 0;
    for (CGAL_Mesh::Face_index fd : faces(CGAL_mesh))
    {
        int j = 0;
        for (CGAL_Mesh::Vertex_index vd : CGAL::vertices_around_face(CGAL_mesh.halfedge(fd), CGAL_mesh))
        {
            f(i, j) = int(vd);
            j++;
        }
        i++;
    }

    return VF;
}

三 绑定CGAL算法示例

3.1 示例函数

在本文中我们尝试绑定[Polygon Mesh Processing](CGAL 5.5.1 - Polygon Mesh Processing: User Manual)中的isotropic_remeshing函数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
struct halfedge2edge
{
    halfedge2edge(const CGAL_Mesh& m, std::vector<edge_descriptor>& edges)
        : m_mesh(m), m_edges(edges)
    {}
    void operator()(const halfedge_descriptor& h) const
    {
        m_edges.push_back(edge(h, m_mesh));
    }
    const CGAL_Mesh& m_mesh;
    std::vector<edge_descriptor>& m_edges;
};

std::pair<py::array_t<double>, py::array_t<int>> isotropic_remeshing(py::array_t<double>& V, py::array_t<int>& F, double target_edge_length = 1, unsigned int nb_iter = 5)
{
    CGAL_Mesh mesh = convert_mesh_from_Numpy_to_CGAL(V, F);
    std::vector<edge_descriptor> border;
    PMP::border_halfedges(faces(mesh), mesh, boost::make_function_output_iterator(halfedge2edge(mesh, border)));
    PMP::split_long_edges(border, target_edge_length, mesh);
    PMP::isotropic_remeshing(faces(mesh), target_edge_length, mesh,
        CGAL::parameters::number_of_iterations(nb_iter)
        .protect_constraints(true));
    return convert_mesh_from_CGAL_to_Numpy(mesh);
}

3.2 绑定部分代码

1
2
3
4
5
6
7
8
PYBIND11_MODULE(numpy_cgal, m) {

    m.doc() = "The CGAL geometry algorithms which called via the Numpy array";
    m.def("isotropic_remeshing", &isotropic_remeshing, "remeshes a triangular mesh",
        py::arg("V"), py::arg("F"), 
        py::arg("target_edge_length") = 1, py::arg("nb_iter") = 5);

}

3.3 示例完整代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Polygon_mesh_processing/remesh.h>
namespace PMP = CGAL::Polygon_mesh_processing;
namespace py = pybind11;
typedef CGAL::Simple_cartesian<double>                          Kernal;
typedef Kernal::Point_3                                         Point;
typedef CGAL::Surface_mesh<Point>                               CGAL_Mesh;
typedef boost::graph_traits<CGAL_Mesh>::halfedge_descriptor     halfedge_descriptor;
typedef boost::graph_traits<CGAL_Mesh>::edge_descriptor         edge_descriptor;


// 将Numpy VF转换成CGAL mesh
CGAL_Mesh convert_mesh_from_Numpy_to_CGAL(py::array_t<double>& V, py::array_t<int>& F)
{
    CGAL_Mesh CGAL_mesh;
    // 获取V, F的信息
    py::buffer_info buf_v = V.request();
    py::buffer_info buf_f = F.request();
    // Validate
    if (buf_v.ndim != 2)
        throw std::runtime_error("Number of dimensions of `V` must be 2");

    if (buf_v.shape[1] != 3)
        throw std::runtime_error("Number of columns in `V` must be 3");

    if (buf_f.ndim != 2)
        throw std::runtime_error("Number of dimensions of `F` must be 2");

    if (buf_f.shape[1] != 3)
        throw std::runtime_error("Number of columns in `F` must be 3");

    auto v = V.unchecked<2>(); // unchecked: can be non - writeable
    auto f = F.unchecked<2>();

    // clear and reserve the mesh
    CGAL_mesh.clear();
    int vn = buf_v.shape[0];    // 顶点个数
    int fn = buf_f.shape[0];    // 面个数
    int e = 0;
    CGAL_mesh.reserve(vn, e, fn);

    //copy the vertices
    double x, y, z;
    for (py::ssize_t i = 0; i < v.shape(0); i++)
    {
        Point p;
        x = v(i, 0);
        y = v(i, 1);
        z = v(i, 2);
        p = Point(x, y, z);
        CGAL_mesh.add_vertex(p);
    }

    // copy the faces
    std::vector <int> vertices;
    for (py::ssize_t i = 0; i < f.shape(0); i++)
    {
        vertices.resize(3);
        vertices[0] = f(i, 0);
        vertices[1] = f(i, 1);
        vertices[2] = f(i, 2);
        CGAL_mesh.add_face(CGAL_Mesh::Vertex_index(vertices[0]),
            CGAL_Mesh::Vertex_index(vertices[1]),
            CGAL_Mesh::Vertex_index(vertices[2]));
    }

    return CGAL_mesh;
}

// 将CGAL mesh转换成Numpy VF
std::pair<py::array_t<double>, py::array_t<int>> convert_mesh_from_CGAL_to_Numpy(CGAL_Mesh& CGAL_mesh)
{
    //申请内存
    py::array_t<double> V = py::array_t<double>(CGAL_mesh.number_of_vertices() * 3);
    py::array_t<int> F = py::array_t<int>(CGAL_mesh.number_of_faces() * 3);
    std::pair<py::array_t<double>, py::array_t<int>> VF(V,F);

    V.resize({ int(CGAL_mesh.number_of_vertices()), 3 });
    F.resize({ int(CGAL_mesh.number_of_faces()), 3 });

    auto v = V.mutable_unchecked<2>();  // mutable_unchecked: can be writeable
    auto f = F.mutable_unchecked<2>();

    int i = 0;
    for (CGAL_Mesh::Vertex_index vd : vertices(CGAL_mesh))
    {
        v(i, 0) = CGAL_mesh.point(vd).x();
        v(i, 1) = CGAL_mesh.point(vd).y();
        v(i, 2) = CGAL_mesh.point(vd).z();
        i++;
    }

    i = 0;
    for (CGAL_Mesh::Face_index fd : faces(CGAL_mesh))
    {
        int j = 0;
        for (CGAL_Mesh::Vertex_index vd : CGAL::vertices_around_face(CGAL_mesh.halfedge(fd), CGAL_mesh))
        {
            f(i, j) = int(vd);
            j++;
        }
        i++;
    }

    return VF;
}

struct halfedge2edge
{
    halfedge2edge(const CGAL_Mesh& m, std::vector<edge_descriptor>& edges)
        : m_mesh(m), m_edges(edges)
    {}
    void operator()(const halfedge_descriptor& h) const
    {
        m_edges.push_back(edge(h, m_mesh));
    }
    const CGAL_Mesh& m_mesh;
    std::vector<edge_descriptor>& m_edges;
};

std::pair<py::array_t<double>, py::array_t<int>> isotropic_remeshing(py::array_t<double>& V, py::array_t<int>& F, double target_edge_length = 1, unsigned int nb_iter = 5)
{
    CGAL_Mesh mesh = convert_mesh_from_Numpy_to_CGAL(V, F);
    std::vector<edge_descriptor> border;
    PMP::border_halfedges(faces(mesh), mesh, boost::make_function_output_iterator(halfedge2edge(mesh, border)));
    PMP::split_long_edges(border, target_edge_length, mesh);
    PMP::isotropic_remeshing(faces(mesh), target_edge_length, mesh,
        CGAL::parameters::number_of_iterations(nb_iter)
        .protect_constraints(true));
    return convert_mesh_from_CGAL_to_Numpy(mesh);
}

// 绑定代码
PYBIND11_MODULE(numpy_cgal, m) 
{
    m.doc() = "The CGAL geometry algorithms which called via the Numpy array";
    m.def("isotropic_remeshing", &isotropic_remeshing, "remeshes a triangular mesh",
        py::arg("V"), py::arg("F"), 
        py::arg("target_edge_length") = 1, py::arg("nb_iter") = 5);
}

这里为了图方便将所有代码全写在了一个源文件里，正常来说应当将绑定代码和绑定对象分开。

四 编译生成和测试

4.1 编译生成pyd文件

在这个示例中，我是直接使用Visual Studio编译生成pyd文件。操作很简单，首先是按照这篇文章：pybind11学习 | VS2022下安装配置在VS中配置好pybind11，然后按照这篇文章：CGAL的安装与在VS中的配置在VS中配置好CGAL。在pybind11和CGAL都配置成功的前提下，生成解决方案即可得到pyd文件。如果电脑上没装Visual Studio，也可以尝试使用CMake进行构建，也是十分简单的，只需在这篇文章：pybind11学习 | 使用CMake构建系统并生成pyd文件的示例基础上在CMakeLists.txt中添加CGAL的库文件和头文件即可。

4.2 Python调用测试

测试代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
import igl
import numpy as np
from CGAL import CGAL_Polygon_mesh_processing
from CGAL.CGAL_Polyhedron_3 import Polyhedron_3
import os.path


def isotropic_remeshing_off(filename, targetedgelength, remeshIter):
    if os.path.isfile(filename + ".off"):
        V, F, _ = igl.read_off(filename + "_iso_remesh.off", False)
        return F, V
    P = Polyhedron_3(filename + ".off")
    flist = []
    for fh in P.facets():
        flist.append(fh)
    CGAL_Polygon_mesh_processing.isotropic_remeshing(flist, targetedgelength, P, remeshIter)
    P.write_to_file(filename + "_iso_remesh.off")
    V, F, _ = igl.read_off(filename + "_iso_remesh.off", False)
    return F, V


def isotropic_remeshing_bindings(V, F, targetedgelength, remeshIter):
    import numpy_cgal
    V, F = numpy_cgal.isotropic_remeshing(V, F, targetedgelength, remeshIter)
    return F, V


if __name__ == "__main__":
    v, f = igl.read_triangle_mesh("test_mesh.off", dtypef = np.float64)
    f_remeshed, v_remeshed = isotropic_remeshing_bindings(v, f, 1, 5)

其中，第一个isotropic_remeshing_off函数是通过文章【CGAL+Python】安装CGAL的Python绑定中提到的绑定加上读写off文件实现的CGAL重新网格化算法调用。第二个isotropic_remeshing_bindings函数是通过调用pybind11生成的Python拓展模块（即本文的方法，numpy_cgal.pyd为上一小节生成的pyd文件）实现的。

调试结果如下：

image-20221230143732408

可以看到，函数输入为ndarray类型，输出仍然为ndarray类型，且成功重新网格化，测试完毕。

五 总结

本文主要介绍一种在Python代码中调用C++第三方库API的思路。主要目的是抛砖引玉，本文的方法不一定是最好的，仅供大家参考学习。

参考和拓展

[1] CGAL 5.5.1 - Surface Mesh: User Manual

[2] CGAL 5.5.1 - Polygon Mesh Processing: User Manual

[3] pybind11 documentation

[4] pybind11—opencv图像处理(numpy数据交换) - 简书 (jianshu.com)

[5] pybind11-wrapped CGAL (github.com)

[6] cmpute/cgal.py: Pybind11 binding of cgal. Aimed at extending and easy use (github.com)