Plotting a polygon in 3D

[sunilpai]

Hello, I am interested in plotting the area inside of polygons in a 3D Surface plot in Python. I was able to successfully extrude the polygon boundary, but now I need to cap the prism at the top and bottom. My question is whether most viable strategy for this would be to use the surface with opacity feature that is currently being developed; in this case I would define a mask of points inside the polygon (so this would basically be an image where I only plot pixels inside the polygon). Or is there a more straightforward feature that does this already (e.g., using plotly shapes)? I tried using Mesh3d with the vertices of the polygon, but the triangulation doesn’t do a great job at showing the polygon:

[Emmanuelle]

Hi @sunilpai Mesh3d seems indeed a good way to go. Could you please share a standalone script (maybe using dummy data) so that we can see what’s wrong or help you improve the plot?

[sunilpai]

I ran this code in a JupyterLab session (requires pip install gdspy). My main issue comes up when dealing with bezier curves which are typically used in CAD models:

import gdspy as gy
import numpy as np
import plotly.graph_objects as go

x = gy.Path(1)
poly = x.bezier([(5, 0), (5, 3), (10, 3)]).polygons[0].T
mesh = go.Mesh3d(x=poly[0], y=poly[1], z=0.1 * np.ones_like(poly[0]), color='red', opacity=0.5, showscale=False)
boundary = go.Scatter3d(x=poly[0], y=poly[1], z=np.ones_like(poly[0]) * 0.1, line = dict(color='red', width=2), mode='lines', showlegend=False)

fig = go.FigureWidget(data=[mesh, boundary])
fig

The above results in:

It seems the convex hull is being triangulated rather than the actual polygon? Not sure…

Would this Mesh3d approach also work for a ring (area between two concentric circles)?

[empet]

Hi @sunilpai,

To plot the polygonal surface bounded by two parallel Bezier curves as opposite sides and two segments on other sides we have to define a parameterization of that surface.

First let us define the trace for the boundary curve:

import gdspy as gy
import numpy as np
import plotly.graph_objects as go

p = gy.Path(1)
poly = p.bezier([(5, 0), (5, 3), (10, 3)]).polygons[0]
poly_boundary = np.concatenate((poly, [poly[0]]), axis=0)
x, y = poly_boundary.T

b_trace = go.Scatter3d(x=x, y=y, z=0.1*np.ones(x.shape), mode='lines', line_width=2, line_color='black')
bfig = go.Figure(data=[b_trace])
bfig.show()

Theoretically if r1(u) is a parameterization of the first Bezier curve , and r2(u) is the parametrerization of the second Bezier curve, with u in [0,1],
then the polygonal surface bounded by the two curves on two opposite sides, and two line segments on the other two sides,
has the parameterization:

s(u, v) = (1-v)*r1(u)+ v* r2(u) , with u, v in [0,1]

i.e. i.e. for u* fixed in [0,1], s(u*, v) is a segment connecting two corresponding points on the two Bezier curves.

The points on the first Bezier curve have as coordinates the rows in poly[:n],
while the corresponding points on the second curve have as coordinates the rows in poly[n:][::-1] (n is half of the number of points stored in the array poly)

The computational version of this theoretical approach is as follows:

N = poly.shape[0]
if N%2:
    raise ValueError("The polygon should have an even number of points")
n = N//2
m=10    
v = np.linspace(0, 1, m)
#curve1 = poly[:n]
curve2 = poly[n:][::-1]
#surface points (xs, ys, zs):
xs = np.outer(1-v, poly[:n,0])+np.outer(v, curve2[:, 0])
ys = np.outer(1-v, poly[:n,1])+np.outer(v, curve2[:, 1])
zs = 0.1*np.ones(xs.shape)   

surf = go.Surface(x=xs, y=ys, z=zs, colorscale= [[0, 'rgb(255,0,0)'], [1, 'rgb(255, 0, 0)']], showscale=False)
fig = go.Figure(data=[surf, b_trace])

[sunilpai]

Hi @empet, thanks for solving this particular case! My ultimate goal, though, is a bit more ambitious: I’d like to achieve this with any polygon as well as be able to do cases where there is a polygon with holes inside. What would you recommend for such cases? It seems here that it is required to be able to split the polygon into two curves that can be interpolated between (to parametrize a surface). However, not all polygons are like this! I was hoping something like Mesh3d would allow me to generalize beyond the bezier case.

[empet]

@sunilpai, I’m not familiar with gdspy. Today I installed it to get your polygon. All your polygons are generated by this library?

By polygon do you mean that you have its sides represented by arbitrary curves?

If the answer to the latter question is positive, then the corresponding surface can be defined as a Coons patch that has the parameterization as in (3) in these slides https://hci.uni-kl.de/fileadmin/user_upload/alggeom_script_ws17_07.pdf.
If you define the u, w meshgrids as masked arrays you can get the polygonal surfaces with holes . Please post here or elsewhere such data and I’ll try to code this method.

[sunilpai]

Hi, we should expect many of our polygons of interest to be simple polygons and bends as in the above example (gdspy). But some methods require arbitrary planar contours and the area in between (not straightforward to generate using gdspy or primitive shapes/curves). These contours will usually be smooth, but I’d ideally like to be able to handle polygons with sharp edges too (essentially, arbitrary!).

For this test, I created a random-ish spiky polygon generator using numpy, and would ideally like to get the area outside of the inner polygon “holes” but inside the big polygon if possible:

import numpy as np
import plotly.graph_objects as go

np.random.seed(100)
def random_poly(avg_r, regularity, spike_amp, spike_freq, num_v, center=(0, 0), ecc=1):
    # generate n angle steps
    angle_steps = np.sort(np.random.dirichlet(regularity * np.ones(num_v + 1)) * 2 * np.pi)
    angles = np.cumsum(angle_steps) + np.random.rand() * 2 * np.pi
    # now generate the points
    rs = spike_amp * np.sin(angles * spike_freq + np.random.rand() * 2 * np.pi) ** 2 + avg_r
    return np.asarray([rs * np.cos(angles), ecc * rs * np.sin(angles)]) + np.asarray(center)[:, np.newaxis]

poly1 = random_poly(100, 0.05, 30, 3, 300, ecc=1, center=(-100, 0))
poly2 = random_poly(100, 0.05, 30, 3, 300, ecc=0.5, center=(100, 200))
bigpoly = random_poly(500, 0.05, 100, 3, 300, ecc=1)

bfig = go.FigureWidget(data=[go.Scatter3d(
    x=bigpoly[0], y=bigpoly[1],
    z=np.ones_like(bigpoly[0]), line = dict(color='black', width=2),
    mode='lines', showlegend=False),
 go.Scatter3d(
    x=poly1[0], y=poly1[1],
    z=np.ones_like(poly1[0]), line = dict(color='black', width=2),
    mode='lines', showlegend=False),
 go.Scatter3d(
    x=poly2[0], y=poly2[1],
    z=np.ones_like(poly2[0]), line = dict(color='black', width=2),
    mode='lines', showlegend=False),
])
bfig

[empet]

@sunilpai Finally I found a simpler solution:

In the 3d space a planar polygon with holes is defined as a go.Surface instance.

First we have to identify the biggest polygon, bigpolygon, i.e. the outer boundary of our surface, and the minimal rectangle with sides parallel with axes, that contains it.

Then we define a dense meshgrid on this rectangle, and for each pair, (bigpolygon, hole_polygon), we identify the meshgrid positions that are outside the bigpolygon , or inside the hole_polygon.

For, we are calling the function skimage.measure.points_in_poly() that returns a mask.

The final mask defined as a chained np.bitwise_and() of the intermediary masks helps us to remove the meshgrid points that are outside the bigpolygon or inside the polygons defining the holes.

import numpy as np
import plotly.graph_objects as go
from skimage.measure import  points_in_poly

np.random.seed(100)
def random_poly(avg_r, regularity, spike_amp, spike_freq, num_v, center=(0, 0), ecc=1):
    # generate n angle steps
    angle_steps = np.sort(np.random.dirichlet(regularity * np.ones(num_v + 1)) * 2 * np.pi)
    angles = np.cumsum(angle_steps) + np.random.rand() * 2 * np.pi
    # now generate the points
    rs = spike_amp * np.sin(angles * spike_freq + np.random.rand() * 2 * np.pi) ** 2 + avg_r
    return np.asarray([rs * np.cos(angles), ecc * rs * np.sin(angles)]) + np.asarray(center)[:, np.newaxis]

holepoly1 = random_poly(100, 0.05, 30, 3, 300, ecc=1, center=(-100, 0))
holepoly2 = random_poly(100, 0.05, 30, 3, 300, ecc=0.5, center=(100, 200))
bigpoly = random_poly(500, 0.05, 100, 3, 300, ecc=1)

# points_in_poly() works with arrays of shape (N,2). Hence we are working with the transpose of the above defined e polygon arrays
polygon1 = np.concatenate((bigpoly.T, holepoly1.T), axis=0)
polygon2 = np.concatenate((bigpoly.T, holepoly2.T), axis=0)
xmin, xmax = bigpoly.T[:, 0].min(), bigpoly.T[:, 0].max()
ymin, ymax = bigpoly.T[:, 1].min(),  bigpoly.T[:, 0].max()

n = 100  #choose m, n such that to get a sufficiently narrow grid/meshgrid
m = 200
X = np.linspace(xmin, xmax, n)
Y = np.linspace(ymin, ymax, m)
x, y = np.meshgrid(X,Y)
pts = np.dstack((x,y))
m, n, _  = pts.shape

pts = pts.reshape((m*n, 2))
mask = np.bitwise_and(points_in_poly(pts, polygon1),  points_in_poly(pts, polygon2))
mask = mask.reshape((m,n))

z = np.ones(mask.shape)
z[np.where(~mask)] = np.nan  # to avoid importing numpy.ma

surf = go.Surface(x= X, y=Y, z=z, colorscale=[[0, 'rgb(255,0,0)'], [1, 'rgb(255,0,0)']], showscale=False)
fig = go.Figure(data=[surf])