Can pytetgen be used for 2D triangulation?

[thomasaarholt]

The name seems to imply 3D, so I'm feeling a bit sheepish having tried the following, where I realised that pytetgen is interpreting my 10_000 2D points as 6666 3D points.

import numpy as np
import pytetgen
from scipy.spatial import Delaunay

N = 10000 # number of points
known_points = np.random.random((N,2))
known_intensities = np.random.random((N,))

tri1 = Delaunay(known_points) 
tri2 = pytetgen.Delaunay(known_points, neighbors=False)

print(tri1.points.shape)
# (10000, 2)
print(tri2.points.shape)
# (6666, 3)

Having done a bit of looking up (stackoverflow), it seems that tetgen is not what I should be looking at instead be using Triangle.

As a cheeky workaround I tried Delaunay on 3D points but with Z=0 for all points, but python crashes:

import numpy as np
import pytetgen
from scipy.interpolate import LinearNDInterpolator
from scipy.spatial import Delaunay

N = 1000 # number of points
known_points1 = np.random.random((N,3)) # 3D for tetgen, with z=0
known_points1[:,2] = 0
known_points2 = known_points1[:,:2] # 2D for scipy

known_intensities = np.random.random((N,))

tri1 = pytetgen.Delaunay(known_points1)
tri2 = Delaunay(known_points2) 

[Marcello-Sega]

Tetgen works in principle also for 2D point sets, but I haven't added that feature into the python wrapper. I'll have a look into that, if it's not too complicated I'll implement it.

[Marcello-Sega]

In fact, it seems that it handles only 3D sets, even though somewhere the opposite is stated, and there are pieces of code that actually (try to) handle the 2D case, although they look more like old code leftovers.

On the other hand, as you also mention, there's the triangle code, which is btw included in the same pdelib suite. I guess it would not be too difficult to add a wrapper to pytetgen. The question is whether it is faster or slower than scipy.spatial.Delaunay (which uses QHull under the hood).

Would you mind checking if that's the case? That would be really helpful: if it's faster, I'd gladly write the wrapper.

[thomasaarholt]

I will - my other comment stating that pytetgen was faster may not be correct, since I was testing it for the 2d case - but I'll be able to check tonight Europe time.

[thomasaarholt]

So I just checked with pytetgen for the 3D case, and it is about 5 times faster!

import numpy as np
import pytetgen
from scipy.spatial import Delaunay

N = 100000 # number of points
known_points = np.random.random((N,3))
known_intensities = np.random.random((N,))

tri1 = Delaunay(known_points) 
tri2 = pytetgen.Delaunay(known_points, neighbors=False)

%timeit Delaunay(known_points) 
# 2.3 sec
%timeit pytetgen.Delaunay(known_points)
# 0.46 sec

[thomasaarholt]

And triangle ( pip install triangle), is about 8 times faster. I run the triangulate with no optional flags, and it looks the same as scipy. I could probably use triangle as-is, but need to look at how scipy's interpolation function expects the delaunay input. Last time I had some trouble with that.

import triangle as tr
import numpy as np
from scipy.spatial import Delaunay

points = np.random.random((100000, 2))
triangle_intput = dict(vertices=points)

%timeit scipy_output = Delaunay(points)
490 ms ± 29.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit triangle_output = tr.triangulate(triangle_intput)
59 ms ± 553 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

The relevant output in the triangle_output dict is triangle_output['triangles']. It has the same length as `scipy., but is in a different order (I think this is fine).

And it looks the same:

# Generates the above figure
import numpy as np
import matplotlib.pyplot as plt
import triangle as tr
from scipy.spatial import Delaunay
points = np.random.random((12, 2))

triangle_input = dict(vertices=points)
triangle_output = tr.triangulate(triangle_input)

scipy_output = Delaunay(points)

fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(4,2))
ax1.set_title('scipy')
ax1.triplot(points[:,0], points[:,1], scipy_output.simplices)
ax1.plot(points[:,0], points[:,1], 'o')
ax1.axis('equal')

ax2.set_title('triangle')
tr.plot(ax2, vertices=triangle_output['vertices'], triangles=triangle_output['triangles'])
ax2.axis('equal');

[Marcello-Sega]

Oh, great, thanks a lot for checking this!

What I missed was, that there was already a triangle cython wrapper. Very good. I will then just use that within pytetgen to handle the 2D case, and provide a homogeneous output.

[thomasaarholt]

That's fantastic. I didn't realise you'd missed the cython wrapper - that made it a lot easier for me to check!

[Marcello-Sega]

I added the wrapper to triangle in https://github.com/Marcello-Sega/pytetgen/tree/2D-triangulation

Could you try it out and tell me if everything works as expected for you? In that case I'd merge to master and publish to PyPi (and maybe also conda forge)

[thomasaarholt]

Yep, that works nicely! Well done!