Instantly share code, notes, and snippets.

# stla/.gitignore

Last active Oct 23, 2021
Stereographic Duoprism with Python (PyVista)
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 pic*
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 # -*- 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.)