Delaunay triangulation on 2D data using meshpy

[thomasaarholt]

Hi! Thanks for writing this package!

I'm wanting to perform Delaunay triangulation on (large, ~10 million) 2D scatter points, for purposes of interpolation of new points. I was intending to benchmark meshpy / triangle's implementation, but I'm not sure how to get it to work!

With Scipy's Delaunay function (based on qhull), one gives it a number of points and it returns (among other things) a list of integers, where each row lists the indices of the three points that are connected to the point associated with the row index (that sentence came out possibly a bit confusing).

But how does one do something similar with meshpy?

points = np.array([[0, 0], [0, 1.1], [1, 0], [1, 1]])
from scipy.spatial import Delaunay
tri = Delaunay(points)
tri.vertices
array([[2, 3, 0],
       [3, 1, 0]], dtype=int32)

I can't seem to get any output when trying the following, and I am not sure how to proceed.

points = np.array([[0, 0], [0, 1.1], [1, 0], [1, 1]])

mesh_info = MeshInfo()
mesh_info.set_points(points)
mesh = build(mesh_info)
print(np.array(mesh.neighbors))
print(np.array(mesh.facets))
print(np.array(mesh.faces))
print(np.array(mesh.elements))
print(np.array(mesh.points))
print(np.array(mesh.element_volumes))
print(np.array(mesh.point_attributes))
# all return []

[inducer]

You may want to adjust the flags passed to build.

[inducer]

Cross-reference here: https://www.cs.cmu.edu/~quake/triangle.switch.html

[wk39]

In advance, thanks for the great library :)

What I did is...
first, set_points() : only 2D points supported
second, set_facet() : facets means point index pairs of boundary line segments (instead of Qhull in scipy, you have to specify boundary line by point indices)
finally, build() : matplotlib triplot by mesh.points and mesh.elements