1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

# https://pyrocko.org - GPLv3 

# 

# The Pyrocko Developers, 21st Century 

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

 

from __future__ import absolute_import, print_function, division 

 

import numpy as num 

 

from pyrocko.guts import Object, Int 

 

from pyrocko import icosphere, geometry as g 

 

guts_prefix = 'sparrow' 

 

and_ = num.logical_and 

 

 

class Level(object): 

def __init__(self, vertices, faces): 

self.vertices = vertices 

self.faces = faces 

self.next = None 

self.nvertices = vertices.shape[0] 

self.nfaces = faces.shape[0] 

self.prepare_projections() 

 

def points_in_face(self, iface, points): 

uv = g.triangle_barycentric_coordinates( 

self.bary_trans[iface, :, :], 

self.vertices[self.faces[iface, 0]], 

points) 

 

return and_( 

and_(uv[:, 0] >= 0., uv[:, 1] >= 0.), uv[:, 0] + uv[:, 1] <= 1.) 

 

def points_to_face(self, points, icandidates=None): 

 

if icandidates is None: 

icandidates = num.arange(self.nfaces) 

 

npoints = points.shape[0] 

ifaces = num.zeros(npoints, dtype=num.int) 

points_flat_all = num.zeros((npoints, 3)) 

for icandidate in icandidates: 

points_flat, mask = self.project(icandidate, points) 

ipoints = num.where(num.logical_and( 

self.points_in_face(icandidate, points_flat), 

mask))[0] 

 

if self.next: 

icandidates_child = num.where( 

self.next.parents == icandidate)[0] 

ifaces[ipoints], points_flat_all[ipoints, :] \ 

= self.next.points_to_face( 

points[ipoints], icandidates_child) 

 

else: 

ifaces[ipoints] = icandidate 

points_flat_all[ipoints, :] = points_flat[ipoints, :] 

 

return ifaces, points_flat_all 

 

def project(self, iface, points): 

denom = g.vdot(points, self.aplanes[iface, :]) 

mask = denom > 0.0 

points = points.copy() 

points[mask, :] /= denom[mask, num.newaxis] 

points[num.logical_not(mask), :] = 0.0 

return points, mask 

 

def prepare_projections(self): 

planes = g.triangle_planes( 

self.vertices, self.faces) 

 

planes /= g.vdot( 

planes, 

self.vertices[self.faces[:, 0]])[:, num.newaxis] 

 

self.aplanes = planes 

 

self.bary_trans = g.triangle_barycentric_transforms( 

self.vertices, self.faces) 

 

 

class HiMesh(Object): 

 

order = Int.T(default=0) 

 

def __init__(self, **kwargs): 

Object.__init__(self, **kwargs) 

self.levels = self.make_levels() 

 

def make_levels(self): 

# root = None 

level = None 

levels = [] 

for vertices, faces in icosphere.iter_icospheres( 

self.order, inflate=False): 

 

new = Level(vertices, faces) 

if not levels: 

new.parents = None 

else: 

new.parents = g.refine_triangle_parents(faces.shape[0]) 

level.next = new 

 

level = new 

levels.append(new) 

 

return levels 

 

def get_vertices(self): 

return g.normalize(self.levels[-1].vertices) 

 

def get_faces(self): 

return self.levels[-1].faces 

 

def points_to_faces(self, points): 

return self.levels[0].points_to_face(points)[0]