1# https://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6from __future__ import absolute_import, print_function, division
8import math
9import numpy as num
10from .geometry import arr_vertices, arr_faces, normalize, refine_triangles, \
11 vdot, vnorm, face_centers
13r2d = 180./math.pi
16def cube():
17 vertices = arr_vertices([
18 [-1, -1, -1],
19 [-1, -1, 1],
20 [-1, 1, -1],
21 [-1, 1, 1],
22 [1, -1, -1],
23 [1, -1, 1],
24 [1, 1, -1],
25 [1, 1, 1]])
27 faces = arr_faces([
28 [0, 1, 3, 2],
29 [0, 4, 5, 1],
30 [4, 6, 7, 5],
31 [2, 3, 7, 6],
32 [1, 5, 7, 3],
33 [0, 2, 6, 4]])
35 return vertices, faces
38def triangles_to_center(vertices, faces):
40 vs = vertices
41 fs = faces
43 vcs = face_centers(vs, fs)
45 nv = vs.shape[0]
46 nf = fs.shape[0]
47 nc = fs.shape[1]
49 intc = (nv + num.arange(nf))[:, num.newaxis]
50 f2s = num.vstack([
51 num.hstack([fs[:, (ic, (ic+1) % nc)], intc])
52 for ic in range(nc)])
54 v2s = num.vstack([vs, vcs])
55 return v2s, f2s
58def tcube():
59 vs, fs = cube()
60 return triangles_to_center(vs, fs)
63def tetrahedron():
64 vertices = arr_vertices([
65 [math.sqrt(8./9.), 0., -1./3.],
66 [-math.sqrt(2./9.), math.sqrt(2./3.), -1./3.],
67 [-math.sqrt(2./9.), -math.sqrt(2./3.), -1./3.],
68 [0., 0., 1.]
69 ])
71 faces = arr_faces([
72 [2, 1, 0],
73 [3, 2, 0],
74 [2, 3, 1],
75 [3, 0, 1]
76 ])
77 return vertices, faces
80def icosahedron():
81 a = 0.5 * (math.sqrt(5) - 1.0)
83 vertices = arr_vertices([
84 [0, 1, a], [a, 0, 1], [1, a, 0],
85 [0, 1, -a], [-a, 0, 1], [1, -a, 0],
86 [0, -1, -a], [-a, 0, -1], [-1, -a, 0],
87 [0, -1, a], [a, 0, -1], [-1, a, 0]
88 ])
90 faces = arr_faces([
91 [6, 5, 9], [9, 8, 6], [8, 7, 6], [7, 10, 6], [10, 5, 6],
92 [5, 10, 2], [2, 1, 5], [5, 1, 9], [9, 1, 4], [4, 8, 9],
93 [4, 11, 8], [8, 11, 7], [7, 11, 3], [3, 10, 7], [3, 2, 10],
94 [3, 0, 2], [0, 1, 2], [0, 4, 1], [0, 11, 4], [0, 3, 11]
95 ])
97 return vertices, faces
100def neighbors(vertices, faces):
102 nv = vertices.shape[0]
103 nf, nc = faces.shape
104 fedges = num.zeros((nf*nc, 2), dtype=int)
105 for ic in range(nc):
106 fedges[ic::nc, :] = faces[:, (ic, (ic+1) % nc)]
108 indface = num.repeat(num.arange(nf), nc)
109 fedges1 = fedges[:, 0] + fedges[:, 1]*nv
110 sortorder = fedges1.argsort()
111 fedges1_rev = fedges[:, 1] + fedges[:, 0]*nv
113 inds = num.searchsorted(fedges1, fedges1_rev, sorter=sortorder)
114 # todo: handle cases without neighbors
115 assert num.all(fedges1[sortorder[inds]] == fedges1_rev)
116 neighbors = indface[sortorder[inds]].reshape((nf, nc))
117 return neighbors
120def adjacent_faces(vertices, faces):
122 nv = vertices.shape[0]
123 nf, nc = faces.shape
125 iverts = faces.reshape(nf*nc)
126 ifaces = num.repeat(num.arange(nf), nc)
127 iverts_order = iverts.argsort()
129 vs = vertices[faces]
130 gs = num.zeros(vs.shape)
131 gs[:, :-1, :] = vs[:, 1:, :] - vs[:, :-1, :]
132 gs[:, -1, :] = vs[:, 0, :] - vs[:, -1, :]
134 vecs = gs.reshape((nf*nc, 3))
135 iverts_ordered = iverts[iverts_order]
136 ifirsts = nextval_indices(iverts_ordered)
137 iselected = iverts_ordered[ifirsts]
139 plane_en = normalize(vertices)
140 plane_e1 = num.zeros((nv, 3))
141 plane_e1[iselected, :] = normalize(
142 project_to_plane_nn(plane_en[iselected, :], vecs[ifirsts, :]))
143 plane_e2 = num.zeros((nv, 3))
144 plane_e2[iselected, :] = num.cross(
145 plane_en[iselected, :], plane_e1[iselected, :])
147 a1 = vdot(vecs[iverts_order, :], plane_e1[iverts_ordered, :])
148 a2 = vdot(vecs[iverts_order, :], plane_e2[iverts_ordered, :])
149 angles = num.arctan2(a1, a2) * r2d + 360. * iverts_ordered
150 iverts_order2 = num.argsort(angles)
152 return iverts_ordered, ifaces[iverts_order][iverts_order2]
155def nextval_indices(a):
156 return num.concatenate(
157 [[0], num.where(num.diff(a) != 0)[0] + 1])
160def project_to_plane(vns, vas):
161 return vas - (vdot(vas, vns) / vdot(vns, vns))[:, num.newaxis] * vns
164def project_to_plane_nn(vns, vas):
165 return vas - (vdot(vas, vns))[:, num.newaxis] * vns
168def corner_handednesses(vertices, faces):
169 vs = vertices[faces]
170 gs = num.zeros(vs.shape)
171 gs[:, :-1, :] = vs[:, 1:, :] - vs[:, :-1, :]
172 gs[:, -1, :] = vs[:, 0, :] - vs[:, -1, :]
173 hs = num.zeros((faces.shape[0], faces.shape[1]))
174 hs[:, 1:] = vdot(num.cross(gs[:, :-1, :], gs[:, 1:, :]), vs[:, 1:, :])
175 hs[:, 0] = vdot(num.cross(gs[:, -1, :], gs[:, 0, :]), vs[:, 0, :])
176 return hs
179def truncate(vertices, faces):
181 nv = vertices.shape[0]
182 iverts, ifaces = adjacent_faces(vertices, faces)
183 ifirsts = nextval_indices(iverts)
184 ilengths = num.zeros(ifirsts.size, dtype=int)
185 ilengths[:-1] = num.diff(ifirsts)
186 ilengths[-1] = iverts.size - ifirsts[-1]
187 nc = num.max(ilengths)
188 nf = nv
189 vertices_new = face_centers(vertices, faces)
190 faces_new = num.zeros((nf, nc), dtype=int)
191 for iface in range(nf):
192 ifirst = ifirsts[iface]
193 ilength = ilengths[iface]
194 faces_new[iface, :ilength] = ifaces[ifirst:ifirst+ilength]
195 faces_new[iface, ilength:] = ifaces[ifirst+ilength-1]
197 faces_new = faces_new[:, ::-1]
199 return vertices_new, faces_new
202def iter_icospheres(order, inflate=True):
203 vertices, faces = icosahedron()
204 vertices /= vnorm(vertices)[:, num.newaxis]
205 yield (vertices, faces)
207 for i in range(order):
208 vertices, faces = refine_triangles(vertices, faces)
209 if inflate:
210 vertices = normalize(vertices)
212 yield (vertices, faces)
215bases = {
216 'icosahedron': icosahedron,
217 'tetrahedron': tetrahedron,
218 'tcube': tcube}
221def sphere(
222 order,
223 base=icosahedron,
224 kind='kind1',
225 inflate=True,
226 radius=1.0,
227 triangulate=True):
229 if isinstance(base, str):
230 base = bases[base]
232 vertices, faces = base()
233 vertices /= vnorm(vertices)[:, num.newaxis]
235 for i in range(order):
236 vertices, faces = refine_triangles(vertices, faces)
237 if inflate:
238 vertices = normalize(vertices)
239 if radius != 1.0:
240 vertices *= radius
242 if kind == 'kind2':
243 vertices, faces = truncate(vertices, faces)
244 if triangulate:
245 vertices, faces = triangles_to_center(vertices, faces)
247 return vertices, faces