文章目录
- 一 前言
- 二 numpy数据交换
- 2.1 pybind11对numpy的支持
- 2.2 Numpy VF(py::array_t)与CGAL mesh(Surface Mesh)之间的转换
- 三 绑定CGAL算法示例
- 3.1 示例函数
- 3.2 绑定部分代码
- 3.3 示例完整代码
- 四 编译生成和测试
- 4.1 编译生成pyd文件
- 4.2 Python调用测试
- 五 总结
- 参考和拓展
一 前言
对于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)