Stereographic Duoprism with Python (PyVista)

.gitignore

pic*

StereoDuoprism.py

# -*- 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()

StereoDuoprism_12-3.gif

StereoDuoprism_3-30.gif

StereoDuoprism_8-4.gif

StereoDuoprism_colors_rotation.py

StereoDuoprism_rotation.py