Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/icosphere.py: 29%

1# https://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

4# ---|P------/S----------~Lg---------- 

5 

6''' 

7Icosphere geometries - polyhedra approximating a sphere. 

8''' 

9 

10import math 

11import numpy as num 

12from .geometry import arr_vertices, arr_faces, normalize, refine_triangles, \ 

13 vdot, vnorm, face_centers 

14 

15r2d = 180./math.pi 

16 

17 

18def cube(): 

19 vertices = arr_vertices([ 

20 [-1, -1, -1], 

21 [-1, -1, 1], 

22 [-1, 1, -1], 

23 [-1, 1, 1], 

24 [1, -1, -1], 

25 [1, -1, 1], 

26 [1, 1, -1], 

27 [1, 1, 1]]) 

28 

29 faces = arr_faces([ 

30 [0, 1, 3, 2], 

31 [0, 4, 5, 1], 

32 [4, 6, 7, 5], 

33 [2, 3, 7, 6], 

34 [1, 5, 7, 3], 

35 [0, 2, 6, 4]]) 

36 

37 return vertices, faces 

38 

39 

40def triangles_to_center(vertices, faces): 

41 

42 vs = vertices 

43 fs = faces 

44 

45 vcs = face_centers(vs, fs) 

46 

47 nv = vs.shape[0] 

48 nf = fs.shape[0] 

49 nc = fs.shape[1] 

50 

51 intc = (nv + num.arange(nf))[:, num.newaxis] 

52 f2s = num.vstack([ 

53 num.hstack([fs[:, (ic, (ic+1) % nc)], intc]) 

54 for ic in range(nc)]) 

55 

56 v2s = num.vstack([vs, vcs]) 

57 return v2s, f2s 

58 

59 

60def tcube(): 

61 vs, fs = cube() 

62 return triangles_to_center(vs, fs) 

63 

64 

65def tetrahedron(): 

66 vertices = arr_vertices([ 

67 [math.sqrt(8./9.), 0., -1./3.], 

68 [-math.sqrt(2./9.), math.sqrt(2./3.), -1./3.], 

69 [-math.sqrt(2./9.), -math.sqrt(2./3.), -1./3.], 

70 [0., 0., 1.] 

71 ]) 

72 

73 faces = arr_faces([ 

74 [2, 1, 0], 

75 [3, 2, 0], 

76 [2, 3, 1], 

77 [3, 0, 1] 

78 ]) 

79 return vertices, faces 

80 

81 

82def icosahedron(): 

83 a = 0.5 * (math.sqrt(5) - 1.0) 

84 

85 vertices = arr_vertices([ 

86 [0, 1, a], [a, 0, 1], [1, a, 0], 

87 [0, 1, -a], [-a, 0, 1], [1, -a, 0], 

88 [0, -1, -a], [-a, 0, -1], [-1, -a, 0], 

89 [0, -1, a], [a, 0, -1], [-1, a, 0] 

90 ]) 

91 

92 faces = arr_faces([ 

93 [6, 5, 9], [9, 8, 6], [8, 7, 6], [7, 10, 6], [10, 5, 6], 

94 [5, 10, 2], [2, 1, 5], [5, 1, 9], [9, 1, 4], [4, 8, 9], 

95 [4, 11, 8], [8, 11, 7], [7, 11, 3], [3, 10, 7], [3, 2, 10], 

96 [3, 0, 2], [0, 1, 2], [0, 4, 1], [0, 11, 4], [0, 3, 11] 

97 ]) 

98 

99 return vertices, faces 

100 

101 

102def neighbors(vertices, faces): 

103 

104 nv = vertices.shape[0] 

105 nf, nc = faces.shape 

106 fedges = num.zeros((nf*nc, 2), dtype=int) 

107 for ic in range(nc): 

108 fedges[ic::nc, :] = faces[:, (ic, (ic+1) % nc)] 

109 

110 indface = num.repeat(num.arange(nf), nc) 

111 fedges1 = fedges[:, 0] + fedges[:, 1]*nv 

112 sortorder = fedges1.argsort() 

113 fedges1_rev = fedges[:, 1] + fedges[:, 0]*nv 

114 

115 inds = num.searchsorted(fedges1, fedges1_rev, sorter=sortorder) 

116 # todo: handle cases without neighbors 

117 assert num.all(fedges1[sortorder[inds]] == fedges1_rev) 

118 neighbors = indface[sortorder[inds]].reshape((nf, nc)) 

119 return neighbors 

120 

121 

122def adjacent_faces(vertices, faces): 

123 

124 nv = vertices.shape[0] 

125 nf, nc = faces.shape 

126 

127 iverts = faces.reshape(nf*nc) 

128 ifaces = num.repeat(num.arange(nf), nc) 

129 iverts_order = iverts.argsort() 

130 

131 vs = vertices[faces] 

132 gs = num.zeros(vs.shape) 

133 gs[:, :-1, :] = vs[:, 1:, :] - vs[:, :-1, :] 

134 gs[:, -1, :] = vs[:, 0, :] - vs[:, -1, :] 

135 

136 vecs = gs.reshape((nf*nc, 3)) 

137 iverts_ordered = iverts[iverts_order] 

138 ifirsts = nextval_indices(iverts_ordered) 

139 iselected = iverts_ordered[ifirsts] 

140 

141 plane_en = normalize(vertices) 

142 plane_e1 = num.zeros((nv, 3)) 

143 plane_e1[iselected, :] = normalize( 

144 project_to_plane_nn(plane_en[iselected, :], vecs[ifirsts, :])) 

145 plane_e2 = num.zeros((nv, 3)) 

146 plane_e2[iselected, :] = num.cross( 

147 plane_en[iselected, :], plane_e1[iselected, :]) 

148 

149 a1 = vdot(vecs[iverts_order, :], plane_e1[iverts_ordered, :]) 

150 a2 = vdot(vecs[iverts_order, :], plane_e2[iverts_ordered, :]) 

151 angles = num.arctan2(a1, a2) * r2d + 360. * iverts_ordered 

152 iverts_order2 = num.argsort(angles) 

153 

154 return iverts_ordered, ifaces[iverts_order][iverts_order2] 

155 

156 

157def nextval_indices(a): 

158 return num.concatenate( 

159 [[0], num.where(num.diff(a) != 0)[0] + 1]) 

160 

161 

162def project_to_plane(vns, vas): 

163 return vas - (vdot(vas, vns) / vdot(vns, vns))[:, num.newaxis] * vns 

164 

165 

166def project_to_plane_nn(vns, vas): 

167 return vas - (vdot(vas, vns))[:, num.newaxis] * vns 

168 

169 

170def corner_handednesses(vertices, faces): 

171 vs = vertices[faces] 

172 gs = num.zeros(vs.shape) 

173 gs[:, :-1, :] = vs[:, 1:, :] - vs[:, :-1, :] 

174 gs[:, -1, :] = vs[:, 0, :] - vs[:, -1, :] 

175 hs = num.zeros((faces.shape[0], faces.shape[1])) 

176 hs[:, 1:] = vdot(num.cross(gs[:, :-1, :], gs[:, 1:, :]), vs[:, 1:, :]) 

177 hs[:, 0] = vdot(num.cross(gs[:, -1, :], gs[:, 0, :]), vs[:, 0, :]) 

178 return hs 

179 

180 

181def truncate(vertices, faces): 

182 

183 nv = vertices.shape[0] 

184 iverts, ifaces = adjacent_faces(vertices, faces) 

185 ifirsts = nextval_indices(iverts) 

186 ilengths = num.zeros(ifirsts.size, dtype=int) 

187 ilengths[:-1] = num.diff(ifirsts) 

188 ilengths[-1] = iverts.size - ifirsts[-1] 

189 nc = num.max(ilengths) 

190 nf = nv 

191 vertices_new = face_centers(vertices, faces) 

192 faces_new = num.zeros((nf, nc), dtype=int) 

193 for iface in range(nf): 

194 ifirst = ifirsts[iface] 

195 ilength = ilengths[iface] 

196 faces_new[iface, :ilength] = ifaces[ifirst:ifirst+ilength] 

197 faces_new[iface, ilength:] = ifaces[ifirst+ilength-1] 

198 

199 faces_new = faces_new[:, ::-1] 

200 

201 return vertices_new, faces_new 

202 

203 

204def iter_icospheres(order, inflate=True): 

205 vertices, faces = icosahedron() 

206 vertices /= vnorm(vertices)[:, num.newaxis] 

207 yield (vertices, faces) 

208 

209 for i in range(order): 

210 vertices, faces = refine_triangles(vertices, faces) 

211 if inflate: 

212 vertices = normalize(vertices) 

213 

214 yield (vertices, faces) 

215 

216 

217bases = { 

218 'icosahedron': icosahedron, 

219 'tetrahedron': tetrahedron, 

220 'tcube': tcube} 

221 

222 

223def sphere( 

224 order, 

225 base=icosahedron, 

226 kind='kind1', 

227 inflate=True, 

228 radius=1.0, 

229 triangulate=True): 

230 

231 if isinstance(base, str): 

232 base = bases[base] 

233 

234 vertices, faces = base() 

235 vertices /= vnorm(vertices)[:, num.newaxis] 

236 

237 for i in range(order): 

238 vertices, faces = refine_triangles(vertices, faces) 

239 if inflate: 

240 vertices = normalize(vertices) 

241 if radius != 1.0: 

242 vertices *= radius 

243 

244 if kind == 'kind2': 

245 vertices, faces = truncate(vertices, faces) 

246 if triangulate: 

247 vertices, faces = triangles_to_center(vertices, faces) 

248 

249 return vertices, faces