Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

How to deal with the (undesired) triangles that form between the edges of my geometry when using Triangulation in matplotlib

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)
like image 656
user190081 Avatar asked Sep 22 '18 15:09

user190081


2 Answers

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

like image 135
ImportanceOfBeingErnest Avatar answered Sep 29 '22 22:09

ImportanceOfBeingErnest


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.

like image 40
Thomas Kühn Avatar answered Sep 30 '22 00:09

Thomas Kühn