1

I wish to take a concave and complex (containing holes) polygon and extrude it 'vertically' into a polyhedron, purely for visualisation. I begin with a shapely Polygon, like below:

poly = Polygon(
    [(0,0), (10,0), (10,10), (5,8), (0,10), (1,7), (0,5), (1,3)], 
    holes=[
        [(2,2),(4,2),(4,4),(2,4)],
        [(6,6), (7,6), (6.5,6.5), (7,7), (6,7), (6.2,6.5)]])

which I correctly plot (reorientating the exterior coordinates to be clockwise, and the hole coordinates to be counterclockwise) in matplotlib as:

enter image description here

I then seek to render this polygon extruded out-of-the-page (along z), using PyVista. There are a few hurdles; PyVista doesn't directly support concave (nor complex) input to its PolyData type. So we first create an extrusion of simple (hole-free) concave polygons, as per this discussion.

def extrude_simple_polygon(xy, z0, z1):

    # force counter-clockwise ordering, so PyVista interprets polygon correctly
    xy = _reorient_coords(xy, clockwise=False)

    # remove duplication of first & last vertex
    xyz0 = [(x,y,z0) for x,y in xy]
    if (xyz0[0] == xyz0[-1]):
        xyz0.pop()

    # explicitly set edge_source
    base_vert = [len(xyz0)] + list(range(len(xyz0)))
    base_data = pyvista.PolyData(xyz0, base_vert)
    base_mesh = base_data.delaunay_2d(edge_source=base_data)
    vol_mesh  = base_mesh.extrude((0, 0, z1-z0), capping=True)

    # force triangulation, so PyVista allows boolean_difference
    return vol_mesh.triangulate()

Observe this works when extruding the outer polygon and each of its internal polygons in-turn:

extrude_simple_polygon(list(poly.exterior.coords), 0, 5).plot()

enter image description here

extrude_simple_polygon(list(poly.interiors[0].coords), 0, 5).plot()
extrude_simple_polygon(list(poly.interiors[1].coords), 0, 5).plot()

enter image description here enter image description here

I reasoned that to create an extrusion of the original complex polygon, I could compute the boolean_difference. Alas, the result of

outer_vol = extrude_simple_polygon(list(poly.exterior.coords), 0, 5)
for hole in poly.interiors:
    hole_vol = extrude_simple_polygon(list(hole.coords), 0, 5)
    outer_vol = outer_vol.boolean_difference(hole_vol)

outer_vol.plot()

is erroneous:

enter image description here

The doc advises to inspect the normals via plot_normals, revealing that all extruded volumes have inward-pointing (or else, unexpected) normals:

enter image description here enter image description here enter image description here

The extrude doc mentions nothing of the extruded surface normals nor the original object (in this case, a polygon) orientation.

We could be forgiven for expecting our polygons must be clockwise, so we set clockwise=True in the first line of extrude_simple_polygon and try again. Alas, PolyData now misinterprets our base polygon; calling base_mesh.plot() reveals (what should look like our original blue outer polygon):

enter image description here

with extrusion

enter image description here

  • Does PyVista always expect counter-clockwise polygons?
  • Why does extrude create volumes with inward-pointing surface normals?
  • How can I correct the extruded surface normals?
  • Otherwise, how can I make PyVista correctly visualise what should be an incredibly simply-extruded concave complex polygon??
Anti Earth
  • 4,671
  • 13
  • 52
  • 83

1 Answers1

4

You're very close. What you have to do is use a single call to delaunay_2d() with all three polygons (i.e., the enclosing one and the two holes) as edge source (loop source?). It's also important to have faces (rather than lines) from each polygon; this is what makes it possible to enforce the holeyness of the holes.

Here's a complete example for your input (where I manually flipped the orientation of the holes; you seem to have a _reorient_coords() helper that you should use instead):

import pyvista as pv

# coordinates of enclosing polygon
poly_points = [
    (0, 0), (10, 0), (10, 10), (5, 8), (0, 10), (1, 7), (0, 5), (1, 3),
]
# hole point order hard-coded here; use your _reorient_coords() function
holes = [
    [(2, 2), (4, 2), (4, 4), (2, 4)][::-1],
    [(6, 6), (7, 6), (6.5, 6.5), (7, 7), (6, 7), (6.2, 6.5)][::-1],
]

z0, z1 = 0.0, 5.0

def points_2d_to_poly(points, z):
    """Convert a sequence of 2d coordinates to a polydata with a polygon."""
    faces = [len(points), *range(len(points))]
    poly = pv.PolyData([p + (z,) for p in points], faces=faces)
    return poly

# bounding polygon
polygon = points_2d_to_poly(poly_points, z0)

# add all holes
for hole_points in holes:
    polygon += points_2d_to_poly(hole_points, z0)

# triangulate poly with all three subpolygons supplying edges
# (relative face orientation is critical here)
polygon_with_holes = polygon.delaunay_2d(edge_source=polygon)

# extrude
holey_solid = polygon_with_holes.extrude((0, 0, z1 - z0), capping=True)
holey_solid.plot()

screenshot of perspective view of the expected extruded "volume"

Here's the top view of the polygon pre-extrusion:

plotter = pv.Plotter()
plotter.add_mesh(polygon_with_holes, show_edges=True, color='cyan')
plotter.view_xy()
plotter.show()

top view of polygon showing holes where expected, and the edges of the triangulation elsewhere

  • Any idea why this breaks for the below simple, convex (crooked cross-shaped) polygon? It seems to discard a vertex, interpolating an edge between to non-adjacent vertices. `poly_points = [(227.0, 119.0), (216.910083959346, 343.61030490499337), (349.0, 348.0), (333.0, 444.0), (212.50063243106501, 441.76853023020493), (204.0, 631.0), (93.0, 632.0), (95.0, 625.0), (100.03766532489936, 439.68588269120187), (9.0, 438.0), (18.0, 337.0), (102.75255979697455, 339.8165503255792), (109.0, 110.0)]` – Anti Earth Jan 23 '23 at 23:41
  • I've resolved the above bug (the superfluous specification of the final vertex as the first, + needing to ensure the coordinate list does not have its own superfluous final vertex repeat), and edited it into your answer – Anti Earth Jan 24 '23 at 17:24
  • 1
    @AntiEarth sorry, I couldn't check this on a laptop sooner. Thanks for the bugfix, the repetition of the first point was because I originally defined lines rather than faces, and you need to come back to the first point if you want a cosed line. But since we need _faces_ this was a mistake. I can confirm that with the updated `faces` definition your cross example works. – Andras Deak -- Слава Україні Jan 24 '23 at 23:30
  • This solution fails for `poly_points = [(389.0, 78.0, 230), (389.0, 784.0, 230), (-1021.0, 784.0, 230), (-1021.0, 78.0, 230)]`, and `hole1_points = [(-396.0, 190.0, 230), (-396.0, 350.0, 230), (-350.0, 350.0, 230), (-350.0, 649.0, 230), (-201.0, 649.0, 230), (-201.0, 350.0, 230), (-152.0, 350.0, 230), (-152.0, 190.0, 230)]`, and `hole2_points = [(-653.0, 201.0, 230), (-653.0, 649.0, 230), (-499.0, 649.0, 230), (-499.0, 201.0, 230)]`. Without extrusion, [this](https://pasteboard.co/nr3sVzUvVFsm.png) is erroneously rendered as [this](https://pasteboard.co/SdTCYjLCzfsW.png) Any idea why? – Anti Earth Feb 21 '23 at 14:52
  • @AntiEarth sorry for the late reply, life's been very busy. If I shift the first hole with +100 along the x axis it starts to work. It looks like the triangulation is ill-conditioned there (there's also a VTK warning in this case). I could fix the problem by making the hole definitions have denser points along the line. One way to do this is with `import numpy as np`, then `scale = np.linspace(0, 1, 2, endpoint=False)[:, None]; holes = [[p for p1, p2 in zip(hole, hole[1:]) for p in p2*scale + p1*(1 - scale)] + [hole[-1]] for hole in holes]`. This doubles the density. – Andras Deak -- Слава Україні Mar 05 '23 at 14:26
  • Using `np.linspace(0, 1, 10)` would give you 10x as many points along each hole contour. To be clear, in the previous comment `holes = [hole1_points, hole2_points]` according to your comment. – Andras Deak -- Слава Україні Mar 05 '23 at 14:27
  • No need to apologise - you've been a fantastic help! Manually refining the mesh worries me, since I use PyVista for 3D visualisation within a user-facing Python library, and need something completely robust. My polyhedra (as described above) are vertically-extruded polygons specified as vertices, and their base polygons render fine in matplotlib. Is PyVista the wrong tool to render their 3D extrusions? In fact, my library creates meticulous conformal meshes and PyVista is used to visualise the "base symbolic scene", *not* the (numerically damaged) mesh. Are there non-mesh-based renderers? – Anti Earth Mar 06 '23 at 15:30
  • 1
    @AntiEarth the extrusions would be fine, where this breaks (and it's VTK doing the breaking, not PyVista) is the triangulation of the 2d base before extrusion. Once you have a proper polygon I'd expect extrusion to Just Work^TM. But in any case I don't really have experience in this regard, and in the broader field. E.g. I'm not sure what "non-mesh-based renderer" might mean. You should consider [asking on the PyVista github discussions page](https://github.com/pyvista/pyvista/discussions) where more eyes and of much broader backgrounds can chime in. – Andras Deak -- Слава Україні Mar 06 '23 at 15:52
  • 1
    Seems so - I've created an [Issue](https://github.com/pyvista/pyvista/issues/4097), and will update this question with any progress. By "non-mesh-based renderer", I mean any facility which lets me specify complex polygons by their coordinates and render them faithfully, without having to perform explicit meshing or tweaking, nor project the coordinates into a (more) approximate subspace. – Anti Earth Mar 07 '23 at 14:37