Polygons, triangulations, and the medial axis
from curvey import Curve, Polygon, Edges
from matplotlib import pyplot as plt
from matplotlib.patches import PathPatch
import numpy as np
plt.rcParams["figure.figsize"] = (3, 3)
plt.rcParams["axes.titlesize"] = 10
plt.rcParams["axes.labelsize"] = 8
plt.rcParams["xtick.labelsize"] = 8
plt.rcParams["ytick.labelsize"] = 8
Polygon construction¶
A Polygon is bounded by Curve.
exterior = Curve.circle(n=100)
poly0 = Polygon(exterior)
_ = poly0.plot_polygon()
Interior boundaries with opposite orientation to the exterior can also be specified:
interiors=(
Curve.circle(r=0.25, n=50).translate([-0.5, 0]).to_cw(),
Curve.circle(r=0.25, n=50).translate([0.5, 0]).to_cw(),
)
poly1 = Polygon(exterior, interiors)
_ = poly1.plot_polygon()
The signed_area property calculates area enclosed by the polygon:
poly0.signed_area, poly1.signed_area
(3.139525976465669, 2.747859621577218)
The combined external and internal boundary edges are available in an 'edge soup' representation as an curvey.Edges object.
poly1.boundary
Edges(n_points=200, n_edges=200)
For convenience, polygons can also be constructed from matplotlib's font rendering modules:
polys = Polygon.from_text("curvey", family="arial", size=72)
for p in polys:
p.plot_polygon()
The Polygon.apply method applies any Curve -> Curve function to its boundarys, returning a new Polygon:
eight = Polygon.from_text("8")[0]
eight = eight.apply(Curve.interpolate_thresh, thresh=5)
eight.plot_edges()
Triangulation¶
The Polygon.triangulate() returns a constrained Delaunay triangulation of the polygon with the bindings to Johnathan Shewchuck's triangle library.
tris0 = eight.triangulate()
_ = tris0.plot_tris()
Triangulate also supports area and angle constraints.
tris = eight.triangulate(max_tri_area=9)
_ = tris.plot_tris()
In order to motivate the use of the medial axis in the next section, we'll introduce here the problem of assigning a z-coordinate to each vertex in the triangulation in order to provide a 'roof' over the boundary of the polygon. A simple way to do this would be to simply set the height of the vertices equal to their distance from the boundary:
_edge_idx, boundary_dist = eight.boundary.closest_edge(tris.points)
fig, ax = plt.subplots(subplot_kw={'projection': '3d'}, figsize=(5, 5))
x, y = tris.points.T
ax.plot_trisurf(x, y, Z=boundary_dist, triangles=tris.faces)
ax.set_zlim([0, 15])
(0.0, 15.0)
The problem is that the edges of the triangulation fall erratically over the midline of the polygon.
The medial axis¶
The medial axis of a polygon is the set of points at the center of all maximally inscribed disks inside the polygon. Polygon.approximate_medial_axis constructs the approximate medial axis of the polygon as per
It returns an Edges object, with vertex distances from the polygon boundary stored in the point metadata property 'distance':
ama = eight.approximate_medial_axis(
dist_thresh=1,
abs_err=0.25,
)
ama
Edges(n_points=52, n_edges=53, point_data={distance})
_ = eight.plot(color='black')
_ = ama.plot_edges(color='grey', directed=False)
_ = ama.plot_points(color=ama.point_data['distance'])
_ = plt.axis('off')
For a nicer triangulation that follows the medial axis, we can combine the edges from the original polygon and the medial axis. We assign a point data 'distance' property of zero to the original edges so that all vertices in the combined edges have a distance property.
Here below we're calling Edges.triangulate instead of Polygon.triangulate so we need to be a little more explicit, including passing points interior to the polygon holes.
all_edges = Edges.concatenate(
eight.boundary.with_point_data(distance=0),
ama,
)
tris = all_edges.triangulate(
polygon=True,
max_tri_area=3,
holes=eight.hole_points,
)
Finally, instead of assigning the vertices a z-coordinate based on their exact distance from the boundary, we instead can interpolate distance between the boundary and the medial axis using scipy's radial basis function interpolator, which should result in a much smoother distance function.
from scipy.interpolate import RBFInterpolator
dist_interpolator = RBFInterpolator(
all_edges.points,
all_edges.point_data['distance'],
kernel='thin_plate_spline',
degree=1,
)
dist = dist_interpolator(tris.points)
dist[tris.is_boundary_vertex] = 0
Finally, instead of setting the z-coordinate directly to their distance, apply a sinusoidal profile:
fig, ax = plt.subplots(subplot_kw={'projection': '3d'}, figsize=(6, 6))
z = dist.max() * np.sin(dist / dist.max())
x, y = tris.points.T
_ = ax.plot_trisurf(x, y, z, triangles=tris.faces)
_ = ax.set_zlim([0, 15])
