pybind11 | 绑定CGAL几何算法(numpy数据交换)

一 前言

对于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为数组所需的维数。

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()对其的缓冲区信息进行访问。

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)之间的转换

// 将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函数。

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 绑定部分代码

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 示例完整代码

#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调用测试

测试代码如下:

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文件)实现的。

调试结果如下:

在这里插入图片描述

可以看到,函数输入为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)

  • 2
    点赞
  • 3
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值