1# https://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import numpy as num
8from pyrocko.guts import Object, Int
10from pyrocko import icosphere, geometry as g
12guts_prefix = 'sparrow'
14and_ = num.logical_and
17class Level(object):
18 def __init__(self, vertices, faces):
19 self.vertices = vertices
20 self.faces = faces
21 self.next = None
22 self.nvertices = vertices.shape[0]
23 self.nfaces = faces.shape[0]
24 self.prepare_projections()
26 def points_in_face(self, iface, points):
27 uv = g.triangle_barycentric_coordinates(
28 self.bary_trans[iface, :, :],
29 self.vertices[self.faces[iface, 0]],
30 points)
32 return and_(
33 and_(uv[:, 0] >= 0., uv[:, 1] >= 0.), uv[:, 0] + uv[:, 1] <= 1.)
35 def points_to_face(self, points, icandidates=None):
37 if icandidates is None:
38 icandidates = num.arange(self.nfaces)
40 npoints = points.shape[0]
41 ifaces = num.zeros(npoints, dtype=int)
42 points_flat_all = num.zeros((npoints, 3))
43 for icandidate in icandidates:
44 points_flat, mask = self.project(icandidate, points)
45 ipoints = num.where(num.logical_and(
46 self.points_in_face(icandidate, points_flat),
47 mask))[0]
49 if self.next:
50 icandidates_child = num.where(
51 self.next.parents == icandidate)[0]
52 ifaces[ipoints], points_flat_all[ipoints, :] \
53 = self.next.points_to_face(
54 points[ipoints], icandidates_child)
56 else:
57 ifaces[ipoints] = icandidate
58 points_flat_all[ipoints, :] = points_flat[ipoints, :]
60 return ifaces, points_flat_all
62 def project(self, iface, points):
63 denom = g.vdot(points, self.aplanes[iface, :])
64 mask = denom > 0.0
65 points = points.copy()
66 points[mask, :] /= denom[mask, num.newaxis]
67 points[num.logical_not(mask), :] = 0.0
68 return points, mask
70 def prepare_projections(self):
71 planes = g.triangle_planes(
72 self.vertices, self.faces)
74 planes /= g.vdot(
75 planes,
76 self.vertices[self.faces[:, 0]])[:, num.newaxis]
78 self.aplanes = planes
80 self.bary_trans = g.triangle_barycentric_transforms(
81 self.vertices, self.faces)
84class HiMesh(Object):
86 order = Int.T(default=0)
88 def __init__(self, **kwargs):
89 Object.__init__(self, **kwargs)
90 self.levels = self.make_levels()
92 def make_levels(self):
93 # root = None
94 level = None
95 levels = []
96 for vertices, faces in icosphere.iter_icospheres(
97 self.order, inflate=False):
99 new = Level(vertices, faces)
100 if not levels:
101 new.parents = None
102 else:
103 new.parents = g.refine_triangle_parents(faces.shape[0])
104 level.next = new
106 level = new
107 levels.append(new)
109 return levels
111 def get_vertices(self):
112 return g.normalize(self.levels[-1].vertices)
114 def get_faces(self):
115 return self.levels[-1].faces
117 def points_to_faces(self, points):
118 return self.levels[0].points_to_face(points)[0]