6

I have a geometry that is defined with a list of (x,y) points in space. I would like to create a triangular mesh out of this data, so I tried the Triangulation function in matplotlib for this purpose. However, since my geometry has some curves, the algorithm is generating undesired triangles between the edges of my part:

Image

Where the red curve is the edge of my geometry.

Is there any way to deal this issue? Maybe the Triangulation function is not what I need, in that case, do you have any recommendation on what to use?

The following code comes from this example. In the example they defined the triangles by explicitly naming three points instead of the Delaunay triangulation that I want to use by calling the function: triang = tri.Triangulation(x, y), and which will give me the same behavior than my original picture.

import matplotlib.pyplot as plt
import matplotlib.tri as tri
import numpy as np

xy = np.asarray([
    [-0.101, 0.872], [-0.080, 0.883], [-0.069, 0.888], [-0.054, 0.890],
    [-0.045, 0.897], [-0.057, 0.895], [-0.073, 0.900], [-0.087, 0.898],
    [-0.090, 0.904], [-0.069, 0.907], [-0.069, 0.921], [-0.080, 0.919],
    [-0.073, 0.928], [-0.052, 0.930], [-0.048, 0.942], [-0.062, 0.949],
    [-0.054, 0.958], [-0.069, 0.954], [-0.087, 0.952], [-0.087, 0.959],
    [-0.080, 0.966], [-0.085, 0.973], [-0.087, 0.965], [-0.097, 0.965],
    [-0.097, 0.975], [-0.092, 0.984], [-0.101, 0.980], [-0.108, 0.980],
    [-0.104, 0.987], [-0.102, 0.993], [-0.115, 1.001], [-0.099, 0.996],
    [-0.101, 1.007], [-0.090, 1.010], [-0.087, 1.021], [-0.069, 1.021],
    [-0.052, 1.022], [-0.052, 1.017], [-0.069, 1.010], [-0.064, 1.005],
    [-0.048, 1.005], [-0.031, 1.005], [-0.031, 0.996], [-0.040, 0.987],
    [-0.045, 0.980], [-0.052, 0.975], [-0.040, 0.973], [-0.026, 0.968],
    [-0.020, 0.954], [-0.006, 0.947], [ 0.003, 0.935], [ 0.006, 0.926],
    [ 0.005, 0.921], [ 0.022, 0.923], [ 0.033, 0.912], [ 0.029, 0.905],
    [ 0.017, 0.900], [ 0.012, 0.895], [ 0.027, 0.893], [ 0.019, 0.886],
    [ 0.001, 0.883], [-0.012, 0.884], [-0.029, 0.883], [-0.038, 0.879],
    [-0.057, 0.881], [-0.062, 0.876], [-0.078, 0.876], [-0.087, 0.872],
    [-0.030, 0.907], [-0.007, 0.905], [-0.057, 0.916], [-0.025, 0.933],
    [-0.077, 0.990], [-0.059, 0.993]])
x = np.degrees(xy[:, 0])
y = np.degrees(xy[:, 1])

triang = tri.Triangulation(x, y)
fig1, ax1 = plt.subplots()
ax1.set_aspect('equal')
ax1.triplot(triang, 'bo-', lw=1)
user190081
  • 222
  • 2
  • 7
  • Can you post the code that you used to produce the picture? – Thomas Kühn Sep 22 '18 at 15:27
  • There seems to be a maximum distance between your points. You may perform a Triangulation and remove all the triangles where one of the side lengths is larger than the maximum distance. Then plot this triangulation with `triplot`. – ImportanceOfBeingErnest Sep 22 '18 at 15:28
  • @ThomasKühn I cannot share my geometry, but I can reproduce the behavior with the example provided in the matplotlib documentation. I have edited my question to include this code. – user190081 Sep 22 '18 at 15:35
  • @ImportanceOfBeingErnest That's seems to be a good solution for me, but I dont know how to ask for the length of a particular edge. Can you guide me on that? – user190081 Sep 22 '18 at 15:39
  • 1
    `triag` should have an attribute `edges` which will be the index of the points of the edges and hence allow to reconstruct the triangles which can again be supplied to `Triangulation`. For now I can only promise to have a closer look at this later today. – ImportanceOfBeingErnest Sep 22 '18 at 16:24
  • You didn't pass any information about the edges to `Triangulation` so it uses the convex hull of the points. I don't know how such information could be provided. There is a `mask` argument but I don't really understand its purpose. I don't think the `edges` attribute can be of any help since `Triangulation` lacks any knowledge about what you call the edges. – Stop harming Monica Sep 22 '18 at 20:37
  • @Goyo yes, sorry `edges` was wrong, I meant the `triangles` attribute as now shown in the answer below. – ImportanceOfBeingErnest Sep 24 '18 at 00:51

2 Answers2

5

If you have the outline of the shape within to plot the triangulation you may apply the answer by @ThomasKühn.

Else, you may have a maximum distance between points above which triangles should not be taken into account. In that case you may mask those triangles out.

import matplotlib.pyplot as plt
import matplotlib.tri as tri
import numpy as np

xy = np.asarray([
    [-0.101, 0.872], [-0.080, 0.883], [-0.069, 0.888], [-0.054, 0.890],
    [-0.045, 0.897], [-0.057, 0.895], [-0.073, 0.900], [-0.087, 0.898],
    [-0.090, 0.904], [-0.069, 0.907], [-0.069, 0.921], [-0.080, 0.919],
    [-0.073, 0.928], [-0.052, 0.930], [-0.048, 0.942], [-0.062, 0.949],
    [-0.054, 0.958], [-0.069, 0.954], [-0.087, 0.952], [-0.087, 0.959],
    [-0.080, 0.966], [-0.085, 0.973], [-0.087, 0.965], [-0.097, 0.965],
    [-0.097, 0.975], [-0.092, 0.984], [-0.101, 0.980], [-0.108, 0.980],
    [-0.104, 0.987], [-0.102, 0.993], [-0.115, 1.001], [-0.099, 0.996],
    [-0.101, 1.007], [-0.090, 1.010], [-0.087, 1.021], [-0.069, 1.021],
    [-0.052, 1.022], [-0.052, 1.017], [-0.069, 1.010], [-0.064, 1.005],
    [-0.048, 1.005], [-0.031, 1.005], [-0.031, 0.996], [-0.040, 0.987],
    [-0.045, 0.980], [-0.052, 0.975], [-0.040, 0.973], [-0.026, 0.968],
    [-0.020, 0.954], [-0.006, 0.947], [ 0.003, 0.935], [ 0.006, 0.926],
    [ 0.005, 0.921], [ 0.022, 0.923], [ 0.033, 0.912], [ 0.029, 0.905],
    [ 0.017, 0.900], [ 0.012, 0.895], [ 0.027, 0.893], [ 0.019, 0.886],
    [ 0.001, 0.883], [-0.012, 0.884], [-0.029, 0.883], [-0.038, 0.879],
    [-0.057, 0.881], [-0.062, 0.876], [-0.078, 0.876], [-0.087, 0.872],
    [-0.030, 0.907], [-0.007, 0.905], [-0.057, 0.916], [-0.025, 0.933],
    [-0.077, 0.990], [-0.059, 0.993]])
x = np.degrees(xy[:, 0])
y = np.degrees(xy[:, 1])

triang = tri.Triangulation(x, y)

fig1, ax1 = plt.subplots()
ax1.set_aspect('equal')

# plot all triangles
ax1.triplot(triang, 'bo-', lw=0.2)

# plot only triangles with sidelength smaller some max_radius
max_radius = 2
triangles = triang.triangles

# Mask off unwanted triangles.
xtri = x[triangles] - np.roll(x[triangles], 1, axis=1)
ytri = y[triangles] - np.roll(y[triangles], 1, axis=1)
maxi = np.max(np.sqrt(xtri**2 + ytri**2), axis=1)
triang.set_mask(maxi > max_radius)

ax1.triplot(triang, color="indigo", lw=2.6)


plt.show()

The narrow lines show all triangles (the convex hull of the points), the bold lines show only those triangles where no side-length is larger then some maximum values (in this case chosen to be 2).

enter image description here

This thread may equally be relevant: matplotlib contour/contourf of **concave** non-gridded data

ImportanceOfBeingErnest
  • 321,279
  • 53
  • 665
  • 712
4

If the shape of the geometry is well defined, say by a curve, one can check for each triangle whether it lies within the shape or not. One can then define a mask and mask out undesired triangles. I found a solution using shapely, where define a polygon for the original shape (outline) and a polygon for each resulting triangle of Triangulation(), for which I then check whether it lies inside outline or not:

import matplotlib.pyplot as plt
import matplotlib.tri as tri
import numpy as np

import shapely
from shapely.geometry import Polygon as sPolygon


fig, (ax1,ax2) = plt.subplots(ncols=2)
ax1.set_aspect('equal')
ax2.set_aspect('equal')

##setting up basic shape
phi = np.linspace(0,2*np.pi,20)
r = 1 + 2*np.sin(phi)**2
x = np.cos(phi)*r
y = np.sin(phi)*r
ax1.plot(x,y,'ro-', lw=3, ms=6, zorder= 1, label='edge')
ax2.plot(x,y,'ro-', lw=3, ms=6, zorder= 1)


##original triangulation
triang1 = tri.Triangulation(x, y)
ax1.triplot(triang1, 'ko--', lw=1, ms=4, zorder=2, label='all')

##masking
outline = sPolygon(zip(x,y))
mask = [
    not outline.contains(sPolygon(zip(x[tri], y[tri])))
    for tri in triang1.get_masked_triangles()
]
triang1.set_mask(mask)
ax1.triplot(triang1, 'b-', lw=1, zorder=3, label='inner')

##adding more points
x_extra = np.random.rand(30)*(x.max()-x.min())+x.min()
y_extra = np.random.rand(30)*(y.max()-y.min())+y.min()

x = np.concatenate([x,x_extra])
y = np.concatenate([y,y_extra])

triang2 = tri.Triangulation(x,y)
ax2.triplot(triang2, 'ko--', lw=1, ms=4,  zorder=2)

##masking
mask = [
    not outline.contains(sPolygon(zip(x[tri], y[tri])))
    for tri in triang2.get_masked_triangles()
]
triang2.set_mask(mask)
ax2.triplot(triang2, 'b-', lw=1, zorder=3)

fig.legend()
plt.show()

The result of the code looks something like this:

result of above code

I wasn't quite sure what the OP wants, so on the left side I use only the points of the edge, while on the right side I added some random extra points for triangulation. In the picture, the outline of the shape is plotted in red, the original result of the Delaunay triangulation is plotted with a black dashed line, and the masked triangulation is plotted in blue.

EDIT:

I just noticed, that apparently one of the outline points is not included on the right side picture after the filtering. This must be due to numerical inaccuracies. One way to get around this is to increase the size of the outline slightly with the buffer() command. Something like this seems appropriate for the problem at hand:

outline = sPolygon(zip(x,y)).buffer(.01)

but the actual amount of buffering probably has to be adjusted.

Thomas Kühn
  • 9,412
  • 3
  • 47
  • 63