2D Irregular Polygon to 3D Prism Region of Interest to Slice STL File

[nickvazz]

Description

Hello! Thank you for this great library!

A little background into my problem, I have a elevation map in 3D (STL):

that I am trying to cut into puzzle pieces (from the top down):

which I have represented as shapely.Polygons.

My thought was that I would be able to extend this 2D polygon into a 3D polygon and slice my elevation map into (in this case) 12 individual puzzle piece meshes using some combination of these examples in the docs:

• Clip with Rotated Box
• Drape 2D Surface From Line.
• Clipping with a Surface
• Boolean Union
• Extending a 2D StructuredGrid to 3D

Am I missing something for making an irregular prism along the z-axis with my irregular polygon as the top and bottom faces?

Example Data

The puzzle polygons could come from the same data structure seen in this sample plot:

import numpy as np
from shapely.geometry import Polygon

points = np.array([
    [0,0],
    [1,1],
    [1,0],
    [2,1],
    [2,-1],
    [0,0]
])

plt.plot(points[:,0], points[:,1])
plt.show()

poly = Polygon(points)
poly

Issue?

When I try to use a 2D shapely.Polygon:

As pyvista.PolyData, it is not connected and not a surface:

Looking at the connectivity

[akaszynski]

Thanks for posting such a thorough explanation!

The issue is that pyvista needs the connectivity array to know which points are connected and in what manner. Since the PolyData class can be used to represent lines or areas, we need to know how these points should be represented and connected, otherwise, they'll just be a set of points.

I adapted your example and set lines based on connectivity. I recommend using numpy to generate the connectivity using np.hstack and np.arange or using a python loop as it's not necessary to hard code it.

import numpy as np
import pyvista as pv

points = np.array([
    [0, 0, 0],
    [1, 1, 0],
    [1, 0, 0],
    [2, 1, 0],
    [2, -1, 0],
    [0, 0, 0]
])

lines = [2, 0, 1,
         2, 1, 2,
         2, 2, 3,
         2, 3, 4,
         2, 4, 5]

pd = pv.PolyData(points)
pd.lines = lines

pd.plot()

[nickvazz]

Thank you for your quick response!

This looks like I'll be able to make it work 😄

lines = [2, 0, 1,
         2, 1, 2,
         2, 2, 3,
         2, 3, 4,
         2, 4, 5]

I think I might be not understanding what the 2's in the first component of each element in the lines mean. Do you think you might be able to clarify what I should be trying to create for the lines?

[nickvazz]

I think I got it! but not sure what those 2's are for still 😛

[akaszynski]

The "2" tells the number of points per segment or face. If you choose a higher number and specify more points per line, you can create complex lines or if you specify faces then you can create polyfaces.

[nickvazz]

Ahh I see! Thank you again.

I was able to make faces look just about right but seems like there is an artifact or something included from (maybe?) the centroid of the volume / starting points for the face to seemingly random edges. Does this have something to do with the origin or something?

[akaszynski]

It looks like you're trying to make a face out of a the a boundary of points. Try taking a look at https://docs.pyvista.org/examples/00-load/create-tri-surface.html, using the delaunay_2d method. There's an option for specifying the bounds that should work for you.

[nickvazz]

Just wanted to close the loop on this, turns out the delaunay_2d would fill in space I had no intended.

My current solution is to use voxelization of the sub-mesh (puzzle piece) I create and find what voxels overlap with the original mesh.

import pyvista as pv

# generate a bunch of pieces (sub mesh - not required to be "water-tight")

vox = pv.voxelize(piece, density=3)
piece = vox.clip_surface(stl_poly_mesh).extract_surface()

Not sure if there is a better way 🤷🏾 but seems to work where speed of generation is heavily dependent on the voxelize's density. That said, it looks like having varied densities in each direction might be possible

def of voxelize

def voxelize(mesh, density=None, check_surface=True):
    """Voxelize mesh to UnstructuredGrid.
    Parameters
    ----------
    density : float
        The uniform size of the voxels. Defaults to 1/100th of the mesh length.
    check_surface : bool
        Specify whether to check the surface for closure. If on, then the
        algorithm first checks to see if the surface is closed and
        manifold. If the surface is not closed and manifold, a runtime
        error is raised.
    """
    if not pyvista.is_pyvista_dataset(mesh):
        mesh = pyvista.wrap(mesh)
    if density is None:
        density = mesh.length / 100
    x_min, x_max, y_min, y_max, z_min, z_max = mesh.bounds
    x = np.arange(x_min, x_max, density)
    y = np.arange(y_min, y_max, density)
    z = np.arange(z_min, z_max, density)
    x, y, z = np.meshgrid(x, y, z)

    # Create unstructured grid from the structured grid
    grid = pyvista.StructuredGrid(x, y, z)
    ugrid = pyvista.UnstructuredGrid(grid)

    # get part of the mesh within the mesh's bounding surface.
    selection = ugrid.select_enclosed_points(mesh.extract_surface(),
                                             tolerance=0.0,
                                             check_surface=check_surface)
    mask = selection.point_arrays['SelectedPoints'].view(np.bool_)

    # extract cells from point indices
    return ugrid.extract_points(mask)

which would potentially reduce unnecessary compute times (in my specific case).

Would you recommend that I open PR related to adding that variable density in or just reimplement the function on my end?

[akaszynski]

Would you recommend that I open PR related to adding that variable density in or just reimplement the function on my end?

A PR would be great. Generally voxelization is uniform, but having support for "stretched" voxels sounds like a feature worth having.

[nickvazz]

Sounds good! Opened that PR-1100