Coverage for trimesh/smoothing.py: 93%

112 statements  

« prev     ^ index     » next       coverage.py v7.14.1, created at 2026-06-24 04:40 +0000

1import numpy as np 

2 

3from .typed import ArrayLike 

4 

5try: 

6 from scipy.sparse import coo_matrix, eye 

7 from scipy.sparse.linalg import spsolve 

8except ImportError as E: 

9 from .exceptions import ExceptionWrapper 

10 

11 wrapper = ExceptionWrapper(E) 

12 eye, spsolve, coo_matrix = wrapper, wrapper, wrapper 

13 

14from . import graph, triangles 

15from .base import Trimesh 

16from .geometry import index_sparse 

17from .triangles import mass_properties 

18from .util import unitize 

19 

20 

21def filter_laplacian( 

22 mesh, 

23 lamb=0.5, 

24 iterations=10, 

25 implicit_time_integration=False, 

26 volume_constraint=True, 

27 laplacian_operator=None, 

28): 

29 """ 

30 Smooth a mesh in-place using laplacian smoothing. 

31 Articles 

32 1 - "Improved Laplacian Smoothing of Noisy Surface Meshes" 

33 J. Vollmer, R. Mencl, and H. Muller 

34 2 - "Implicit Fairing of Irregular Meshes using Diffusion 

35 and Curvature Flow". M. Desbrun, M. Meyer, 

36 P. Schroder, A.H.B. Caltech 

37 Parameters 

38 ------------ 

39 mesh : trimesh.Trimesh 

40 Mesh to be smoothed in place 

41 lamb : float 

42 Diffusion speed constant 

43 If 0.0, no diffusion 

44 If > 0.0, diffusion occurs 

45 implicit_time_integration: boolean 

46 if False: explicit time integration 

47 -lamb <= 1.0 - Stability Limit (Article 1) 

48 if True: implicit time integration 

49 -lamb no limit (Article 2) 

50 iterations : int 

51 Number of passes to run filter 

52 laplacian_operator : None or scipy.sparse.coo.coo_matrix 

53 Sparse matrix laplacian operator 

54 Will be autogenerated if None 

55 """ 

56 

57 # if the laplacian operator was not passed create it here 

58 if laplacian_operator is None: 

59 laplacian_operator = laplacian_calculation(mesh) 

60 

61 # save initial volume 

62 if volume_constraint: 

63 vol_ini = mesh.volume 

64 

65 # get mesh vertices and faces as vanilla numpy array 

66 vertices = mesh.vertices.copy().view(np.ndarray) 

67 faces = mesh.faces.copy().view(np.ndarray) 

68 

69 # Set matrix for linear system of equations 

70 if implicit_time_integration: 

71 dlap = laplacian_operator.shape[0] 

72 AA = eye(dlap) + lamb * (eye(dlap) - laplacian_operator) 

73 

74 # Number of passes 

75 for _index in range(iterations): 

76 # Classic Explicit Time Integration - Article 1 

77 if not implicit_time_integration: 

78 dot = laplacian_operator.dot(vertices) - vertices 

79 vertices += lamb * dot 

80 

81 # Implicit Time Integration - Article 2 

82 else: 

83 vertices = spsolve(AA, vertices) 

84 

85 # volume constraint 

86 if volume_constraint: 

87 # find the volume with new vertex positions 

88 vol_new = triangles.mass_properties(vertices[faces], skip_inertia=True)[ 

89 "volume" 

90 ] 

91 # scale by volume ratio 

92 vertices *= (vol_ini / vol_new) ** (1.0 / 3.0) 

93 

94 # assign modified vertices back to mesh 

95 mesh.vertices = vertices 

96 return mesh 

97 

98 

99def filter_humphrey(mesh, alpha=0.1, beta=0.5, iterations=10, laplacian_operator=None): 

100 """ 

101 Smooth a mesh in-place using laplacian smoothing 

102 and Humphrey filtering. 

103 Articles 

104 "Improved Laplacian Smoothing of Noisy Surface Meshes" 

105 J. Vollmer, R. Mencl, and H. Muller 

106 Parameters 

107 ------------ 

108 mesh : trimesh.Trimesh 

109 Mesh to be smoothed in place 

110 alpha : float 

111 Controls shrinkage, range is 0.0 - 1.0 

112 If 0.0, not considered 

113 If 1.0, no smoothing 

114 beta : float 

115 Controls how aggressive smoothing is 

116 If 0.0, no smoothing 

117 If 1.0, full aggressiveness 

118 iterations : int 

119 Number of passes to run filter 

120 laplacian_operator : None or scipy.sparse.coo.coo_matrix 

121 Sparse matrix laplacian operator 

122 Will be autogenerated if None 

123 """ 

124 # if the laplacian operator was not passed create it here 

125 if laplacian_operator is None: 

126 laplacian_operator = laplacian_calculation(mesh) 

127 

128 # get mesh vertices as vanilla numpy array 

129 vertices = mesh.vertices.copy().view(np.ndarray) 

130 # save original unmodified vertices 

131 original = vertices.copy() 

132 

133 # run through iterations of filter 

134 for _index in range(iterations): 

135 vert_q = vertices.copy() 

136 vertices = laplacian_operator.dot(vertices) 

137 vert_b = vertices - (alpha * original + (1.0 - alpha) * vert_q) 

138 vertices -= beta * vert_b + (1.0 - beta) * laplacian_operator.dot(vert_b) 

139 

140 # assign modified vertices back to mesh 

141 mesh.vertices = vertices 

142 return mesh 

143 

144 

145def filter_taubin(mesh, lamb=0.5, nu=0.5, iterations=10, laplacian_operator=None): 

146 """ 

147 Smooth a mesh in-place using laplacian smoothing 

148 and taubin filtering. 

149 Articles 

150 "Improved Laplacian Smoothing of Noisy Surface Meshes" 

151 J. Vollmer, R. Mencl, and H. Muller 

152 Parameters 

153 ------------ 

154 mesh : trimesh.Trimesh 

155 Mesh to be smoothed in place. 

156 lamb : float 

157 Controls shrinkage, range is 0.0 - 1.0 

158 nu : float 

159 Controls dilation, range is 0.0 - 1.0 

160 Nu shall be between 0.0 < 1.0/lambda - 1.0/nu < 0.1 

161 iterations : int 

162 Number of passes to run the filter 

163 laplacian_operator : None or scipy.sparse.coo.coo_matrix 

164 Sparse matrix laplacian operator 

165 Will be autogenerated if None 

166 """ 

167 # if the laplacian operator was not passed create it here 

168 if laplacian_operator is None: 

169 laplacian_operator = laplacian_calculation(mesh) 

170 

171 # get mesh vertices as vanilla numpy array 

172 vertices = mesh.vertices.copy().view(np.ndarray) 

173 

174 # run through multiple passes of the filter 

175 for index in range(iterations): 

176 # do a sparse dot product on the vertices 

177 dot = laplacian_operator.dot(vertices) - vertices 

178 # alternate shrinkage and dilation 

179 if index % 2 == 0: 

180 vertices += lamb * dot 

181 else: 

182 vertices -= nu * dot 

183 

184 # assign updated vertices back to mesh 

185 mesh.vertices = vertices 

186 return mesh 

187 

188 

189def filter_mut_dif_laplacian( 

190 mesh, lamb=0.5, iterations=10, volume_constraint=True, laplacian_operator=None 

191): 

192 """ 

193 Smooth a mesh in-place using laplacian smoothing using a 

194 mutable diffusion laplacian. 

195 

196 Articles 

197 Barroqueiro, B., Andrade-Campos, A., Dias-de-Oliveira, 

198 J., and Valente, R. (January 21, 2021). 

199 "Bridging between topology optimization and additive 

200 manufacturing via Laplacian smoothing." ASME. J. Mech. Des. 

201 

202 

203 Parameters 

204 ------------ 

205 mesh : trimesh.Trimesh 

206 Mesh to be smoothed in place 

207 lamb : float 

208 Diffusion speed constant 

209 If 0.0, no diffusion 

210 If > 0.0, diffusion occurs 

211 iterations : int 

212 Number of passes to run filter 

213 laplacian_operator : None or scipy.sparse.coo.coo_matrix 

214 Sparse matrix laplacian operator 

215 Will be autogenerated if None 

216 """ 

217 

218 # if the laplacian operator was not passed create it here 

219 if laplacian_operator is None: 

220 laplacian_operator = laplacian_calculation(mesh) 

221 

222 # Set volume constraint 

223 if volume_constraint: 

224 v_ini = mesh.volume 

225 

226 # get mesh vertices as vanilla numpy array 

227 vertices = mesh.vertices.copy().view(np.ndarray) 

228 faces = mesh.faces.copy().view(np.ndarray) 

229 eps = 0.01 * (np.max(mesh.area_faces) ** 0.5) 

230 

231 # Number of passes 

232 for _index in range(iterations): 

233 # Mutable diffusion 

234 normals = get_vertices_normals(mesh) 

235 qi = laplacian_operator.dot(vertices) 

236 pi_qi = vertices - qi 

237 

238 adil = np.abs((normals * pi_qi).dot(np.ones((3, 1)))) 

239 adil = 1.0 / np.maximum(1e-12, adil) 

240 lamber = np.maximum(0.2 * lamb, np.minimum(1.0, lamb * adil / np.mean(adil))) 

241 

242 # Filter 

243 dot = laplacian_operator.dot(vertices) 

244 vertices += lamber * (dot - vertices) 

245 

246 # Volume constraint 

247 if volume_constraint: 

248 vol = mass_properties(vertices[faces], skip_inertia=True)["volume"] 

249 if _index == 0: 

250 slope = dilate_slope(vertices, faces, normals, vol, eps) 

251 vertices += normals * slope * (v_ini - vol) 

252 

253 # assign modified vertices back to mesh 

254 mesh.vertices = vertices 

255 

256 return mesh 

257 

258 

259def laplacian_calculation( 

260 mesh: Trimesh, 

261 equal_weight: bool = True, 

262 pinned_vertices: ArrayLike | None = None, 

263): 

264 """ 

265 Calculate a sparse matrix for laplacian operations. 

266 

267 Note that setting equal_weight to False significantly hampers performance. 

268 

269 Parameters 

270 ------------- 

271 mesh : trimesh.Trimesh 

272 Input geometry 

273 equal_weight : bool 

274 If True, all neighbors will be considered equally 

275 If False, all neighbors will be weighted by inverse distance 

276 pinned_vertices : None or list of ints 

277 If None, no vertices are pinned 

278 If list, vertices will be pinned, such that they will not be moved 

279 Returns 

280 ---------- 

281 laplacian : scipy.sparse.coo.coo_matrix 

282 Laplacian operator 

283 """ 

284 if equal_weight: 

285 laplacian = graph.edges_to_coo(mesh.edges) 

286 

287 if pinned_vertices is not None: 

288 # Set pinned vertices to only have themselves as neighbours 

289 laplacian.data[np.isin(laplacian.row, pinned_vertices)] = False 

290 laplacian.row = np.concatenate((laplacian.row, pinned_vertices)) 

291 laplacian.col = np.concatenate((laplacian.col, pinned_vertices)) 

292 laplacian.data = np.concatenate( 

293 (laplacian.data, np.ones(len(pinned_vertices), dtype=bool)) 

294 ) 

295 

296 laplacian = laplacian / laplacian.sum(axis=1) 

297 

298 else: 

299 # get the vertex neighbors from the cache 

300 neighbors = mesh.vertex_neighbors 

301 

302 # if a node is pinned, it will average his coordinates by himself 

303 # in practice it will not move 

304 if pinned_vertices is not None: 

305 for i in pinned_vertices: 

306 neighbors[i] = [i] 

307 

308 # avoid hitting crc checks in loops 

309 vertices = mesh.vertices.view(np.ndarray) 

310 

311 # stack neighbors to 1D arrays 

312 col = np.concatenate(neighbors) 

313 row = np.concatenate([[i] * len(n) for i, n in enumerate(neighbors)]) 

314 

315 # umbrella weights, distance-weighted 

316 # use dot product of ones to replace array.sum(axis=1) 

317 ones = np.ones(3) 

318 # the distance from verticesex to neighbors 

319 norms = [ 

320 1.0 

321 / np.maximum(1e-6, np.sqrt(np.dot((vertices[i] - vertices[n]) ** 2, ones))) 

322 for i, n in enumerate(neighbors) 

323 ] 

324 # normalize group and stack into single array 

325 data = np.concatenate([i / i.sum() for i in norms]) 

326 

327 # create the sparse matrix 

328 laplacian = coo_matrix((data, (row, col)), shape=[len(vertices)] * 2) 

329 

330 return laplacian 

331 

332 

333def get_vertices_normals(mesh): 

334 """ 

335 Compute Vertex normals using equal weighting of neighbors faces. 

336 Parameters 

337 ------------- 

338 mesh : trimesh.Trimesh 

339 Input geometry 

340 Returns 

341 ---------- 

342 vertices_normals: array 

343 Vertices normals 

344 """ 

345 

346 # get mesh vertices and faces 

347 vertices = mesh.vertices 

348 faces = mesh.faces 

349 

350 # get face normals 

351 face_normals = mesh.face_normals 

352 

353 # Compute Vert normals 

354 vert_normals = index_sparse(len(vertices), faces).dot(face_normals) 

355 

356 return unitize(vert_normals) 

357 

358 

359def dilate_slope(vertices, faces, normals, v, eps): 

360 """ 

361 Get the derivate of dilation scalar by the volume variation by finite differences 

362 Thus, Vertices += vertex_normals*dilate_slope*(Initial_Volume - Srinked_Volume) 

363 Parameters 

364 ------------- 

365 mesh : trimesh.Trimesh 

366 Input geometry 

367 vertices: mesh.vertices 

368 faces: mesh.faces 

369 normals: array 

370 vertices normals 

371 Returns 

372 ---------- 

373 dilate_slope: float 

374 derivative 

375 """ 

376 

377 # finite difference derivative 

378 vertices2 = vertices + normals * eps 

379 v2 = mass_properties(vertices2[faces], skip_inertia=True)["volume"] 

380 

381 return (eps) / (v2 - v)