Coverage for trimesh/curvature.py: 79%

1""" 

2curvature.py 

3--------------- 

4 

5Query mesh curvature. 

6""" 

7 

8import numpy as np 

9 

10from . import util 

11from .util import diagonal_dot 

12 

13try: 

14 from scipy.sparse import coo_matrix 

15except ImportError as E: 

16 from . import exceptions 

17 

18 coo_matrix = exceptions.ExceptionWrapper(E) 

19 

20 

21def face_angles_sparse(mesh): 

22 """ 

23 A sparse matrix representation of the face angles. 

24 

25 Returns 

26 ---------- 

27 sparse : scipy.sparse.coo_matrix 

28 matrix is float shaped (len(vertices), len(faces)) 

29 """ 

30 matrix = coo_matrix( 

31 (mesh.face_angles.flatten(), (mesh.faces_sparse.row, mesh.faces_sparse.col)), 

32 mesh.faces_sparse.shape, 

33 ) 

34 return matrix 

35 

36 

37def vertex_defects(mesh): 

38 """ 

39 Return the vertex defects, or (2*pi) minus the sum of the 

40 angles of every face that includes that vertex. 

41 

42 If a vertex is only included by coplanar triangles, this 

43 will be zero. For convex regions this is positive, and 

44 concave negative. 

45 

46 Returns 

47 -------- 

48 vertex_defect : (len(self.vertices), ) float 

49 Vertex defect at the every vertex 

50 """ 

51 angle_sum = np.array(mesh.face_angles_sparse.sum(axis=1)).flatten() 

52 defect = (2 * np.pi) - angle_sum 

53 return defect 

54 

55 

56def discrete_gaussian_curvature_measure(mesh, points, radius): 

57 """ 

58 Return the discrete gaussian curvature measure of a sphere 

59 centered at a point as detailed in 'Restricted Delaunay 

60 triangulations and normal cycle'- Cohen-Steiner and Morvan. 

61 

62 This is the sum of the vertex defects at all vertices 

63 within the radius for each point. 

64 

65 Parameters 

66 ---------- 

67 points : (n, 3) float 

68 Points in space 

69 radius : float , 

70 The sphere radius, which can be zero if vertices 

71 passed are points. 

72 

73 Returns 

74 -------- 

75 gaussian_curvature: (n,) float 

76 Discrete gaussian curvature measure. 

77 """ 

78 

79 points = np.asanyarray(points, dtype=np.float64) 

80 if not util.is_shape(points, (-1, 3)): 

81 raise ValueError("points must be (n,3)!") 

82 

83 nearest = mesh.kdtree.query_ball_point(points, radius) 

84 gauss_curv = [mesh.vertex_defects[vertices].sum() for vertices in nearest] 

85 

86 return np.asarray(gauss_curv) 

87 

88 

89def discrete_mean_curvature_measure(mesh, points, radius): 

90 """ 

91 Return the discrete mean curvature measure of a sphere 

92 centered at a point as detailed in 'Restricted Delaunay 

93 triangulations and normal cycle'- Cohen-Steiner and Morvan. 

94 

95 This is the sum of the angle at all edges contained in the 

96 sphere for each point. 

97 

98 Parameters 

99 ---------- 

100 points : (n, 3) float 

101 Points in space 

102 radius : float 

103 Sphere radius which should typically be greater than zero 

104 

105 Returns 

106 -------- 

107 mean_curvature : (n,) float 

108 Discrete mean curvature measure. 

109 """ 

110 

111 points = np.asanyarray(points, dtype=np.float64) 

112 if not util.is_shape(points, (-1, 3)): 

113 raise ValueError("points must be (n,3)!") 

114 

115 # resolve the cached properties once rather than once per query point 

116 vertices = np.asarray(mesh.vertices) 

117 adjacency_edges = np.asarray(mesh.face_adjacency_edges) 

118 adjacency_angles = np.asarray(mesh.face_adjacency_angles) 

119 adjacency_convex = np.asarray(mesh.face_adjacency_convex) 

120 

121 tree = mesh.face_adjacency_tree 

122 # axis aligned bounds 

123 mins = points - radius 

124 maxs = points + radius 

125 

126 try: 

127 # use the batch API added in 1.4.0 and fixed to actually work in 1.4.1 

128 hit_ids, hit_counts = tree.intersection_v(mins, maxs) 

129 candidates = np.asarray(hit_ids, dtype=np.int64) 

130 counts = np.asarray(hit_counts, dtype=np.int64) 

131 except BaseException: 

132 # fall back to a list comprehension 

133 per_point = [list(tree.intersection(b)) for b in np.column_stack((mins, maxs))] 

134 counts = np.array([len(c) for c in per_point], dtype=np.int64) 

135 candidates = np.fromiter( 

136 (i for c in per_point for i in c), dtype=np.int64, count=int(counts.sum()) 

137 ) 

138 

139 if len(candidates) == 0: 

140 return np.zeros(len(points)) 

141 

142 # the index of the query point each candidate edge belongs to 

143 owner = np.repeat(np.arange(len(points)), counts) 

144 endpoints = vertices[adjacency_edges[candidates]] 

145 

146 # `line_ball_intersection` broadcasts a per-row center already 

147 lengths = line_ball_intersection( 

148 endpoints[:, 0], endpoints[:, 1], center=points[owner], radius=radius 

149 ) 

150 signs = np.where(adjacency_convex[candidates], 1, -1) 

151 

152 # sum the contribution of every candidate edge into its own query point 

153 return ( 

154 np.bincount( 

155 owner, 

156 weights=lengths * adjacency_angles[candidates] * signs, 

157 minlength=len(points), 

158 ) 

159 / 2 

160 ) 

161 

162 

163def line_ball_intersection(start_points, end_points, center, radius): 

164 """ 

165 Compute the length of the intersection of a line segment with a ball. 

166 

167 Parameters 

168 ---------- 

169 start_points : (n,3) float, list of points in space 

170 end_points : (n,3) float, list of points in space 

171 center : (3,) float, the sphere center 

172 radius : float, the sphere radius 

173 

174 Returns 

175 -------- 

176 lengths: (n,) float, the lengths. 

177 

178 """ 

179 

180 # We solve for the intersection of |x-c|**2 = r**2 and 

181 # x = o + dL. This yields 

182 # d = (-l.(o-c) +- sqrt[ l.(o-c)**2 - l.l((o-c).(o-c) - r^**2) ]) / l.l 

183 L = end_points - start_points 

184 oc = start_points - center # o-c 

185 r = radius 

186 ldotl = diagonal_dot(L, L) 

187 ldotoc = diagonal_dot(L, oc) 

188 ocdotoc = diagonal_dot(oc, oc) 

189 discrims = ldotoc**2 - ldotl * (ocdotoc - r**2) 

190 

191 # If discriminant is non-positive, then we have zero length 

192 lengths = np.zeros(len(start_points)) 

193 # Otherwise we solve for the solns with d2 > d1. 

194 m = discrims > 0 # mask 

195 d1 = (-ldotoc[m] - np.sqrt(discrims[m])) / ldotl[m] 

196 d2 = (-ldotoc[m] + np.sqrt(discrims[m])) / ldotl[m] 

197 

198 # Line segment means we have 0 <= d <= 1 

199 d1 = np.clip(d1, 0, 1) 

200 d2 = np.clip(d2, 0, 1) 

201 

202 # Length is |o + d2 l - o + d1 l| = (d2 - d1) |l| 

203 lengths[m] = (d2 - d1) * np.sqrt(ldotl[m]) 

204 

205 return lengths 

206 

207 

208def sphere_ball_intersection(R, r): 

209 """ 

210 Compute the surface area of the intersection of sphere of radius R centered 

211 at (0, 0, 0) with a ball of radius r centered at (R, 0, 0). 

212 

213 Parameters 

214 ---------- 

215 R : float, sphere radius 

216 r : float, ball radius 

217 

218 Returns 

219 -------- 

220 area: float, the surface are. 

221 """ 

222 x = (2 * R**2 - r**2) / (2 * R) # x coord of plane 

223 if x >= -R: 

224 return 2 * np.pi * R * (R - x) 

225 if x < -R: 

226 return 4 * np.pi * R**2