Skip to content

Instantly share code, notes, and snippets.

@stla
Last active Oct 23, 2021
Embed
What would you like to do?
Stereographic Duoprism with Python (PyVista)
# -*- coding: utf-8 -*-
from math import sqrt, sin, cos, pi
import pyvista as pv
import numpy as np
A = 8 # number of sides of the first polygon
B = 4 # number of sides of the second polygon
# construction of the vertices
vertices = np.empty((A, B, 4))
for i in range(A):
v1 = [cos(i/A*2*pi), sin(i/A*2*pi)]
for j in range(B):
v2 = [cos(j/B*2*pi), sin(j/B*2*pi)]
vertices[i, j, :] = np.array(v1 + v2)
# construction of the edges
edges = np.empty((2, 2, 2*A*B), dtype=int)
def dominates(c1, c2):
return c2[0]>c1[0] or (c2[0]==c1[0] and c2[1]>c1[1])
counter = 0
for i in range(A):
for j in range(B):
c1 = np.array([i, j])
candidate = np.array([i, (j-1) % B])
if dominates(c1, candidate):
edges[:, :, counter] = np.column_stack((c1, candidate))
counter += 1
candidate = np.array([i, (j+1) % B])
if dominates(c1, candidate):
edges[:, :, counter] = np.column_stack((c1, candidate))
counter += 1
candidate = np.array([(i-1) % A, j])
if dominates(c1, candidate):
edges[:, :, counter] = np.column_stack((c1, candidate))
counter += 1
candidate = np.array([(i+1) % A, j])
if dominates(c1, candidate):
edges[:, :, counter] = np.column_stack((c1, candidate))
counter += 1
# stereographic projection
def stereog(v):
return v[0:3] / (sqrt(2) - v[3])
# spherical segment
def sphericalSegment(P, Q, n):
out = np.empty((n+1, 4))
for i in range(n+1):
pt = P + (i/n)*(Q-P)
out[i, :] = sqrt(2)/np.linalg.norm(pt) * pt
return out
# stereographic edge
def polyline_from_points(points):
poly = pv.PolyData()
poly.points = points
the_cell = np.arange(0, len(points), dtype=np.int_)
the_cell = np.insert(the_cell, 0, len(points))
poly.lines = the_cell
return poly
def stereoEdge(verts, v1, v2):
P = verts[v1[0], v1[1], :]
Q = verts[v2[0], v2[1], :]
PQ = sphericalSegment(P, Q, 100)
pq = np.apply_along_axis(stereog, 1, PQ)
dists = np.sqrt(np.apply_along_axis(lambda x: np.vdot(x, x), 1, pq))
radii = dists / 15
r = radii.min()
rfactor = radii.max() / r
polyline = polyline_from_points(pq)
polyline["R"] = radii
tube = polyline.tube(
radius=r, scalars="R", radius_factor=rfactor, n_sides=200
)
return tube
# projected vertices
vs = np.apply_along_axis(stereog, 2, vertices)
####~~~~ plot ~~~~####
pltr = pv.Plotter(window_size=[512, 512])
pltr.set_background("#363940")
## plot the edges
for k in range(2*A*B):
v1 = edges[:, 0, k]
v2 = edges[:, 1, k]
edge = stereoEdge(vertices, v1, v2)
pltr.add_mesh(edge, color = "gold", specular=10, smooth_shading=True)
## plot the vertices
for i in range(A):
for j in range(B):
v = vs[i, j, :]
ball = pv.Sphere(np.linalg.norm(v)/10, center = v)
pltr.add_mesh(ball, color="#EEC900", specular=10, smooth_shading=True)
pltr.show()
View raw

(Sorry about that, but we can’t show files that are this big right now.)

View raw

(Sorry about that, but we can’t show files that are this big right now.)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment