Coverage for trimesh/curvature.py: 83%
58 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-24 04:40 +0000
1"""
2curvature.py
3---------------
5Query mesh curvature.
6"""
8import numpy as np
10from . import util
11from .util import diagonal_dot
13try:
14 from scipy.sparse import coo_matrix
15except ImportError as E:
16 from . import exceptions
18 coo_matrix = exceptions.ExceptionWrapper(E)
21def face_angles_sparse(mesh):
22 """
23 A sparse matrix representation of the face angles.
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
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.
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.
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
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.
62 This is the sum of the vertex defects at all vertices
63 within the radius for each point.
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.
73 Returns
74 --------
75 gaussian_curvature: (n,) float
76 Discrete gaussian curvature measure.
77 """
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)!")
83 nearest = mesh.kdtree.query_ball_point(points, radius)
84 gauss_curv = [mesh.vertex_defects[vertices].sum() for vertices in nearest]
86 return np.asarray(gauss_curv)
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.
95 This is the sum of the angle at all edges contained in the
96 sphere for each point.
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
105 Returns
106 --------
107 mean_curvature : (n,) float
108 Discrete mean curvature measure.
109 """
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)!")
115 # axis aligned bounds
116 bounds = np.column_stack((points - radius, points + radius))
118 # line segments that intersect axis aligned bounding box
119 candidates = [list(mesh.face_adjacency_tree.intersection(b)) for b in bounds]
121 mean_curv = np.zeros(len(points))
122 for i, (x, x_candidates) in enumerate(zip(points, candidates)):
123 endpoints = mesh.vertices[mesh.face_adjacency_edges[x_candidates]]
124 lengths = line_ball_intersection(
125 endpoints[:, 0], endpoints[:, 1], center=x, radius=radius
126 )
127 angles = mesh.face_adjacency_angles[x_candidates]
128 signs = np.where(mesh.face_adjacency_convex[x_candidates], 1, -1)
129 mean_curv[i] = (lengths * angles * signs).sum() / 2
131 return mean_curv
134def line_ball_intersection(start_points, end_points, center, radius):
135 """
136 Compute the length of the intersection of a line segment with a ball.
138 Parameters
139 ----------
140 start_points : (n,3) float, list of points in space
141 end_points : (n,3) float, list of points in space
142 center : (3,) float, the sphere center
143 radius : float, the sphere radius
145 Returns
146 --------
147 lengths: (n,) float, the lengths.
149 """
151 # We solve for the intersection of |x-c|**2 = r**2 and
152 # x = o + dL. This yields
153 # d = (-l.(o-c) +- sqrt[ l.(o-c)**2 - l.l((o-c).(o-c) - r^**2) ]) / l.l
154 L = end_points - start_points
155 oc = start_points - center # o-c
156 r = radius
157 ldotl = diagonal_dot(L, L)
158 ldotoc = diagonal_dot(L, oc)
159 ocdotoc = diagonal_dot(oc, oc)
160 discrims = ldotoc**2 - ldotl * (ocdotoc - r**2)
162 # If discriminant is non-positive, then we have zero length
163 lengths = np.zeros(len(start_points))
164 # Otherwise we solve for the solns with d2 > d1.
165 m = discrims > 0 # mask
166 d1 = (-ldotoc[m] - np.sqrt(discrims[m])) / ldotl[m]
167 d2 = (-ldotoc[m] + np.sqrt(discrims[m])) / ldotl[m]
169 # Line segment means we have 0 <= d <= 1
170 d1 = np.clip(d1, 0, 1)
171 d2 = np.clip(d2, 0, 1)
173 # Length is |o + d2 l - o + d1 l| = (d2 - d1) |l|
174 lengths[m] = (d2 - d1) * np.sqrt(ldotl[m])
176 return lengths
179def sphere_ball_intersection(R, r):
180 """
181 Compute the surface area of the intersection of sphere of radius R centered
182 at (0, 0, 0) with a ball of radius r centered at (R, 0, 0).
184 Parameters
185 ----------
186 R : float, sphere radius
187 r : float, ball radius
189 Returns
190 --------
191 area: float, the surface are.
192 """
193 x = (2 * R**2 - r**2) / (2 * R) # x coord of plane
194 if x >= -R:
195 return 2 * np.pi * R * (R - x)
196 if x < -R:
197 return 4 * np.pi * R**2