1# https://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import math
7import numpy as num
8from .geometry import arr_vertices, arr_faces, normalize, refine_triangles, \
9 vdot, vnorm, face_centers
11r2d = 180./math.pi
14def cube():
15 vertices = arr_vertices([
16 [-1, -1, -1],
17 [-1, -1, 1],
18 [-1, 1, -1],
19 [-1, 1, 1],
20 [1, -1, -1],
21 [1, -1, 1],
22 [1, 1, -1],
23 [1, 1, 1]])
25 faces = arr_faces([
26 [0, 1, 3, 2],
27 [0, 4, 5, 1],
28 [4, 6, 7, 5],
29 [2, 3, 7, 6],
30 [1, 5, 7, 3],
31 [0, 2, 6, 4]])
33 return vertices, faces
36def triangles_to_center(vertices, faces):
38 vs = vertices
39 fs = faces
41 vcs = face_centers(vs, fs)
43 nv = vs.shape[0]
44 nf = fs.shape[0]
45 nc = fs.shape[1]
47 intc = (nv + num.arange(nf))[:, num.newaxis]
48 f2s = num.vstack([
49 num.hstack([fs[:, (ic, (ic+1) % nc)], intc])
50 for ic in range(nc)])
52 v2s = num.vstack([vs, vcs])
53 return v2s, f2s
56def tcube():
57 vs, fs = cube()
58 return triangles_to_center(vs, fs)
61def tetrahedron():
62 vertices = arr_vertices([
63 [math.sqrt(8./9.), 0., -1./3.],
64 [-math.sqrt(2./9.), math.sqrt(2./3.), -1./3.],
65 [-math.sqrt(2./9.), -math.sqrt(2./3.), -1./3.],
66 [0., 0., 1.]
67 ])
69 faces = arr_faces([
70 [2, 1, 0],
71 [3, 2, 0],
72 [2, 3, 1],
73 [3, 0, 1]
74 ])
75 return vertices, faces
78def icosahedron():
79 a = 0.5 * (math.sqrt(5) - 1.0)
81 vertices = arr_vertices([
82 [0, 1, a], [a, 0, 1], [1, a, 0],
83 [0, 1, -a], [-a, 0, 1], [1, -a, 0],
84 [0, -1, -a], [-a, 0, -1], [-1, -a, 0],
85 [0, -1, a], [a, 0, -1], [-1, a, 0]
86 ])
88 faces = arr_faces([
89 [6, 5, 9], [9, 8, 6], [8, 7, 6], [7, 10, 6], [10, 5, 6],
90 [5, 10, 2], [2, 1, 5], [5, 1, 9], [9, 1, 4], [4, 8, 9],
91 [4, 11, 8], [8, 11, 7], [7, 11, 3], [3, 10, 7], [3, 2, 10],
92 [3, 0, 2], [0, 1, 2], [0, 4, 1], [0, 11, 4], [0, 3, 11]
93 ])
95 return vertices, faces
98def neighbors(vertices, faces):
100 nv = vertices.shape[0]
101 nf, nc = faces.shape
102 fedges = num.zeros((nf*nc, 2), dtype=int)
103 for ic in range(nc):
104 fedges[ic::nc, :] = faces[:, (ic, (ic+1) % nc)]
106 indface = num.repeat(num.arange(nf), nc)
107 fedges1 = fedges[:, 0] + fedges[:, 1]*nv
108 sortorder = fedges1.argsort()
109 fedges1_rev = fedges[:, 1] + fedges[:, 0]*nv
111 inds = num.searchsorted(fedges1, fedges1_rev, sorter=sortorder)
112 # todo: handle cases without neighbors
113 assert num.all(fedges1[sortorder[inds]] == fedges1_rev)
114 neighbors = indface[sortorder[inds]].reshape((nf, nc))
115 return neighbors
118def adjacent_faces(vertices, faces):
120 nv = vertices.shape[0]
121 nf, nc = faces.shape
123 iverts = faces.reshape(nf*nc)
124 ifaces = num.repeat(num.arange(nf), nc)
125 iverts_order = iverts.argsort()
127 vs = vertices[faces]
128 gs = num.zeros(vs.shape)
129 gs[:, :-1, :] = vs[:, 1:, :] - vs[:, :-1, :]
130 gs[:, -1, :] = vs[:, 0, :] - vs[:, -1, :]
132 vecs = gs.reshape((nf*nc, 3))
133 iverts_ordered = iverts[iverts_order]
134 ifirsts = nextval_indices(iverts_ordered)
135 iselected = iverts_ordered[ifirsts]
137 plane_en = normalize(vertices)
138 plane_e1 = num.zeros((nv, 3))
139 plane_e1[iselected, :] = normalize(
140 project_to_plane_nn(plane_en[iselected, :], vecs[ifirsts, :]))
141 plane_e2 = num.zeros((nv, 3))
142 plane_e2[iselected, :] = num.cross(
143 plane_en[iselected, :], plane_e1[iselected, :])
145 a1 = vdot(vecs[iverts_order, :], plane_e1[iverts_ordered, :])
146 a2 = vdot(vecs[iverts_order, :], plane_e2[iverts_ordered, :])
147 angles = num.arctan2(a1, a2) * r2d + 360. * iverts_ordered
148 iverts_order2 = num.argsort(angles)
150 return iverts_ordered, ifaces[iverts_order][iverts_order2]
153def nextval_indices(a):
154 return num.concatenate(
155 [[0], num.where(num.diff(a) != 0)[0] + 1])
158def project_to_plane(vns, vas):
159 return vas - (vdot(vas, vns) / vdot(vns, vns))[:, num.newaxis] * vns
162def project_to_plane_nn(vns, vas):
163 return vas - (vdot(vas, vns))[:, num.newaxis] * vns
166def corner_handednesses(vertices, faces):
167 vs = vertices[faces]
168 gs = num.zeros(vs.shape)
169 gs[:, :-1, :] = vs[:, 1:, :] - vs[:, :-1, :]
170 gs[:, -1, :] = vs[:, 0, :] - vs[:, -1, :]
171 hs = num.zeros((faces.shape[0], faces.shape[1]))
172 hs[:, 1:] = vdot(num.cross(gs[:, :-1, :], gs[:, 1:, :]), vs[:, 1:, :])
173 hs[:, 0] = vdot(num.cross(gs[:, -1, :], gs[:, 0, :]), vs[:, 0, :])
174 return hs
177def truncate(vertices, faces):
179 nv = vertices.shape[0]
180 iverts, ifaces = adjacent_faces(vertices, faces)
181 ifirsts = nextval_indices(iverts)
182 ilengths = num.zeros(ifirsts.size, dtype=int)
183 ilengths[:-1] = num.diff(ifirsts)
184 ilengths[-1] = iverts.size - ifirsts[-1]
185 nc = num.max(ilengths)
186 nf = nv
187 vertices_new = face_centers(vertices, faces)
188 faces_new = num.zeros((nf, nc), dtype=int)
189 for iface in range(nf):
190 ifirst = ifirsts[iface]
191 ilength = ilengths[iface]
192 faces_new[iface, :ilength] = ifaces[ifirst:ifirst+ilength]
193 faces_new[iface, ilength:] = ifaces[ifirst+ilength-1]
195 faces_new = faces_new[:, ::-1]
197 return vertices_new, faces_new
200def iter_icospheres(order, inflate=True):
201 vertices, faces = icosahedron()
202 vertices /= vnorm(vertices)[:, num.newaxis]
203 yield (vertices, faces)
205 for i in range(order):
206 vertices, faces = refine_triangles(vertices, faces)
207 if inflate:
208 vertices = normalize(vertices)
210 yield (vertices, faces)
213bases = {
214 'icosahedron': icosahedron,
215 'tetrahedron': tetrahedron,
216 'tcube': tcube}
219def sphere(
220 order,
221 base=icosahedron,
222 kind='kind1',
223 inflate=True,
224 radius=1.0,
225 triangulate=True):
227 if isinstance(base, str):
228 base = bases[base]
230 vertices, faces = base()
231 vertices /= vnorm(vertices)[:, num.newaxis]
233 for i in range(order):
234 vertices, faces = refine_triangles(vertices, faces)
235 if inflate:
236 vertices = normalize(vertices)
237 if radius != 1.0:
238 vertices *= radius
240 if kind == 'kind2':
241 vertices, faces = truncate(vertices, faces)
242 if triangulate:
243 vertices, faces = triangles_to_center(vertices, faces)
245 return vertices, faces