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 numpy as num
10from pyrocko.guts import Object, Int
12from pyrocko import icosphere, geometry as g
14guts_prefix = 'sparrow'
16and_ = num.logical_and
19class Level(object):
20 def __init__(self, vertices, faces):
21 self.vertices = vertices
22 self.faces = faces
23 self.next = None
24 self.nvertices = vertices.shape[0]
25 self.nfaces = faces.shape[0]
26 self.prepare_projections()
28 def points_in_face(self, iface, points):
29 uv = g.triangle_barycentric_coordinates(
30 self.bary_trans[iface, :, :],
31 self.vertices[self.faces[iface, 0]],
32 points)
34 return and_(
35 and_(uv[:, 0] >= 0., uv[:, 1] >= 0.), uv[:, 0] + uv[:, 1] <= 1.)
37 def points_to_face(self, points, icandidates=None):
39 if icandidates is None:
40 icandidates = num.arange(self.nfaces)
42 npoints = points.shape[0]
43 ifaces = num.zeros(npoints, dtype=int)
44 points_flat_all = num.zeros((npoints, 3))
45 for icandidate in icandidates:
46 points_flat, mask = self.project(icandidate, points)
47 ipoints = num.where(num.logical_and(
48 self.points_in_face(icandidate, points_flat),
49 mask))[0]
51 if self.next:
52 icandidates_child = num.where(
53 self.next.parents == icandidate)[0]
54 ifaces[ipoints], points_flat_all[ipoints, :] \
55 = self.next.points_to_face(
56 points[ipoints], icandidates_child)
58 else:
59 ifaces[ipoints] = icandidate
60 points_flat_all[ipoints, :] = points_flat[ipoints, :]
62 return ifaces, points_flat_all
64 def project(self, iface, points):
65 denom = g.vdot(points, self.aplanes[iface, :])
66 mask = denom > 0.0
67 points = points.copy()
68 points[mask, :] /= denom[mask, num.newaxis]
69 points[num.logical_not(mask), :] = 0.0
70 return points, mask
72 def prepare_projections(self):
73 planes = g.triangle_planes(
74 self.vertices, self.faces)
76 planes /= g.vdot(
77 planes,
78 self.vertices[self.faces[:, 0]])[:, num.newaxis]
80 self.aplanes = planes
82 self.bary_trans = g.triangle_barycentric_transforms(
83 self.vertices, self.faces)
86class HiMesh(Object):
88 order = Int.T(default=0)
90 def __init__(self, **kwargs):
91 Object.__init__(self, **kwargs)
92 self.levels = self.make_levels()
94 def make_levels(self):
95 # root = None
96 level = None
97 levels = []
98 for vertices, faces in icosphere.iter_icospheres(
99 self.order, inflate=False):
101 new = Level(vertices, faces)
102 if not levels:
103 new.parents = None
104 else:
105 new.parents = g.refine_triangle_parents(faces.shape[0])
106 level.next = new
108 level = new
109 levels.append(new)
111 return levels
113 def get_vertices(self):
114 return g.normalize(self.levels[-1].vertices)
116 def get_faces(self):
117 return self.levels[-1].faces
119 def points_to_faces(self, points):
120 return self.levels[0].points_to_face(points)[0]