Last active
February 27, 2019 13:42
-
-
Save jwass/ea225c99676adb2a2be9 to your computer and use it in GitHub Desktop.
Example using proposed Delaunay triangulation method in Shapely
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
import matplotlib.pyplot as plt | |
import numpy as np | |
from shapely.geometry import LineString | |
from shapely.ops import triangulate | |
def uniform_sample(poly, n=100): | |
""" | |
Uniformly sample the Delaunay triangulation of a polygon. If the polygon | |
is convex, this will uniformly sample its area. | |
Parameters | |
---------- | |
poly: Shapely Polygon | |
n: Number of points | |
Returns | |
------- | |
[n x 2] numpy array of x-y coordinates that are uniformly distributed | |
over the Delaunay triangulation. | |
""" | |
polys = triangulate(poly) | |
# Normalize the areas | |
areas = np.array([p.area for p in polys]) | |
areas /= areas.sum() | |
# Randomly select a triangle weighted by area | |
# t_inds is the index of the chosen triangle | |
t_inds = np.searchsorted(np.cumsum(areas), np.random.random(n)) | |
# Randomly sample the area of each triangle according to | |
# P = (1-sqrt(r1))A + (sqrt(r1)(1-r2))B + (sqrt(r1)r2)C | |
# where r1, r2 are sampled uniform [0, 1] and A, B, C are the triangle | |
# vertices | |
# http://math.stackexchange.com/questions/18686/uniform-random-point-in-triangle | |
# Compute the coefficients | |
sr1 = np.sqrt(np.random.random(n)) # sr1 is sqrt(r1) above | |
r2 = np.random.random(n) | |
c0 = 1 - sr1 | |
c1 = sr1 * (1 - r2) | |
c2 = sr1 * r2 | |
# Grab the triangle vertices. | |
# v is a 3-element list where each element is [len(polys) x 2] | |
# array of each triangle vertex. It represents, A, B, C above. | |
v = [np.array([p.exterior.coords[i] for p in polys]) | |
for i in range(3)] | |
# Compute the points. v[i] is [N x 2] and the coefficients are [N x 1] | |
P = (c0[:, np.newaxis] * v[0][t_inds,:] + | |
c1[:, np.newaxis] * v[1][t_inds,:] + | |
c2[:, np.newaxis] * v[2][t_inds,:]) | |
return P | |
# Generate a random convex polygon by taking the convex hull | |
# of a bunch of random points | |
n = 8 | |
points = np.random.random((n,2)) | |
line = LineString(points) | |
poly = line.convex_hull | |
P = uniform_sample(poly, n=1000) | |
# Plot the polygon, the Delaunay triangulation and the points | |
# randomly sampled across its face. | |
shape = np.array(poly.exterior.coords) | |
plt.hold(True) | |
for p in triangulate(poly): | |
coords = np.array(p.exterior.coords) | |
plt.plot(coords[:,0], coords[:,1], 'g--') | |
plt.plot(shape[:,0], shape[:,1], 'k') | |
plt.plot(P[:,0], P[:,1], 'b.') |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment