vedo.mesh
Submodule to work with polygonal meshes
1#!/usr/bin/env python3 2# -*- coding: utf-8 -*- 3import numpy as np 4from typing import List, Tuple, Union, MutableSequence, Any 5from typing_extensions import Self 6 7import vedo.vtkclasses as vtki # a wrapper for lazy imports 8 9import vedo 10from vedo.colors import get_color 11from vedo.pointcloud import Points 12from vedo.utils import buildPolyData, is_sequence, mag, precision 13from vedo.utils import numpy2vtk, vtk2numpy, OperationNode 14from vedo.visual import MeshVisual 15 16__docformat__ = "google" 17 18__doc__ = """ 19Submodule to work with polygonal meshes 20 21 22""" 23 24__all__ = ["Mesh"] 25 26 27#################################################### 28class Mesh(MeshVisual, Points): 29 """ 30 Build an instance of object `Mesh` derived from `vedo.PointCloud`. 31 """ 32 33 def __init__(self, inputobj=None, c="gold", alpha=1): 34 """ 35 Initialize a ``Mesh`` object. 36 37 Arguments: 38 inputobj : (str, vtkPolyData, vtkActor, vedo.Mesh) 39 If inputobj is `None` an empty mesh is created. 40 If inputobj is a `str` then it is interpreted as the name of a file to load as mesh. 41 If inputobj is an `vtkPolyData` or `vtkActor` or `vedo.Mesh` 42 then a shallow copy of it is created. 43 If inputobj is a `vedo.Mesh` then a shallow copy of it is created. 44 45 Examples: 46 - [buildmesh.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/buildmesh.py) 47 (and many others!) 48 49  50 """ 51 # print("INIT MESH", super()) 52 super().__init__() 53 54 self.name = "Mesh" 55 56 if inputobj is None: 57 # self.dataset = vtki.vtkPolyData() 58 pass 59 60 elif isinstance(inputobj, str): 61 self.dataset = vedo.file_io.load(inputobj).dataset 62 self.filename = inputobj 63 64 elif isinstance(inputobj, vtki.vtkPolyData): 65 # self.dataset.DeepCopy(inputobj) # NO 66 self.dataset = inputobj 67 if self.dataset.GetNumberOfCells() == 0: 68 carr = vtki.vtkCellArray() 69 for i in range(inputobj.GetNumberOfPoints()): 70 carr.InsertNextCell(1) 71 carr.InsertCellPoint(i) 72 self.dataset.SetVerts(carr) 73 74 elif isinstance(inputobj, Mesh): 75 self.dataset = inputobj.dataset 76 77 elif is_sequence(inputobj): 78 ninp = len(inputobj) 79 if ninp == 4: # assume input is [vertices, faces, lines, strips] 80 self.dataset = buildPolyData(inputobj[0], inputobj[1], inputobj[2], inputobj[3]) 81 elif ninp == 3: # assume input is [vertices, faces, lines] 82 self.dataset = buildPolyData(inputobj[0], inputobj[1], inputobj[2]) 83 elif ninp == 2: # assume input is [vertices, faces] 84 self.dataset = buildPolyData(inputobj[0], inputobj[1]) 85 elif ninp == 1: # assume input is [vertices] 86 self.dataset = buildPolyData(inputobj[0]) 87 else: 88 vedo.logger.error("input must be a list of max 4 elements.") 89 raise ValueError() 90 91 elif isinstance(inputobj, vtki.vtkActor): 92 self.dataset.DeepCopy(inputobj.GetMapper().GetInput()) 93 v = inputobj.GetMapper().GetScalarVisibility() 94 self.mapper.SetScalarVisibility(v) 95 pr = vtki.vtkProperty() 96 pr.DeepCopy(inputobj.GetProperty()) 97 self.actor.SetProperty(pr) 98 self.properties = pr 99 100 elif isinstance(inputobj, (vtki.vtkStructuredGrid, vtki.vtkRectilinearGrid)): 101 gf = vtki.new("GeometryFilter") 102 gf.SetInputData(inputobj) 103 gf.Update() 104 self.dataset = gf.GetOutput() 105 106 elif "meshlab" in str(type(inputobj)): 107 self.dataset = vedo.utils.meshlab2vedo(inputobj).dataset 108 109 elif "meshlib" in str(type(inputobj)): 110 import meshlib.mrmeshnumpy as mrmeshnumpy 111 self.dataset = buildPolyData( 112 mrmeshnumpy.getNumpyVerts(inputobj), 113 mrmeshnumpy.getNumpyFaces(inputobj.topology), 114 ) 115 116 elif "trimesh" in str(type(inputobj)): 117 self.dataset = vedo.utils.trimesh2vedo(inputobj).dataset 118 119 elif "meshio" in str(type(inputobj)): 120 # self.dataset = vedo.utils.meshio2vedo(inputobj) ##TODO 121 if len(inputobj.cells) > 0: 122 mcells = [] 123 for cellblock in inputobj.cells: 124 if cellblock.type in ("triangle", "quad"): 125 mcells += cellblock.data.tolist() 126 self.dataset = buildPolyData(inputobj.points, mcells) 127 else: 128 self.dataset = buildPolyData(inputobj.points, None) 129 # add arrays: 130 try: 131 if len(inputobj.point_data) > 0: 132 for k in inputobj.point_data.keys(): 133 vdata = numpy2vtk(inputobj.point_data[k]) 134 vdata.SetName(str(k)) 135 self.dataset.GetPointData().AddArray(vdata) 136 except AssertionError: 137 print("Could not add meshio point data, skip.") 138 139 else: 140 try: 141 gf = vtki.new("GeometryFilter") 142 gf.SetInputData(inputobj) 143 gf.Update() 144 self.dataset = gf.GetOutput() 145 except: 146 vedo.logger.error(f"cannot build mesh from type {type(inputobj)}") 147 raise RuntimeError() 148 149 self.mapper.SetInputData(self.dataset) 150 self.actor.SetMapper(self.mapper) 151 152 self.properties.SetInterpolationToPhong() 153 self.properties.SetColor(get_color(c)) 154 155 if alpha is not None: 156 self.properties.SetOpacity(alpha) 157 158 self.mapper.SetInterpolateScalarsBeforeMapping( 159 vedo.settings.interpolate_scalars_before_mapping 160 ) 161 162 if vedo.settings.use_polygon_offset: 163 self.mapper.SetResolveCoincidentTopologyToPolygonOffset() 164 pof = vedo.settings.polygon_offset_factor 165 pou = vedo.settings.polygon_offset_units 166 self.mapper.SetResolveCoincidentTopologyPolygonOffsetParameters(pof, pou) 167 168 n = self.dataset.GetNumberOfPoints() 169 self.pipeline = OperationNode(self, comment=f"#pts {n}") 170 171 def _repr_html_(self): 172 """ 173 HTML representation of the Mesh object for Jupyter Notebooks. 174 175 Returns: 176 HTML text with the image and some properties. 177 """ 178 import io 179 import base64 180 from PIL import Image 181 182 library_name = "vedo.mesh.Mesh" 183 help_url = "https://vedo.embl.es/docs/vedo/mesh.html#Mesh" 184 185 arr = self.thumbnail() 186 im = Image.fromarray(arr) 187 buffered = io.BytesIO() 188 im.save(buffered, format="PNG", quality=100) 189 encoded = base64.b64encode(buffered.getvalue()).decode("utf-8") 190 url = "data:image/png;base64," + encoded 191 image = f"<img src='{url}'></img>" 192 193 bounds = "<br/>".join( 194 [ 195 precision(min_x, 4) + " ... " + precision(max_x, 4) 196 for min_x, max_x in zip(self.bounds()[::2], self.bounds()[1::2]) 197 ] 198 ) 199 average_size = "{size:.3f}".format(size=self.average_size()) 200 201 help_text = "" 202 if self.name: 203 help_text += f"<b> {self.name}:   </b>" 204 help_text += '<b><a href="' + help_url + '" target="_blank">' + library_name + "</a></b>" 205 if self.filename: 206 dots = "" 207 if len(self.filename) > 30: 208 dots = "..." 209 help_text += f"<br/><code><i>({dots}{self.filename[-30:]})</i></code>" 210 211 pdata = "" 212 if self.dataset.GetPointData().GetScalars(): 213 if self.dataset.GetPointData().GetScalars().GetName(): 214 name = self.dataset.GetPointData().GetScalars().GetName() 215 pdata = "<tr><td><b> point data array </b></td><td>" + name + "</td></tr>" 216 217 cdata = "" 218 if self.dataset.GetCellData().GetScalars(): 219 if self.dataset.GetCellData().GetScalars().GetName(): 220 name = self.dataset.GetCellData().GetScalars().GetName() 221 cdata = "<tr><td><b> cell data array </b></td><td>" + name + "</td></tr>" 222 223 allt = [ 224 "<table>", 225 "<tr>", 226 "<td>", 227 image, 228 "</td>", 229 "<td style='text-align: center; vertical-align: center;'><br/>", 230 help_text, 231 "<table>", 232 "<tr><td><b> bounds </b> <br/> (x/y/z) </td><td>" + str(bounds) + "</td></tr>", 233 "<tr><td><b> center of mass </b></td><td>" 234 + precision(self.center_of_mass(), 3) 235 + "</td></tr>", 236 "<tr><td><b> average size </b></td><td>" + str(average_size) + "</td></tr>", 237 "<tr><td><b> nr. points / faces </b></td><td>" 238 + str(self.npoints) 239 + " / " 240 + str(self.ncells) 241 + "</td></tr>", 242 pdata, 243 cdata, 244 "</table>", 245 "</table>", 246 ] 247 return "\n".join(allt) 248 249 def faces(self, ids=()): 250 """DEPRECATED. Use property `mesh.cells` instead.""" 251 vedo.printc("WARNING: use property mesh.cells instead of mesh.faces()",c='y') 252 return self.cells 253 254 @property 255 def edges(self): 256 """Return an array containing the edges connectivity.""" 257 extractEdges = vtki.new("ExtractEdges") 258 extractEdges.SetInputData(self.dataset) 259 # eed.UseAllPointsOn() 260 extractEdges.Update() 261 lpoly = extractEdges.GetOutput() 262 263 arr1d = vtk2numpy(lpoly.GetLines().GetData()) 264 # [nids1, id0 ... idn, niids2, id0 ... idm, etc]. 265 266 i = 0 267 conn = [] 268 n = len(arr1d) 269 for _ in range(n): 270 cell = [arr1d[i + k + 1] for k in range(arr1d[i])] 271 conn.append(cell) 272 i += arr1d[i] + 1 273 if i >= n: 274 break 275 return conn # cannot always make a numpy array of it! 276 277 @property 278 def cell_normals(self): 279 """ 280 Retrieve face normals as a numpy array. 281 Check out also `compute_normals(cells=True)` and `compute_normals_with_pca()`. 282 """ 283 vtknormals = self.dataset.GetCellData().GetNormals() 284 numpy_normals = vtk2numpy(vtknormals) 285 if len(numpy_normals) == 0 and len(self.cells) != 0: 286 raise ValueError("VTK failed to return any normal vectors. You may need to call `Mesh.compute_normals()` before accessing `Mesh.cell_normals`.") 287 return numpy_normals 288 289 def compute_normals(self, points=True, cells=True, feature_angle=None, consistency=True) -> Self: 290 """ 291 Compute cell and vertex normals for the mesh. 292 293 Arguments: 294 points : (bool) 295 do the computation for the vertices too 296 cells : (bool) 297 do the computation for the cells too 298 feature_angle : (float) 299 specify the angle that defines a sharp edge. 300 If the difference in angle across neighboring polygons is greater than this value, 301 the shared edge is considered "sharp" and it is split. 302 consistency : (bool) 303 turn on/off the enforcement of consistent polygon ordering. 304 305 .. warning:: 306 If `feature_angle` is set then the Mesh can be modified, and it 307 can have a different nr. of vertices from the original. 308 309 Note that the appearance of the mesh may change if the normals are computed, 310 as shading is automatically enabled when such information is present. 311 Use `mesh.flat()` to avoid smoothing effects. 312 """ 313 pdnorm = vtki.new("PolyDataNormals") 314 pdnorm.SetInputData(self.dataset) 315 pdnorm.SetComputePointNormals(points) 316 pdnorm.SetComputeCellNormals(cells) 317 pdnorm.SetConsistency(consistency) 318 pdnorm.FlipNormalsOff() 319 if feature_angle: 320 pdnorm.SetSplitting(True) 321 pdnorm.SetFeatureAngle(feature_angle) 322 else: 323 pdnorm.SetSplitting(False) 324 pdnorm.Update() 325 out = pdnorm.GetOutput() 326 self._update(out, reset_locators=False) 327 return self 328 329 def reverse(self, cells=True, normals=False) -> Self: 330 """ 331 Reverse the order of polygonal cells 332 and/or reverse the direction of point and cell normals. 333 334 Two flags are used to control these operations: 335 - `cells=True` reverses the order of the indices in the cell connectivity list. 336 If cell is a list of IDs only those cells will be reversed. 337 - `normals=True` reverses the normals by multiplying the normal vector by -1 338 (both point and cell normals, if present). 339 """ 340 poly = self.dataset 341 342 if is_sequence(cells): 343 for cell in cells: 344 poly.ReverseCell(cell) 345 poly.GetCellData().Modified() 346 return self ############## 347 348 rev = vtki.new("ReverseSense") 349 if cells: 350 rev.ReverseCellsOn() 351 else: 352 rev.ReverseCellsOff() 353 if normals: 354 rev.ReverseNormalsOn() 355 else: 356 rev.ReverseNormalsOff() 357 rev.SetInputData(poly) 358 rev.Update() 359 self._update(rev.GetOutput(), reset_locators=False) 360 self.pipeline = OperationNode("reverse", parents=[self]) 361 return self 362 363 def volume(self) -> float: 364 """ 365 Compute the volume occupied by mesh. 366 The mesh must be triangular for this to work. 367 To triangulate a mesh use `mesh.triangulate()`. 368 """ 369 mass = vtki.new("MassProperties") 370 mass.SetGlobalWarningDisplay(0) 371 mass.SetInputData(self.dataset) 372 mass.Update() 373 mass.SetGlobalWarningDisplay(1) 374 return mass.GetVolume() 375 376 def area(self) -> float: 377 """ 378 Compute the surface area of the mesh. 379 The mesh must be triangular for this to work. 380 To triangulate a mesh use `mesh.triangulate()`. 381 """ 382 mass = vtki.new("MassProperties") 383 mass.SetGlobalWarningDisplay(0) 384 mass.SetInputData(self.dataset) 385 mass.Update() 386 mass.SetGlobalWarningDisplay(1) 387 return mass.GetSurfaceArea() 388 389 def is_closed(self) -> bool: 390 """ 391 Return `True` if the mesh is watertight. 392 Note that if the mesh contains coincident points the result may be flase. 393 Use in this case `mesh.clean()` to merge coincident points. 394 """ 395 fe = vtki.new("FeatureEdges") 396 fe.BoundaryEdgesOn() 397 fe.FeatureEdgesOff() 398 fe.NonManifoldEdgesOn() 399 fe.SetInputData(self.dataset) 400 fe.Update() 401 ne = fe.GetOutput().GetNumberOfCells() 402 return not bool(ne) 403 404 def is_manifold(self) -> bool: 405 """Return `True` if the mesh is manifold.""" 406 fe = vtki.new("FeatureEdges") 407 fe.BoundaryEdgesOff() 408 fe.FeatureEdgesOff() 409 fe.NonManifoldEdgesOn() 410 fe.SetInputData(self.dataset) 411 fe.Update() 412 ne = fe.GetOutput().GetNumberOfCells() 413 return not bool(ne) 414 415 def non_manifold_faces(self, remove=True, tol="auto") -> Self: 416 """ 417 Detect and (try to) remove non-manifold faces of a triangular mesh: 418 419 - set `remove` to `False` to mark cells without removing them. 420 - set `tol=0` for zero-tolerance, the result will be manifold but with holes. 421 - set `tol>0` to cut off non-manifold faces, and try to recover the good ones. 422 - set `tol="auto"` to make an automatic choice of the tolerance. 423 """ 424 # mark original point and cell ids 425 self.add_ids() 426 toremove = self.boundaries( 427 boundary_edges=False, 428 non_manifold_edges=True, 429 cell_edge=True, 430 return_cell_ids=True, 431 ) 432 if len(toremove) == 0: # type: ignore 433 return self 434 435 points = self.vertices 436 faces = self.cells 437 centers = self.cell_centers 438 439 copy = self.clone() 440 copy.delete_cells(toremove).clean() 441 copy.compute_normals(cells=False) 442 normals = copy.vertex_normals 443 deltas, deltas_i = [], [] 444 445 for i in vedo.utils.progressbar(toremove, delay=3, title="recover faces"): 446 pids = copy.closest_point(centers[i], n=3, return_point_id=True) 447 norms = normals[pids] 448 n = np.mean(norms, axis=0) 449 dn = np.linalg.norm(n) 450 if not dn: 451 continue 452 n = n / dn 453 454 p0, p1, p2 = points[faces[i]][:3] 455 v = np.cross(p1 - p0, p2 - p0) 456 lv = np.linalg.norm(v) 457 if not lv: 458 continue 459 v = v / lv 460 461 cosa = 1 - np.dot(n, v) 462 deltas.append(cosa) 463 deltas_i.append(i) 464 465 recover = [] 466 if len(deltas) > 0: 467 mean_delta = np.mean(deltas) 468 err_delta = np.std(deltas) 469 txt = "" 470 if tol == "auto": # automatic choice 471 tol = mean_delta / 5 472 txt = f"\n Automatic tol. : {tol: .4f}" 473 for i, cosa in zip(deltas_i, deltas): 474 if cosa < tol: 475 recover.append(i) 476 477 vedo.logger.info( 478 f"\n --------- Non manifold faces ---------" 479 f"\n Average tol. : {mean_delta: .4f} +- {err_delta: .4f}{txt}" 480 f"\n Removed faces : {len(toremove)}" # type: ignore 481 f"\n Recovered faces: {len(recover)}" 482 ) 483 484 toremove = list(set(toremove) - set(recover)) # type: ignore 485 486 if not remove: 487 mark = np.zeros(self.ncells, dtype=np.uint8) 488 mark[recover] = 1 489 mark[toremove] = 2 490 self.celldata["NonManifoldCell"] = mark 491 else: 492 self.delete_cells(toremove) # type: ignore 493 494 self.pipeline = OperationNode( 495 "non_manifold_faces", 496 parents=[self], 497 comment=f"#cells {self.dataset.GetNumberOfCells()}", 498 ) 499 return self 500 501 502 def euler_characteristic(self) -> int: 503 """ 504 Compute the Euler characteristic of the mesh. 505 The Euler characteristic is a topological invariant for surfaces. 506 """ 507 return self.npoints - len(self.edges) + self.ncells 508 509 def genus(self) -> int: 510 """ 511 Compute the genus of the mesh. 512 The genus is a topological invariant for surfaces. 513 """ 514 nb = len(self.boundaries().split()) - 1 515 return (2 - self.euler_characteristic() - nb ) / 2 516 517 def to_reeb_graph(self, field_id=0): 518 """ 519 Convert the mesh into a Reeb graph. 520 The Reeb graph is a topological structure that captures the evolution 521 of the level sets of a scalar field. 522 523 Arguments: 524 field_id : (int) 525 the id of the scalar field to use. 526 527 Example: 528 ```python 529 from vedo import * 530 mesh = Mesh("https://discourse.paraview.org/uploads/short-url/qVuZ1fiRjwhE1qYtgGE2HGXybgo.stl") 531 mesh.rotate_x(10).rotate_y(15).alpha(0.5) 532 mesh.pointdata["scalars"] = mesh.vertices[:, 2] 533 534 printc("is_closed :", mesh.is_closed()) 535 printc("is_manifold:", mesh.is_manifold()) 536 printc("euler_char :", mesh.euler_characteristic()) 537 printc("genus :", mesh.genus()) 538 539 reeb = mesh.to_reeb_graph() 540 ids = reeb[0].pointdata["Vertex Ids"] 541 pts = Points(mesh.vertices[ids], r=10) 542 543 show([[mesh, pts], reeb], N=2, sharecam=False) 544 ``` 545 """ 546 rg = vtki.new("PolyDataToReebGraphFilter") 547 rg.SetInputData(self.dataset) 548 rg.SetFieldId(field_id) 549 rg.Update() 550 gr = vedo.pyplot.DirectedGraph() 551 gr.mdg = rg.GetOutput() 552 gr.build() 553 return gr 554 555 556 def shrink(self, fraction=0.85) -> Self: 557 """ 558 Shrink the triangle polydata in the representation of the input mesh. 559 560 Examples: 561 - [shrink.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/shrink.py) 562 563  564 """ 565 # Overriding base class method core.shrink() 566 shrink = vtki.new("ShrinkPolyData") 567 shrink.SetInputData(self.dataset) 568 shrink.SetShrinkFactor(fraction) 569 shrink.Update() 570 self._update(shrink.GetOutput()) 571 self.pipeline = OperationNode("shrink", parents=[self]) 572 return self 573 574 def cap(self, return_cap=False) -> Self: 575 """ 576 Generate a "cap" on a clipped mesh, or caps sharp edges. 577 578 Examples: 579 - [cut_and_cap.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/cut_and_cap.py) 580 581  582 583 See also: `join()`, `join_segments()`, `slice()`. 584 """ 585 fe = vtki.new("FeatureEdges") 586 fe.SetInputData(self.dataset) 587 fe.BoundaryEdgesOn() 588 fe.FeatureEdgesOff() 589 fe.NonManifoldEdgesOff() 590 fe.ManifoldEdgesOff() 591 fe.Update() 592 593 stripper = vtki.new("Stripper") 594 stripper.SetInputData(fe.GetOutput()) 595 stripper.JoinContiguousSegmentsOn() 596 stripper.Update() 597 598 boundary_poly = vtki.vtkPolyData() 599 boundary_poly.SetPoints(stripper.GetOutput().GetPoints()) 600 boundary_poly.SetPolys(stripper.GetOutput().GetLines()) 601 602 rev = vtki.new("ReverseSense") 603 rev.ReverseCellsOn() 604 rev.SetInputData(boundary_poly) 605 rev.Update() 606 607 tf = vtki.new("TriangleFilter") 608 tf.SetInputData(rev.GetOutput()) 609 tf.Update() 610 611 if return_cap: 612 m = Mesh(tf.GetOutput()) 613 m.pipeline = OperationNode( 614 "cap", parents=[self], comment=f"#pts {m.dataset.GetNumberOfPoints()}" 615 ) 616 m.name = "MeshCap" 617 return m 618 619 polyapp = vtki.new("AppendPolyData") 620 polyapp.AddInputData(self.dataset) 621 polyapp.AddInputData(tf.GetOutput()) 622 polyapp.Update() 623 624 self._update(polyapp.GetOutput()) 625 self.clean() 626 627 self.pipeline = OperationNode( 628 "capped", parents=[self], comment=f"#pts {self.dataset.GetNumberOfPoints()}" 629 ) 630 return self 631 632 def join(self, polys=True, reset=False) -> Self: 633 """ 634 Generate triangle strips and/or polylines from 635 input polygons, triangle strips, and lines. 636 637 Input polygons are assembled into triangle strips only if they are triangles; 638 other types of polygons are passed through to the output and not stripped. 639 Use mesh.triangulate() to triangulate non-triangular polygons prior to running 640 this filter if you need to strip all the data. 641 642 Also note that if triangle strips or polylines are present in the input 643 they are passed through and not joined nor extended. 644 If you wish to strip these use mesh.triangulate() to fragment the input 645 into triangles and lines prior to applying join(). 646 647 Arguments: 648 polys : (bool) 649 polygonal segments will be joined if they are contiguous 650 reset : (bool) 651 reset points ordering 652 653 Warning: 654 If triangle strips or polylines exist in the input data 655 they will be passed through to the output data. 656 This filter will only construct triangle strips if triangle polygons 657 are available; and will only construct polylines if lines are available. 658 659 Example: 660 ```python 661 from vedo import * 662 c1 = Cylinder(pos=(0,0,0), r=2, height=3, axis=(1,.0,0), alpha=.1).triangulate() 663 c2 = Cylinder(pos=(0,0,2), r=1, height=2, axis=(0,.3,1), alpha=.1).triangulate() 664 intersect = c1.intersect_with(c2).join(reset=True) 665 spline = Spline(intersect).c('blue').lw(5) 666 show(c1, c2, spline, intersect.labels('id'), axes=1).close() 667 ``` 668  669 """ 670 sf = vtki.new("Stripper") 671 sf.SetPassThroughCellIds(True) 672 sf.SetPassThroughPointIds(True) 673 sf.SetJoinContiguousSegments(polys) 674 sf.SetInputData(self.dataset) 675 sf.Update() 676 if reset: 677 poly = sf.GetOutput() 678 cpd = vtki.new("CleanPolyData") 679 cpd.PointMergingOn() 680 cpd.ConvertLinesToPointsOn() 681 cpd.ConvertPolysToLinesOn() 682 cpd.ConvertStripsToPolysOn() 683 cpd.SetInputData(poly) 684 cpd.Update() 685 poly = cpd.GetOutput() 686 vpts = poly.GetCell(0).GetPoints().GetData() 687 poly.GetPoints().SetData(vpts) 688 else: 689 poly = sf.GetOutput() 690 691 self._update(poly) 692 693 self.pipeline = OperationNode( 694 "join", parents=[self], comment=f"#pts {self.dataset.GetNumberOfPoints()}" 695 ) 696 return self 697 698 def join_segments(self, closed=True, tol=1e-03) -> list: 699 """ 700 Join line segments into contiguous lines. 701 Useful to call with `triangulate()` method. 702 703 Returns: 704 list of `shapes.Lines` 705 706 Example: 707 ```python 708 from vedo import * 709 msh = Torus().alpha(0.1).wireframe() 710 intersection = msh.intersect_with_plane(normal=[1,1,1]).c('purple5') 711 slices = [s.triangulate() for s in intersection.join_segments()] 712 show(msh, intersection, merge(slices), axes=1, viewup='z') 713 ``` 714  715 """ 716 vlines = [] 717 for ipiece, outline in enumerate(self.split(must_share_edge=False)): # type: ignore 718 719 outline.clean() 720 pts = outline.vertices 721 if len(pts) < 3: 722 continue 723 avesize = outline.average_size() 724 lines = outline.lines 725 # print("---lines", lines, "in piece", ipiece) 726 tol = avesize / pts.shape[0] * tol 727 728 k = 0 729 joinedpts = [pts[k]] 730 for _ in range(len(pts)): 731 pk = pts[k] 732 for j, line in enumerate(lines): 733 734 id0, id1 = line[0], line[-1] 735 p0, p1 = pts[id0], pts[id1] 736 737 if np.linalg.norm(p0 - pk) < tol: 738 n = len(line) 739 for m in range(1, n): 740 joinedpts.append(pts[line[m]]) 741 # joinedpts.append(p1) 742 k = id1 743 lines.pop(j) 744 break 745 746 elif np.linalg.norm(p1 - pk) < tol: 747 n = len(line) 748 for m in reversed(range(0, n - 1)): 749 joinedpts.append(pts[line[m]]) 750 # joinedpts.append(p0) 751 k = id0 752 lines.pop(j) 753 break 754 755 if len(joinedpts) > 1: 756 newline = vedo.shapes.Line(joinedpts, closed=closed) 757 newline.clean() 758 newline.actor.SetProperty(self.properties) 759 newline.properties = self.properties 760 newline.pipeline = OperationNode( 761 "join_segments", 762 parents=[self], 763 comment=f"#pts {newline.dataset.GetNumberOfPoints()}", 764 ) 765 vlines.append(newline) 766 767 return vlines 768 769 def join_with_strips(self, b1, closed=True) -> Self: 770 """ 771 Join booundary lines by creating a triangle strip between them. 772 773 Example: 774 ```python 775 from vedo import * 776 m1 = Cylinder(cap=False).boundaries() 777 m2 = Cylinder(cap=False).boundaries().pos(0.2,0,1) 778 strips = m1.join_with_strips(m2) 779 show(m1, m2, strips, axes=1).close() 780 ``` 781 """ 782 b0 = self.clone().join() 783 b1 = b1.clone().join() 784 785 vertices0 = b0.vertices.tolist() 786 vertices1 = b1.vertices.tolist() 787 788 lines0 = b0.lines 789 lines1 = b1.lines 790 m = len(lines0) 791 assert m == len(lines1), ( 792 "lines must have the same number of points\n" 793 f"line has {m} points in b0 and {len(lines1)} in b1" 794 ) 795 796 strips = [] 797 points: List[Any] = [] 798 799 for j in range(m): 800 801 ids0j = list(lines0[j]) 802 ids1j = list(lines1[j]) 803 804 n = len(ids0j) 805 assert n == len(ids1j), ( 806 "lines must have the same number of points\n" 807 f"line {j} has {n} points in b0 and {len(ids1j)} in b1" 808 ) 809 810 if closed: 811 ids0j.append(ids0j[0]) 812 ids1j.append(ids1j[0]) 813 vertices0.append(vertices0[ids0j[0]]) 814 vertices1.append(vertices1[ids1j[0]]) 815 n = n + 1 816 817 strip = [] # create a triangle strip 818 npt = len(points) 819 for ipt in range(n): 820 points.append(vertices0[ids0j[ipt]]) 821 points.append(vertices1[ids1j[ipt]]) 822 823 strip = list(range(npt, npt + 2*n)) 824 strips.append(strip) 825 826 return Mesh([points, [], [], strips], c="k6") 827 828 def split_polylines(self) -> Self: 829 """Split polylines into separate segments.""" 830 tf = vtki.new("TriangleFilter") 831 tf.SetPassLines(True) 832 tf.SetPassVerts(False) 833 tf.SetInputData(self.dataset) 834 tf.Update() 835 self._update(tf.GetOutput(), reset_locators=False) 836 self.lw(0).lighting("default").pickable() 837 self.pipeline = OperationNode( 838 "split_polylines", parents=[self], 839 comment=f"#lines {self.dataset.GetNumberOfLines()}" 840 ) 841 return self 842 843 def slice(self, origin=(0, 0, 0), normal=(1, 0, 0)) -> Self: 844 """ 845 Slice a mesh with a plane and fill the contour. 846 847 Example: 848 ```python 849 from vedo import * 850 msh = Mesh(dataurl+"bunny.obj").alpha(0.1).wireframe() 851 mslice = msh.slice(normal=[0,1,0.3], origin=[0,0.16,0]) 852 mslice.c('purple5') 853 show(msh, mslice, axes=1) 854 ``` 855  856 857 See also: `join()`, `join_segments()`, `cap()`, `cut_with_plane()`. 858 """ 859 intersection = self.intersect_with_plane(origin=origin, normal=normal) 860 slices = [s.triangulate() for s in intersection.join_segments()] 861 mslices = vedo.pointcloud.merge(slices) 862 if mslices: 863 mslices.name = "MeshSlice" 864 mslices.pipeline = OperationNode("slice", parents=[self], comment=f"normal = {normal}") 865 return mslices 866 867 def triangulate(self, verts=True, lines=True) -> Self: 868 """ 869 Converts mesh polygons into triangles. 870 871 If the input mesh is only made of 2D lines (no faces) the output will be a triangulation 872 that fills the internal area. The contours may be concave, and may even contain holes, 873 i.e. a contour may contain an internal contour winding in the opposite 874 direction to indicate that it is a hole. 875 876 Arguments: 877 verts : (bool) 878 if True, break input vertex cells into individual vertex cells (one point per cell). 879 If False, the input vertex cells will be ignored. 880 lines : (bool) 881 if True, break input polylines into line segments. 882 If False, input lines will be ignored and the output will have no lines. 883 """ 884 if self.dataset.GetNumberOfPolys() or self.dataset.GetNumberOfStrips(): 885 # print("Using vtkTriangleFilter") 886 tf = vtki.new("TriangleFilter") 887 tf.SetPassLines(lines) 888 tf.SetPassVerts(verts) 889 890 elif self.dataset.GetNumberOfLines(): 891 # print("Using vtkContourTriangulator") 892 tf = vtki.new("ContourTriangulator") 893 tf.TriangulationErrorDisplayOn() 894 895 else: 896 vedo.logger.debug("input in triangulate() seems to be void! Skip.") 897 return self 898 899 tf.SetInputData(self.dataset) 900 tf.Update() 901 self._update(tf.GetOutput(), reset_locators=False) 902 self.lw(0).lighting("default").pickable() 903 904 self.pipeline = OperationNode( 905 "triangulate", parents=[self], comment=f"#cells {self.dataset.GetNumberOfCells()}" 906 ) 907 return self 908 909 def compute_cell_vertex_count(self) -> Self: 910 """ 911 Add to this mesh a cell data array containing the nr of vertices that a polygonal face has. 912 """ 913 csf = vtki.new("CellSizeFilter") 914 csf.SetInputData(self.dataset) 915 csf.SetComputeArea(False) 916 csf.SetComputeVolume(False) 917 csf.SetComputeLength(False) 918 csf.SetComputeVertexCount(True) 919 csf.SetVertexCountArrayName("VertexCount") 920 csf.Update() 921 self.dataset.GetCellData().AddArray( 922 csf.GetOutput().GetCellData().GetArray("VertexCount") 923 ) 924 return self 925 926 def compute_quality(self, metric=6) -> Self: 927 """ 928 Calculate metrics of quality for the elements of a triangular mesh. 929 This method adds to the mesh a cell array named "Quality". 930 See class 931 [vtkMeshQuality](https://vtk.org/doc/nightly/html/classvtkMeshQuality.html). 932 933 Arguments: 934 metric : (int) 935 type of available estimators are: 936 - EDGE RATIO, 0 937 - ASPECT RATIO, 1 938 - RADIUS RATIO, 2 939 - ASPECT FROBENIUS, 3 940 - MED ASPECT FROBENIUS, 4 941 - MAX ASPECT FROBENIUS, 5 942 - MIN_ANGLE, 6 943 - COLLAPSE RATIO, 7 944 - MAX ANGLE, 8 945 - CONDITION, 9 946 - SCALED JACOBIAN, 10 947 - SHEAR, 11 948 - RELATIVE SIZE SQUARED, 12 949 - SHAPE, 13 950 - SHAPE AND SIZE, 14 951 - DISTORTION, 15 952 - MAX EDGE RATIO, 16 953 - SKEW, 17 954 - TAPER, 18 955 - VOLUME, 19 956 - STRETCH, 20 957 - DIAGONAL, 21 958 - DIMENSION, 22 959 - ODDY, 23 960 - SHEAR AND SIZE, 24 961 - JACOBIAN, 25 962 - WARPAGE, 26 963 - ASPECT GAMMA, 27 964 - AREA, 28 965 - ASPECT BETA, 29 966 967 Examples: 968 - [meshquality.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/meshquality.py) 969 970  971 """ 972 qf = vtki.new("MeshQuality") 973 qf.SetInputData(self.dataset) 974 qf.SetTriangleQualityMeasure(metric) 975 qf.SaveCellQualityOn() 976 qf.Update() 977 self._update(qf.GetOutput(), reset_locators=False) 978 self.mapper.SetScalarModeToUseCellData() 979 self.pipeline = OperationNode("compute_quality", parents=[self]) 980 return self 981 982 def count_vertices(self) -> np.ndarray: 983 """Count the number of vertices each cell has and return it as a numpy array""" 984 vc = vtki.new("CountVertices") 985 vc.SetInputData(self.dataset) 986 vc.SetOutputArrayName("VertexCount") 987 vc.Update() 988 varr = vc.GetOutput().GetCellData().GetArray("VertexCount") 989 return vtk2numpy(varr) 990 991 def check_validity(self, tol=0) -> np.ndarray: 992 """ 993 Return a numpy array of possible problematic faces following this convention: 994 - Valid = 0 995 - WrongNumberOfPoints = 1 996 - IntersectingEdges = 2 997 - IntersectingFaces = 4 998 - NoncontiguousEdges = 8 999 - Nonconvex = 10 1000 - OrientedIncorrectly = 20 1001 1002 Arguments: 1003 tol : (float) 1004 value is used as an epsilon for floating point 1005 equality checks throughout the cell checking process. 1006 """ 1007 vald = vtki.new("CellValidator") 1008 if tol: 1009 vald.SetTolerance(tol) 1010 vald.SetInputData(self.dataset) 1011 vald.Update() 1012 varr = vald.GetOutput().GetCellData().GetArray("ValidityState") 1013 return vtk2numpy(varr) 1014 1015 def compute_curvature(self, method=0) -> Self: 1016 """ 1017 Add scalars to `Mesh` that contains the curvature calculated in three different ways. 1018 1019 Variable `method` can be: 1020 - 0 = gaussian 1021 - 1 = mean curvature 1022 - 2 = max curvature 1023 - 3 = min curvature 1024 1025 Example: 1026 ```python 1027 from vedo import Torus 1028 Torus().compute_curvature().add_scalarbar().show().close() 1029 ``` 1030  1031 """ 1032 curve = vtki.new("Curvatures") 1033 curve.SetInputData(self.dataset) 1034 curve.SetCurvatureType(method) 1035 curve.Update() 1036 self._update(curve.GetOutput(), reset_locators=False) 1037 self.mapper.ScalarVisibilityOn() 1038 return self 1039 1040 def compute_elevation(self, low=(0, 0, 0), high=(0, 0, 1), vrange=(0, 1)) -> Self: 1041 """ 1042 Add to `Mesh` a scalar array that contains distance along a specified direction. 1043 1044 Arguments: 1045 low : (list) 1046 one end of the line (small scalar values) 1047 high : (list) 1048 other end of the line (large scalar values) 1049 vrange : (list) 1050 set the range of the scalar 1051 1052 Example: 1053 ```python 1054 from vedo import Sphere 1055 s = Sphere().compute_elevation(low=(0,0,0), high=(1,1,1)) 1056 s.add_scalarbar().show(axes=1).close() 1057 ``` 1058  1059 """ 1060 ef = vtki.new("ElevationFilter") 1061 ef.SetInputData(self.dataset) 1062 ef.SetLowPoint(low) 1063 ef.SetHighPoint(high) 1064 ef.SetScalarRange(vrange) 1065 ef.Update() 1066 self._update(ef.GetOutput(), reset_locators=False) 1067 self.mapper.ScalarVisibilityOn() 1068 return self 1069 1070 1071 def laplacian_diffusion(self, array_name, dt, num_steps) -> Self: 1072 """ 1073 Apply a diffusion process to a scalar array defined on the points of a mesh. 1074 1075 Arguments: 1076 array_name : (str) 1077 name of the array to diffuse. 1078 dt : (float) 1079 time step. 1080 num_steps : (int) 1081 number of iterations. 1082 """ 1083 try: 1084 import scipy.sparse 1085 import scipy.sparse.linalg 1086 except ImportError: 1087 vedo.logger.error("scipy not found. Cannot run laplacian_diffusion()") 1088 return self 1089 1090 def build_laplacian(): 1091 rows = [] 1092 cols = [] 1093 data = [] 1094 n_points = points.shape[0] 1095 avg_area = np.mean(areas) * 10000 1096 # print("avg_area", avg_area) 1097 1098 for triangle in cells: 1099 for i in range(3): 1100 for j in range(i + 1, 3): 1101 u = triangle[i] 1102 v = triangle[j] 1103 rows.append(u) 1104 cols.append(v) 1105 rows.append(v) 1106 cols.append(u) 1107 data.append(-1/avg_area) 1108 data.append(-1/avg_area) 1109 1110 L = scipy.sparse.coo_matrix( 1111 (data, (rows, cols)), shape=(n_points, n_points) 1112 ).tocsc() 1113 1114 degree = -np.array(L.sum(axis=1)).flatten() # adjust the diagonal 1115 # print("degree", degree) 1116 L.setdiag(degree) 1117 return L 1118 1119 def _diffuse(u0, L, dt, num_steps): 1120 # mean_area = np.mean(areas) * 10000 1121 # print("mean_area", mean_area) 1122 mean_area = 1 1123 I = scipy.sparse.eye(L.shape[0], format="csc") 1124 A = I - (dt/mean_area) * L 1125 u = u0 1126 for _ in range(int(num_steps)): 1127 u = A.dot(u) 1128 return u 1129 1130 self.compute_cell_size() 1131 areas = self.celldata["Area"] 1132 points = self.vertices 1133 cells = self.cells 1134 u0 = self.pointdata[array_name] 1135 1136 # Simulate diffusion 1137 L = build_laplacian() 1138 u = _diffuse(u0, L, dt, num_steps) 1139 self.pointdata[array_name] = u 1140 return self 1141 1142 1143 def subdivide(self, n=1, method=0, mel=None) -> Self: 1144 """ 1145 Increase the number of vertices of a surface mesh. 1146 1147 Arguments: 1148 n : (int) 1149 number of subdivisions. 1150 method : (int) 1151 Loop(0), Linear(1), Adaptive(2), Butterfly(3), Centroid(4) 1152 mel : (float) 1153 Maximum Edge Length (applicable to Adaptive method only). 1154 """ 1155 triangles = vtki.new("TriangleFilter") 1156 triangles.SetInputData(self.dataset) 1157 triangles.Update() 1158 tri_mesh = triangles.GetOutput() 1159 if method == 0: 1160 sdf = vtki.new("LoopSubdivisionFilter") 1161 elif method == 1: 1162 sdf = vtki.new("LinearSubdivisionFilter") 1163 elif method == 2: 1164 sdf = vtki.new("AdaptiveSubdivisionFilter") 1165 if mel is None: 1166 mel = self.diagonal_size() / np.sqrt(self.dataset.GetNumberOfPoints()) / n 1167 sdf.SetMaximumEdgeLength(mel) 1168 elif method == 3: 1169 sdf = vtki.new("ButterflySubdivisionFilter") 1170 elif method == 4: 1171 sdf = vtki.new("DensifyPolyData") 1172 else: 1173 vedo.logger.error(f"in subdivide() unknown method {method}") 1174 raise RuntimeError() 1175 1176 if method != 2: 1177 sdf.SetNumberOfSubdivisions(n) 1178 1179 sdf.SetInputData(tri_mesh) 1180 sdf.Update() 1181 1182 self._update(sdf.GetOutput()) 1183 1184 self.pipeline = OperationNode( 1185 "subdivide", 1186 parents=[self], 1187 comment=f"#pts {self.dataset.GetNumberOfPoints()}", 1188 ) 1189 return self 1190 1191 1192 def decimate(self, fraction=0.5, n=None, preserve_volume=True, regularization=0.0) -> Self: 1193 """ 1194 Downsample the number of vertices in a mesh to `fraction`. 1195 1196 This filter preserves the `pointdata` of the input dataset. In previous versions 1197 of vedo, this decimation algorithm was referred to as quadric decimation. 1198 1199 Arguments: 1200 fraction : (float) 1201 the desired target of reduction. 1202 n : (int) 1203 the desired number of final points 1204 (`fraction` is recalculated based on it). 1205 preserve_volume : (bool) 1206 Decide whether to activate volume preservation which greatly 1207 reduces errors in triangle normal direction. 1208 regularization : (float) 1209 regularize the point finding algorithm so as to have better quality 1210 mesh elements at the cost of a slightly lower precision on the 1211 geometry potentially (mostly at sharp edges). 1212 Can be useful for decimating meshes that have been triangulated on noisy data. 1213 1214 Note: 1215 Setting `fraction=0.1` leaves 10% of the original number of vertices. 1216 Internally the VTK class 1217 [vtkQuadricDecimation](https://vtk.org/doc/nightly/html/classvtkQuadricDecimation.html) 1218 is used for this operation. 1219 1220 See also: `decimate_binned()` and `decimate_pro()`. 1221 """ 1222 poly = self.dataset 1223 if n: # N = desired number of points 1224 npt = poly.GetNumberOfPoints() 1225 fraction = n / npt 1226 if fraction >= 1: 1227 return self 1228 1229 decimate = vtki.new("QuadricDecimation") 1230 decimate.SetVolumePreservation(preserve_volume) 1231 # decimate.AttributeErrorMetricOn() 1232 if regularization: 1233 decimate.SetRegularize(True) 1234 decimate.SetRegularization(regularization) 1235 1236 try: 1237 decimate.MapPointDataOn() 1238 except AttributeError: 1239 pass 1240 1241 decimate.SetTargetReduction(1 - fraction) 1242 decimate.SetInputData(poly) 1243 decimate.Update() 1244 1245 self._update(decimate.GetOutput()) 1246 self.metadata["decimate_actual_fraction"] = 1 - decimate.GetActualReduction() 1247 1248 self.pipeline = OperationNode( 1249 "decimate", 1250 parents=[self], 1251 comment=f"#pts {self.dataset.GetNumberOfPoints()}", 1252 ) 1253 return self 1254 1255 def decimate_pro( 1256 self, 1257 fraction=0.5, 1258 n=None, 1259 preserve_topology=True, 1260 preserve_boundaries=True, 1261 splitting=False, 1262 splitting_angle=75, 1263 feature_angle=0, 1264 inflection_point_ratio=10, 1265 vertex_degree=0, 1266 ) -> Self: 1267 """ 1268 Downsample the number of vertices in a mesh to `fraction`. 1269 1270 This filter preserves the `pointdata` of the input dataset. 1271 1272 Arguments: 1273 fraction : (float) 1274 The desired target of reduction. 1275 Setting `fraction=0.1` leaves 10% of the original number of vertices. 1276 n : (int) 1277 the desired number of final points (`fraction` is recalculated based on it). 1278 preserve_topology : (bool) 1279 If on, mesh splitting and hole elimination will not occur. 1280 This may limit the maximum reduction that may be achieved. 1281 preserve_boundaries : (bool) 1282 Turn on/off the deletion of vertices on the boundary of a mesh. 1283 Control whether mesh boundaries are preserved during decimation. 1284 feature_angle : (float) 1285 Specify the angle that defines a feature. 1286 This angle is used to define what an edge is 1287 (i.e., if the surface normal between two adjacent triangles 1288 is >= FeatureAngle, an edge exists). 1289 splitting : (bool) 1290 Turn on/off the splitting of the mesh at corners, 1291 along edges, at non-manifold points, or anywhere else a split is required. 1292 Turning splitting off will better preserve the original topology of the mesh, 1293 but you may not obtain the requested reduction. 1294 splitting_angle : (float) 1295 Specify the angle that defines a sharp edge. 1296 This angle is used to control the splitting of the mesh. 1297 A split line exists when the surface normals between two edge connected triangles 1298 are >= `splitting_angle`. 1299 inflection_point_ratio : (float) 1300 An inflection point occurs when the ratio of reduction error between two iterations 1301 is greater than or equal to the `inflection_point_ratio` value. 1302 vertex_degree : (int) 1303 If the number of triangles connected to a vertex exceeds it then the vertex will be split. 1304 1305 Note: 1306 Setting `fraction=0.1` leaves 10% of the original number of vertices 1307 1308 See also: 1309 `decimate()` and `decimate_binned()`. 1310 """ 1311 poly = self.dataset 1312 if n: # N = desired number of points 1313 npt = poly.GetNumberOfPoints() 1314 fraction = n / npt 1315 if fraction >= 1: 1316 return self 1317 1318 decimate = vtki.new("DecimatePro") 1319 decimate.SetPreserveTopology(preserve_topology) 1320 decimate.SetBoundaryVertexDeletion(preserve_boundaries) 1321 if feature_angle: 1322 decimate.SetFeatureAngle(feature_angle) 1323 decimate.SetSplitting(splitting) 1324 decimate.SetSplitAngle(splitting_angle) 1325 decimate.SetInflectionPointRatio(inflection_point_ratio) 1326 if vertex_degree: 1327 decimate.SetDegree(vertex_degree) 1328 1329 decimate.SetTargetReduction(1 - fraction) 1330 decimate.SetInputData(poly) 1331 decimate.Update() 1332 self._update(decimate.GetOutput()) 1333 1334 self.pipeline = OperationNode( 1335 "decimate_pro", 1336 parents=[self], 1337 comment=f"#pts {self.dataset.GetNumberOfPoints()}", 1338 ) 1339 return self 1340 1341 def decimate_binned(self, divisions=(), use_clustering=False) -> Self: 1342 """ 1343 Downsample the number of vertices in a mesh. 1344 1345 This filter preserves the `celldata` of the input dataset, 1346 if `use_clustering=True` also the `pointdata` will be preserved in the result. 1347 1348 Arguments: 1349 divisions : (list) 1350 number of divisions along x, y and z axes. 1351 auto_adjust : (bool) 1352 if True, the number of divisions is automatically adjusted to 1353 create more uniform cells. 1354 use_clustering : (bool) 1355 use [vtkQuadricClustering](https://vtk.org/doc/nightly/html/classvtkQuadricClustering.html) 1356 instead of 1357 [vtkBinnedDecimation](https://vtk.org/doc/nightly/html/classvtkBinnedDecimation.html). 1358 1359 See also: `decimate()` and `decimate_pro()`. 1360 """ 1361 if use_clustering: 1362 decimate = vtki.new("QuadricClustering") 1363 decimate.CopyCellDataOn() 1364 else: 1365 decimate = vtki.new("BinnedDecimation") 1366 decimate.ProducePointDataOn() 1367 decimate.ProduceCellDataOn() 1368 1369 decimate.SetInputData(self.dataset) 1370 1371 if len(divisions) == 0: 1372 decimate.SetAutoAdjustNumberOfDivisions(1) 1373 else: 1374 decimate.SetAutoAdjustNumberOfDivisions(0) 1375 decimate.SetNumberOfDivisions(divisions) 1376 decimate.Update() 1377 1378 self._update(decimate.GetOutput()) 1379 self.metadata["decimate_binned_divisions"] = decimate.GetNumberOfDivisions() 1380 self.pipeline = OperationNode( 1381 "decimate_binned", 1382 parents=[self], 1383 comment=f"#pts {self.dataset.GetNumberOfPoints()}", 1384 ) 1385 return self 1386 1387 def generate_random_points(self, n: int, min_radius=0.0) -> "Points": 1388 """ 1389 Generate `n` uniformly distributed random points 1390 inside the polygonal mesh. 1391 1392 A new point data array is added to the output points 1393 called "OriginalCellID" which contains the index of 1394 the cell ID in which the point was generated. 1395 1396 Arguments: 1397 n : (int) 1398 number of points to generate. 1399 min_radius: (float) 1400 impose a minimum distance between points. 1401 If `min_radius` is set to 0, the points are 1402 generated uniformly at random inside the mesh. 1403 If `min_radius` is set to a positive value, 1404 the points are generated uniformly at random 1405 inside the mesh, but points closer than `min_radius` 1406 to any other point are discarded. 1407 1408 Returns a `vedo.Points` object. 1409 1410 Note: 1411 Consider using `points.probe(msh)` or 1412 `points.interpolate_data_from(msh)` 1413 to interpolate existing mesh data onto the new points. 1414 1415 Example: 1416 ```python 1417 from vedo import * 1418 msh = Mesh(dataurl + "panther.stl").lw(2) 1419 pts = msh.generate_random_points(20000, min_radius=0.5) 1420 print("Original cell ids:", pts.pointdata["OriginalCellID"]) 1421 show(pts, msh, axes=1).close() 1422 ``` 1423 """ 1424 cmesh = self.clone().clean().triangulate().compute_cell_size() 1425 triangles = cmesh.cells 1426 vertices = cmesh.vertices 1427 cumul = np.cumsum(cmesh.celldata["Area"]) 1428 1429 out_pts = [] 1430 orig_cell = [] 1431 for _ in range(n): 1432 # choose a triangle based on area 1433 random_area = np.random.random() * cumul[-1] 1434 it = np.searchsorted(cumul, random_area) 1435 A, B, C = vertices[triangles[it]] 1436 # calculate the random point in the triangle 1437 r1, r2 = np.random.random(2) 1438 if r1 + r2 > 1: 1439 r1 = 1 - r1 1440 r2 = 1 - r2 1441 out_pts.append((1 - r1 - r2) * A + r1 * B + r2 * C) 1442 orig_cell.append(it) 1443 nporig_cell = np.array(orig_cell, dtype=np.uint32) 1444 1445 vpts = Points(out_pts) 1446 vpts.pointdata["OriginalCellID"] = nporig_cell 1447 1448 if min_radius > 0: 1449 vpts.subsample(min_radius, absolute=True) 1450 1451 vpts.point_size(5).color("k1") 1452 vpts.name = "RandomPoints" 1453 vpts.pipeline = OperationNode( 1454 "generate_random_points", c="#edabab", parents=[self]) 1455 return vpts 1456 1457 def delete_cells(self, ids: List[int]) -> Self: 1458 """ 1459 Remove cells from the mesh object by their ID. 1460 Points (vertices) are not removed (you may use `clean()` to remove those). 1461 """ 1462 self.dataset.BuildLinks() 1463 for cid in ids: 1464 self.dataset.DeleteCell(cid) 1465 self.dataset.RemoveDeletedCells() 1466 self.dataset.Modified() 1467 self.mapper.Modified() 1468 self.pipeline = OperationNode( 1469 "delete_cells", 1470 parents=[self], 1471 comment=f"#cells {self.dataset.GetNumberOfCells()}", 1472 ) 1473 return self 1474 1475 def delete_cells_by_point_index(self, indices: List[int]) -> Self: 1476 """ 1477 Delete a list of vertices identified by any of their vertex index. 1478 1479 See also `delete_cells()`. 1480 1481 Examples: 1482 - [delete_mesh_pts.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/delete_mesh_pts.py) 1483 1484  1485 """ 1486 cell_ids = vtki.vtkIdList() 1487 self.dataset.BuildLinks() 1488 n = 0 1489 for i in np.unique(indices): 1490 self.dataset.GetPointCells(i, cell_ids) 1491 for j in range(cell_ids.GetNumberOfIds()): 1492 self.dataset.DeleteCell(cell_ids.GetId(j)) # flag cell 1493 n += 1 1494 1495 self.dataset.RemoveDeletedCells() 1496 self.dataset.Modified() 1497 self.pipeline = OperationNode("delete_cells_by_point_index", parents=[self]) 1498 return self 1499 1500 def collapse_edges(self, distance: float, iterations=1) -> Self: 1501 """ 1502 Collapse mesh edges so that are all above `distance`. 1503 1504 Example: 1505 ```python 1506 from vedo import * 1507 np.random.seed(2) 1508 grid1 = Grid().add_gaussian_noise(0.8).triangulate().lw(1) 1509 grid1.celldata['scalar'] = grid1.cell_centers[:,1] 1510 grid2 = grid1.clone().collapse_edges(0.1) 1511 show(grid1, grid2, N=2, axes=1) 1512 ``` 1513 """ 1514 for _ in range(iterations): 1515 medges = self.edges 1516 pts = self.vertices 1517 newpts = np.array(pts) 1518 moved = [] 1519 for e in medges: 1520 if len(e) == 2: 1521 id0, id1 = e 1522 p0, p1 = pts[id0], pts[id1] 1523 if (np.linalg.norm(p1-p0) < distance 1524 and id0 not in moved 1525 and id1 not in moved 1526 ): 1527 p = (p0 + p1) / 2 1528 newpts[id0] = p 1529 newpts[id1] = p 1530 moved += [id0, id1] 1531 self.vertices = newpts 1532 cpd = vtki.new("CleanPolyData") 1533 cpd.ConvertLinesToPointsOff() 1534 cpd.ConvertPolysToLinesOff() 1535 cpd.ConvertStripsToPolysOff() 1536 cpd.SetInputData(self.dataset) 1537 cpd.Update() 1538 self._update(cpd.GetOutput()) 1539 1540 self.pipeline = OperationNode( 1541 "collapse_edges", 1542 parents=[self], 1543 comment=f"#pts {self.dataset.GetNumberOfPoints()}", 1544 ) 1545 return self 1546 1547 def adjacency_list(self) -> List[set]: 1548 """ 1549 Computes the adjacency list for mesh edge-graph. 1550 1551 Returns: 1552 a list with i-th entry being the set if indices of vertices connected by an edge to i-th vertex 1553 """ 1554 inc = [set()] * self.nvertices 1555 for cell in self.cells: 1556 nc = len(cell) 1557 if nc > 1: 1558 for i in range(nc-1): 1559 ci = cell[i] 1560 inc[ci] = inc[ci].union({cell[i-1], cell[i+1]}) 1561 return inc 1562 1563 def graph_ball(self, index, n: int) -> set: 1564 """ 1565 Computes the ball of radius `n` in the mesh' edge-graph metric centred in vertex `index`. 1566 1567 Arguments: 1568 index : (int) 1569 index of the vertex 1570 n : (int) 1571 radius in the graph metric 1572 1573 Returns: 1574 the set of indices of the vertices which are at most `n` edges from vertex `index`. 1575 """ 1576 if n == 0: 1577 return {index} 1578 else: 1579 al = self.adjacency_list() 1580 ball = {index} 1581 i = 0 1582 while i < n and len(ball) < self.nvertices: 1583 for v in ball: 1584 ball = ball.union(al[v]) 1585 i += 1 1586 return ball 1587 1588 def smooth(self, niter=15, pass_band=0.1, edge_angle=15, feature_angle=60, boundary=False) -> Self: 1589 """ 1590 Adjust mesh point positions using the so-called "Windowed Sinc" method. 1591 1592 Arguments: 1593 niter : (int) 1594 number of iterations. 1595 pass_band : (float) 1596 set the pass_band value for the windowed sinc filter. 1597 edge_angle : (float) 1598 edge angle to control smoothing along edges (either interior or boundary). 1599 feature_angle : (float) 1600 specifies the feature angle for sharp edge identification. 1601 boundary : (bool) 1602 specify if boundary should also be smoothed or kept unmodified 1603 1604 Examples: 1605 - [mesh_smoother1.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/mesh_smoother1.py) 1606 1607  1608 """ 1609 cl = vtki.new("CleanPolyData") 1610 cl.SetInputData(self.dataset) 1611 cl.Update() 1612 smf = vtki.new("WindowedSincPolyDataFilter") 1613 smf.SetInputData(cl.GetOutput()) 1614 smf.SetNumberOfIterations(niter) 1615 smf.SetEdgeAngle(edge_angle) 1616 smf.SetFeatureAngle(feature_angle) 1617 smf.SetPassBand(pass_band) 1618 smf.NormalizeCoordinatesOn() 1619 smf.NonManifoldSmoothingOn() 1620 smf.FeatureEdgeSmoothingOn() 1621 smf.SetBoundarySmoothing(boundary) 1622 smf.Update() 1623 1624 self._update(smf.GetOutput()) 1625 1626 self.pipeline = OperationNode( 1627 "smooth", parents=[self], comment=f"#pts {self.dataset.GetNumberOfPoints()}" 1628 ) 1629 return self 1630 1631 def fill_holes(self, size=None) -> Self: 1632 """ 1633 Identifies and fills holes in the input mesh. 1634 Holes are identified by locating boundary edges, linking them together 1635 into loops, and then triangulating the resulting loops. 1636 1637 Arguments: 1638 size : (float) 1639 Approximate limit to the size of the hole that can be filled. 1640 1641 Examples: 1642 - [fillholes.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/fillholes.py) 1643 """ 1644 fh = vtki.new("FillHolesFilter") 1645 if not size: 1646 mb = self.diagonal_size() 1647 size = mb / 10 1648 fh.SetHoleSize(size) 1649 fh.SetInputData(self.dataset) 1650 fh.Update() 1651 1652 self._update(fh.GetOutput()) 1653 1654 self.pipeline = OperationNode( 1655 "fill_holes", 1656 parents=[self], 1657 comment=f"#pts {self.dataset.GetNumberOfPoints()}", 1658 ) 1659 return self 1660 1661 def contains(self, point: tuple, tol=1e-05) -> bool: 1662 """ 1663 Return True if point is inside a polydata closed surface. 1664 1665 Note: 1666 if you have many points to check use `inside_points()` instead. 1667 1668 Example: 1669 ```python 1670 from vedo import * 1671 s = Sphere().c('green5').alpha(0.5) 1672 pt = [0.1, 0.2, 0.3] 1673 print("Sphere contains", pt, s.contains(pt)) 1674 show(s, Point(pt), axes=1).close() 1675 ``` 1676 """ 1677 points = vtki.vtkPoints() 1678 points.InsertNextPoint(point) 1679 poly = vtki.vtkPolyData() 1680 poly.SetPoints(points) 1681 sep = vtki.new("SelectEnclosedPoints") 1682 sep.SetTolerance(tol) 1683 sep.CheckSurfaceOff() 1684 sep.SetInputData(poly) 1685 sep.SetSurfaceData(self.dataset) 1686 sep.Update() 1687 return bool(sep.IsInside(0)) 1688 1689 def inside_points(self, pts: Union["Points", list], invert=False, tol=1e-05, return_ids=False) -> Union["Points", np.ndarray]: 1690 """ 1691 Return the point cloud that is inside mesh surface as a new Points object. 1692 1693 If return_ids is True a list of IDs is returned and in addition input points 1694 are marked by a pointdata array named "IsInside". 1695 1696 Example: 1697 `print(pts.pointdata["IsInside"])` 1698 1699 Examples: 1700 - [pca_ellipsoid.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/pca_ellipsoid.py) 1701 1702  1703 """ 1704 if isinstance(pts, Points): 1705 poly = pts.dataset 1706 ptsa = pts.vertices 1707 else: 1708 ptsa = np.asarray(pts) 1709 vpoints = vtki.vtkPoints() 1710 vpoints.SetData(numpy2vtk(ptsa, dtype=np.float32)) 1711 poly = vtki.vtkPolyData() 1712 poly.SetPoints(vpoints) 1713 1714 sep = vtki.new("SelectEnclosedPoints") 1715 # sep = vtki.new("ExtractEnclosedPoints() 1716 sep.SetTolerance(tol) 1717 sep.SetInputData(poly) 1718 sep.SetSurfaceData(self.dataset) 1719 sep.SetInsideOut(invert) 1720 sep.Update() 1721 1722 varr = sep.GetOutput().GetPointData().GetArray("SelectedPoints") 1723 mask = vtk2numpy(varr).astype(bool) 1724 ids = np.array(range(len(ptsa)), dtype=int)[mask] 1725 1726 if isinstance(pts, Points): 1727 varr.SetName("IsInside") 1728 pts.dataset.GetPointData().AddArray(varr) 1729 1730 if return_ids: 1731 return ids 1732 1733 pcl = Points(ptsa[ids]) 1734 pcl.name = "InsidePoints" 1735 1736 pcl.pipeline = OperationNode( 1737 "inside_points", 1738 parents=[self, ptsa], 1739 comment=f"#pts {pcl.dataset.GetNumberOfPoints()}", 1740 ) 1741 return pcl 1742 1743 def boundaries( 1744 self, 1745 boundary_edges=True, 1746 manifold_edges=False, 1747 non_manifold_edges=False, 1748 feature_angle=None, 1749 return_point_ids=False, 1750 return_cell_ids=False, 1751 cell_edge=False, 1752 ) -> Union[Self, np.ndarray]: 1753 """ 1754 Return the boundary lines of an input mesh. 1755 Check also `vedo.core.CommonAlgorithms.mark_boundaries()` method. 1756 1757 Arguments: 1758 boundary_edges : (bool) 1759 Turn on/off the extraction of boundary edges. 1760 manifold_edges : (bool) 1761 Turn on/off the extraction of manifold edges. 1762 non_manifold_edges : (bool) 1763 Turn on/off the extraction of non-manifold edges. 1764 feature_angle : (bool) 1765 Specify the min angle btw 2 faces for extracting edges. 1766 return_point_ids : (bool) 1767 return a numpy array of point indices 1768 return_cell_ids : (bool) 1769 return a numpy array of cell indices 1770 cell_edge : (bool) 1771 set to `True` if a cell need to share an edge with 1772 the boundary line, or `False` if a single vertex is enough 1773 1774 Examples: 1775 - [boundaries.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/boundaries.py) 1776 1777  1778 """ 1779 fe = vtki.new("FeatureEdges") 1780 fe.SetBoundaryEdges(boundary_edges) 1781 fe.SetNonManifoldEdges(non_manifold_edges) 1782 fe.SetManifoldEdges(manifold_edges) 1783 try: 1784 fe.SetPassLines(True) # vtk9.2 1785 except AttributeError: 1786 pass 1787 fe.ColoringOff() 1788 fe.SetFeatureEdges(False) 1789 if feature_angle is not None: 1790 fe.SetFeatureEdges(True) 1791 fe.SetFeatureAngle(feature_angle) 1792 1793 if return_point_ids or return_cell_ids: 1794 idf = vtki.new("IdFilter") 1795 idf.SetInputData(self.dataset) 1796 idf.SetPointIdsArrayName("BoundaryIds") 1797 idf.SetPointIds(True) 1798 idf.Update() 1799 1800 fe.SetInputData(idf.GetOutput()) 1801 fe.Update() 1802 1803 vid = fe.GetOutput().GetPointData().GetArray("BoundaryIds") 1804 npid = vtk2numpy(vid).astype(int) 1805 1806 if return_point_ids: 1807 return npid 1808 1809 if return_cell_ids: 1810 n = 1 if cell_edge else 0 1811 inface = [] 1812 for i, face in enumerate(self.cells): 1813 # isin = np.any([vtx in npid for vtx in face]) 1814 isin = 0 1815 for vtx in face: 1816 isin += int(vtx in npid) 1817 if isin > n: 1818 break 1819 if isin > n: 1820 inface.append(i) 1821 return np.array(inface).astype(int) 1822 1823 return self 1824 1825 else: 1826 1827 fe.SetInputData(self.dataset) 1828 fe.Update() 1829 msh = Mesh(fe.GetOutput(), c="p").lw(5).lighting("off") 1830 msh.name = "MeshBoundaries" 1831 1832 msh.pipeline = OperationNode( 1833 "boundaries", 1834 parents=[self], 1835 shape="octagon", 1836 comment=f"#pts {msh.dataset.GetNumberOfPoints()}", 1837 ) 1838 return msh 1839 1840 def imprint(self, loopline, tol=0.01) -> Self: 1841 """ 1842 Imprint the contact surface of one object onto another surface. 1843 1844 Arguments: 1845 loopline : (vedo.Line) 1846 a Line object to be imprinted onto the mesh. 1847 tol : (float) 1848 projection tolerance which controls how close the imprint 1849 surface must be to the target. 1850 1851 Example: 1852 ```python 1853 from vedo import * 1854 grid = Grid()#.triangulate() 1855 circle = Circle(r=0.3, res=24).pos(0.11,0.12) 1856 line = Line(circle, closed=True, lw=4, c='r4') 1857 grid.imprint(line) 1858 show(grid, line, axes=1).close() 1859 ``` 1860  1861 """ 1862 loop = vtki.new("ContourLoopExtraction") 1863 loop.SetInputData(loopline.dataset) 1864 loop.Update() 1865 1866 clean_loop = vtki.new("CleanPolyData") 1867 clean_loop.SetInputData(loop.GetOutput()) 1868 clean_loop.Update() 1869 1870 imp = vtki.new("ImprintFilter") 1871 imp.SetTargetData(self.dataset) 1872 imp.SetImprintData(clean_loop.GetOutput()) 1873 imp.SetTolerance(tol) 1874 imp.BoundaryEdgeInsertionOn() 1875 imp.TriangulateOutputOn() 1876 imp.Update() 1877 1878 self._update(imp.GetOutput()) 1879 1880 self.pipeline = OperationNode( 1881 "imprint", 1882 parents=[self], 1883 comment=f"#pts {self.dataset.GetNumberOfPoints()}", 1884 ) 1885 return self 1886 1887 def connected_vertices(self, index: int) -> List[int]: 1888 """Find all vertices connected to an input vertex specified by its index. 1889 1890 Examples: 1891 - [connected_vtx.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/connected_vtx.py) 1892 1893  1894 """ 1895 poly = self.dataset 1896 1897 cell_idlist = vtki.vtkIdList() 1898 poly.GetPointCells(index, cell_idlist) 1899 1900 idxs = [] 1901 for i in range(cell_idlist.GetNumberOfIds()): 1902 point_idlist = vtki.vtkIdList() 1903 poly.GetCellPoints(cell_idlist.GetId(i), point_idlist) 1904 for j in range(point_idlist.GetNumberOfIds()): 1905 idj = point_idlist.GetId(j) 1906 if idj == index: 1907 continue 1908 if idj in idxs: 1909 continue 1910 idxs.append(idj) 1911 1912 return idxs 1913 1914 def extract_cells(self, ids: List[int]) -> Self: 1915 """ 1916 Extract a subset of cells from a mesh and return it as a new mesh. 1917 """ 1918 selectCells = vtki.new("SelectionNode") 1919 selectCells.SetFieldType(vtki.get_class("SelectionNode").CELL) 1920 selectCells.SetContentType(vtki.get_class("SelectionNode").INDICES) 1921 idarr = vtki.vtkIdTypeArray() 1922 idarr.SetNumberOfComponents(1) 1923 idarr.SetNumberOfValues(len(ids)) 1924 for i, v in enumerate(ids): 1925 idarr.SetValue(i, v) 1926 selectCells.SetSelectionList(idarr) 1927 1928 selection = vtki.new("Selection") 1929 selection.AddNode(selectCells) 1930 1931 extractSelection = vtki.new("ExtractSelection") 1932 extractSelection.SetInputData(0, self.dataset) 1933 extractSelection.SetInputData(1, selection) 1934 extractSelection.Update() 1935 1936 gf = vtki.new("GeometryFilter") 1937 gf.SetInputData(extractSelection.GetOutput()) 1938 gf.Update() 1939 msh = Mesh(gf.GetOutput()) 1940 msh.copy_properties_from(self) 1941 return msh 1942 1943 def connected_cells(self, index: int, return_ids=False) -> Union[Self, List[int]]: 1944 """Find all cellls connected to an input vertex specified by its index.""" 1945 1946 # Find all cells connected to point index 1947 dpoly = self.dataset 1948 idlist = vtki.vtkIdList() 1949 dpoly.GetPointCells(index, idlist) 1950 1951 ids = vtki.vtkIdTypeArray() 1952 ids.SetNumberOfComponents(1) 1953 rids = [] 1954 for k in range(idlist.GetNumberOfIds()): 1955 cid = idlist.GetId(k) 1956 ids.InsertNextValue(cid) 1957 rids.append(int(cid)) 1958 if return_ids: 1959 return rids 1960 1961 selection_node = vtki.new("SelectionNode") 1962 selection_node.SetFieldType(vtki.get_class("SelectionNode").CELL) 1963 selection_node.SetContentType(vtki.get_class("SelectionNode").INDICES) 1964 selection_node.SetSelectionList(ids) 1965 selection = vtki.new("Selection") 1966 selection.AddNode(selection_node) 1967 extractSelection = vtki.new("ExtractSelection") 1968 extractSelection.SetInputData(0, dpoly) 1969 extractSelection.SetInputData(1, selection) 1970 extractSelection.Update() 1971 gf = vtki.new("GeometryFilter") 1972 gf.SetInputData(extractSelection.GetOutput()) 1973 gf.Update() 1974 return Mesh(gf.GetOutput()).lw(1) 1975 1976 def silhouette(self, direction=None, border_edges=True, feature_angle=False) -> Self: 1977 """ 1978 Return a new line `Mesh` which corresponds to the outer `silhouette` 1979 of the input as seen along a specified `direction`, this can also be 1980 a `vtkCamera` object. 1981 1982 Arguments: 1983 direction : (list) 1984 viewpoint direction vector. 1985 If `None` this is guessed by looking at the minimum 1986 of the sides of the bounding box. 1987 border_edges : (bool) 1988 enable or disable generation of border edges 1989 feature_angle : (float) 1990 minimal angle for sharp edges detection. 1991 If set to `False` the functionality is disabled. 1992 1993 Examples: 1994 - [silhouette1.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/silhouette1.py) 1995 1996  1997 """ 1998 sil = vtki.new("PolyDataSilhouette") 1999 sil.SetInputData(self.dataset) 2000 sil.SetBorderEdges(border_edges) 2001 if feature_angle is False: 2002 sil.SetEnableFeatureAngle(0) 2003 else: 2004 sil.SetEnableFeatureAngle(1) 2005 sil.SetFeatureAngle(feature_angle) 2006 2007 if direction is None and vedo.plotter_instance and vedo.plotter_instance.camera: 2008 sil.SetCamera(vedo.plotter_instance.camera) 2009 m = Mesh() 2010 m.mapper.SetInputConnection(sil.GetOutputPort()) 2011 2012 elif isinstance(direction, vtki.vtkCamera): 2013 sil.SetCamera(direction) 2014 m = Mesh() 2015 m.mapper.SetInputConnection(sil.GetOutputPort()) 2016 2017 elif direction == "2d": 2018 sil.SetVector(3.4, 4.5, 5.6) # random 2019 sil.SetDirectionToSpecifiedVector() 2020 sil.Update() 2021 m = Mesh(sil.GetOutput()) 2022 2023 elif is_sequence(direction): 2024 sil.SetVector(direction) 2025 sil.SetDirectionToSpecifiedVector() 2026 sil.Update() 2027 m = Mesh(sil.GetOutput()) 2028 else: 2029 vedo.logger.error(f"in silhouette() unknown direction type {type(direction)}") 2030 vedo.logger.error("first render the scene with show() or specify camera/direction") 2031 return self 2032 2033 m.lw(2).c((0, 0, 0)).lighting("off") 2034 m.mapper.SetResolveCoincidentTopologyToPolygonOffset() 2035 m.pipeline = OperationNode("silhouette", parents=[self]) 2036 m.name = "Silhouette" 2037 return m 2038 2039 def isobands(self, n=10, vmin=None, vmax=None) -> Self: 2040 """ 2041 Return a new `Mesh` representing the isobands of the active scalars. 2042 This is a new mesh where the scalar is now associated to cell faces and 2043 used to colorize the mesh. 2044 2045 Arguments: 2046 n : (int) 2047 number of isobands in the range 2048 vmin : (float) 2049 minimum of the range 2050 vmax : (float) 2051 maximum of the range 2052 2053 Examples: 2054 - [isolines.py](https://github.com/marcomusy/vedo/tree/master/examples/pyplot/isolines.py) 2055 """ 2056 r0, r1 = self.dataset.GetScalarRange() 2057 if vmin is None: 2058 vmin = r0 2059 if vmax is None: 2060 vmax = r1 2061 2062 # -------------------------------- 2063 bands = [] 2064 dx = (vmax - vmin) / float(n) 2065 b = [vmin, vmin + dx / 2.0, vmin + dx] 2066 i = 0 2067 while i < n: 2068 bands.append(b) 2069 b = [b[0] + dx, b[1] + dx, b[2] + dx] 2070 i += 1 2071 2072 # annotate, use the midpoint of the band as the label 2073 lut = self.mapper.GetLookupTable() 2074 labels = [] 2075 for b in bands: 2076 labels.append("{:4.2f}".format(b[1])) 2077 values = vtki.vtkVariantArray() 2078 for la in labels: 2079 values.InsertNextValue(vtki.vtkVariant(la)) 2080 for i in range(values.GetNumberOfTuples()): 2081 lut.SetAnnotation(i, values.GetValue(i).ToString()) 2082 2083 bcf = vtki.new("BandedPolyDataContourFilter") 2084 bcf.SetInputData(self.dataset) 2085 # Use either the minimum or maximum value for each band. 2086 for i, band in enumerate(bands): 2087 bcf.SetValue(i, band[2]) 2088 # We will use an indexed lookup table. 2089 bcf.SetScalarModeToIndex() 2090 bcf.GenerateContourEdgesOff() 2091 bcf.Update() 2092 bcf.GetOutput().GetCellData().GetScalars().SetName("IsoBands") 2093 2094 m1 = Mesh(bcf.GetOutput()).compute_normals(cells=True) 2095 m1.mapper.SetLookupTable(lut) 2096 m1.mapper.SetScalarRange(lut.GetRange()) 2097 m1.pipeline = OperationNode("isobands", parents=[self]) 2098 m1.name = "IsoBands" 2099 return m1 2100 2101 def isolines(self, n=10, vmin=None, vmax=None) -> Self: 2102 """ 2103 Return a new `Mesh` representing the isolines of the active scalars. 2104 2105 Arguments: 2106 n : (int) 2107 number of isolines in the range 2108 vmin : (float) 2109 minimum of the range 2110 vmax : (float) 2111 maximum of the range 2112 2113 Examples: 2114 - [isolines.py](https://github.com/marcomusy/vedo/tree/master/examples/pyplot/isolines.py) 2115 2116  2117 """ 2118 bcf = vtki.new("ContourFilter") 2119 bcf.SetInputData(self.dataset) 2120 r0, r1 = self.dataset.GetScalarRange() 2121 if vmin is None: 2122 vmin = r0 2123 if vmax is None: 2124 vmax = r1 2125 bcf.GenerateValues(n, vmin, vmax) 2126 bcf.Update() 2127 sf = vtki.new("Stripper") 2128 sf.SetJoinContiguousSegments(True) 2129 sf.SetInputData(bcf.GetOutput()) 2130 sf.Update() 2131 cl = vtki.new("CleanPolyData") 2132 cl.SetInputData(sf.GetOutput()) 2133 cl.Update() 2134 msh = Mesh(cl.GetOutput(), c="k").lighting("off") 2135 msh.mapper.SetResolveCoincidentTopologyToPolygonOffset() 2136 msh.pipeline = OperationNode("isolines", parents=[self]) 2137 msh.name = "IsoLines" 2138 return msh 2139 2140 def extrude(self, zshift=1.0, direction=(), rotation=0.0, dr=0.0, cap=True, res=1) -> Self: 2141 """ 2142 Sweep a polygonal data creating a "skirt" from free edges and lines, and lines from vertices. 2143 The input dataset is swept around the z-axis to create new polygonal primitives. 2144 For example, sweeping a line results in a cylindrical shell, and sweeping a circle creates a torus. 2145 2146 You can control whether the sweep of a 2D object (i.e., polygon or triangle strip) 2147 is capped with the generating geometry. 2148 Also, you can control the angle of rotation, and whether translation along the z-axis 2149 is performed along with the rotation. (Translation is useful for creating "springs"). 2150 You also can adjust the radius of the generating geometry using the "dR" keyword. 2151 2152 The skirt is generated by locating certain topological features. 2153 Free edges (edges of polygons or triangle strips only used by one polygon or triangle strips) 2154 generate surfaces. This is true also of lines or polylines. Vertices generate lines. 2155 2156 This filter can be used to model axisymmetric objects like cylinders, bottles, and wine glasses; 2157 or translational/rotational symmetric objects like springs or corkscrews. 2158 2159 Arguments: 2160 zshift : (float) 2161 shift along z axis. 2162 direction : (list) 2163 extrusion direction in the xy plane. 2164 note that zshift is forced to be the 3rd component of direction, 2165 which is therefore ignored. 2166 rotation : (float) 2167 set the angle of rotation. 2168 dr : (float) 2169 set the radius variation in absolute units. 2170 cap : (bool) 2171 enable or disable capping. 2172 res : (int) 2173 set the resolution of the generating geometry. 2174 2175 Warning: 2176 Some polygonal objects have no free edges (e.g., sphere). When swept, this will result 2177 in two separate surfaces if capping is on, or no surface if capping is off. 2178 2179 Examples: 2180 - [extrude.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/extrude.py) 2181 2182  2183 """ 2184 rf = vtki.new("RotationalExtrusionFilter") 2185 # rf = vtki.new("LinearExtrusionFilter") 2186 rf.SetInputData(self.dataset) # must not be transformed 2187 rf.SetResolution(res) 2188 rf.SetCapping(cap) 2189 rf.SetAngle(rotation) 2190 rf.SetTranslation(zshift) 2191 rf.SetDeltaRadius(dr) 2192 rf.Update() 2193 2194 # convert triangle strips to polygonal data 2195 tris = vtki.new("TriangleFilter") 2196 tris.SetInputData(rf.GetOutput()) 2197 tris.Update() 2198 2199 m = Mesh(tris.GetOutput()) 2200 2201 if len(direction) > 1: 2202 p = self.pos() 2203 LT = vedo.LinearTransform() 2204 LT.translate(-p) 2205 LT.concatenate([ 2206 [1, 0, direction[0]], 2207 [0, 1, direction[1]], 2208 [0, 0, 1] 2209 ]) 2210 LT.translate(p) 2211 m.apply_transform(LT) 2212 2213 m.copy_properties_from(self).flat().lighting("default") 2214 m.pipeline = OperationNode( 2215 "extrude", parents=[self], 2216 comment=f"#pts {m.dataset.GetNumberOfPoints()}" 2217 ) 2218 m.name = "ExtrudedMesh" 2219 return m 2220 2221 def extrude_and_trim_with( 2222 self, 2223 surface: "Mesh", 2224 direction=(), 2225 strategy="all", 2226 cap=True, 2227 cap_strategy="max", 2228 ) -> Self: 2229 """ 2230 Extrude a Mesh and trim it with an input surface mesh. 2231 2232 Arguments: 2233 surface : (Mesh) 2234 the surface mesh to trim with. 2235 direction : (list) 2236 extrusion direction in the xy plane. 2237 strategy : (str) 2238 either "boundary_edges" or "all_edges". 2239 cap : (bool) 2240 enable or disable capping. 2241 cap_strategy : (str) 2242 either "intersection", "minimum_distance", "maximum_distance", "average_distance". 2243 2244 The input Mesh is swept along a specified direction forming a "skirt" 2245 from the boundary edges 2D primitives (i.e., edges used by only one polygon); 2246 and/or from vertices and lines. 2247 The extent of the sweeping is limited by a second input: defined where 2248 the sweep intersects a user-specified surface. 2249 2250 Capping of the extrusion can be enabled. 2251 In this case the input, generating primitive is copied inplace as well 2252 as to the end of the extrusion skirt. 2253 (See warnings below on what happens if the intersecting sweep does not 2254 intersect, or partially intersects the trim surface.) 2255 2256 Note that this method operates in two fundamentally different modes 2257 based on the extrusion strategy. 2258 If the strategy is "boundary_edges", then only the boundary edges of the input's 2259 2D primitives are extruded (verts and lines are extruded to generate lines and quads). 2260 However, if the extrusions strategy is "all_edges", then every edge of the 2D primitives 2261 is used to sweep out a quadrilateral polygon (again verts and lines are swept to produce lines and quads). 2262 2263 Warning: 2264 The extrusion direction is assumed to define an infinite line. 2265 The intersection with the trim surface is along a ray from the - to + direction, 2266 however only the first intersection is taken. 2267 Some polygonal objects have no free edges (e.g., sphere). When swept, this will result in two separate 2268 surfaces if capping is on and "boundary_edges" enabled, 2269 or no surface if capping is off and "boundary_edges" is enabled. 2270 If all the extrusion lines emanating from an extruding primitive do not intersect the trim surface, 2271 then no output for that primitive will be generated. In extreme cases, it is possible that no output 2272 whatsoever will be generated. 2273 2274 Example: 2275 ```python 2276 from vedo import * 2277 sphere = Sphere([-1,0,4]).rotate_x(25).wireframe().color('red5') 2278 circle = Circle([0,0,0], r=2, res=100).color('b6') 2279 extruded_circle = circle.extrude_and_trim_with( 2280 sphere, 2281 direction=[0,-0.2,1], 2282 strategy="bound", 2283 cap=True, 2284 cap_strategy="intersection", 2285 ) 2286 circle.lw(3).color("tomato").shift(dz=-0.1) 2287 show(circle, sphere, extruded_circle, axes=1).close() 2288 ``` 2289 """ 2290 trimmer = vtki.new("TrimmedExtrusionFilter") 2291 trimmer.SetInputData(self.dataset) 2292 trimmer.SetCapping(cap) 2293 trimmer.SetExtrusionDirection(direction) 2294 trimmer.SetTrimSurfaceData(surface.dataset) 2295 if "bound" in strategy: 2296 trimmer.SetExtrusionStrategyToBoundaryEdges() 2297 elif "all" in strategy: 2298 trimmer.SetExtrusionStrategyToAllEdges() 2299 else: 2300 vedo.logger.warning(f"extrude_and_trim(): unknown strategy {strategy}") 2301 # print (trimmer.GetExtrusionStrategy()) 2302 2303 if "intersect" in cap_strategy: 2304 trimmer.SetCappingStrategyToIntersection() 2305 elif "min" in cap_strategy: 2306 trimmer.SetCappingStrategyToMinimumDistance() 2307 elif "max" in cap_strategy: 2308 trimmer.SetCappingStrategyToMaximumDistance() 2309 elif "ave" in cap_strategy: 2310 trimmer.SetCappingStrategyToAverageDistance() 2311 else: 2312 vedo.logger.warning(f"extrude_and_trim(): unknown cap_strategy {cap_strategy}") 2313 # print (trimmer.GetCappingStrategy()) 2314 2315 trimmer.Update() 2316 2317 m = Mesh(trimmer.GetOutput()) 2318 m.copy_properties_from(self).flat().lighting("default") 2319 m.pipeline = OperationNode( 2320 "extrude_and_trim", parents=[self, surface], 2321 comment=f"#pts {m.dataset.GetNumberOfPoints()}" 2322 ) 2323 m.name = "ExtrudedAndTrimmedMesh" 2324 return m 2325 2326 def split( 2327 self, maxdepth=1000, flag=False, must_share_edge=False, sort_by_area=True 2328 ) -> List[Self]: 2329 """ 2330 Split a mesh by connectivity and order the pieces by increasing area. 2331 2332 Arguments: 2333 maxdepth : (int) 2334 only consider this maximum number of mesh parts. 2335 flag : (bool) 2336 if set to True return the same single object, 2337 but add a "RegionId" array to flag the mesh subparts 2338 must_share_edge : (bool) 2339 if True, mesh regions that only share single points will be split. 2340 sort_by_area : (bool) 2341 if True, sort the mesh parts by decreasing area. 2342 2343 Examples: 2344 - [splitmesh.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/splitmesh.py) 2345 2346  2347 """ 2348 pd = self.dataset 2349 if must_share_edge: 2350 if pd.GetNumberOfPolys() == 0: 2351 vedo.logger.warning("in split(): no polygons found. Skip.") 2352 return [self] 2353 cf = vtki.new("PolyDataEdgeConnectivityFilter") 2354 cf.BarrierEdgesOff() 2355 else: 2356 cf = vtki.new("PolyDataConnectivityFilter") 2357 2358 cf.SetInputData(pd) 2359 cf.SetExtractionModeToAllRegions() 2360 cf.SetColorRegions(True) 2361 cf.Update() 2362 out = cf.GetOutput() 2363 2364 if not out.GetNumberOfPoints(): 2365 return [self] 2366 2367 if flag: 2368 self.pipeline = OperationNode("split mesh", parents=[self]) 2369 self._update(out) 2370 return [self] 2371 2372 msh = Mesh(out) 2373 if must_share_edge: 2374 arr = msh.celldata["RegionId"] 2375 on = "cells" 2376 else: 2377 arr = msh.pointdata["RegionId"] 2378 on = "points" 2379 2380 alist = [] 2381 for t in range(max(arr) + 1): 2382 if t == maxdepth: 2383 break 2384 suba = msh.clone().threshold("RegionId", t, t, on=on) 2385 if sort_by_area: 2386 area = suba.area() 2387 else: 2388 area = 0 # dummy 2389 suba.name = "MeshRegion" + str(t) 2390 alist.append([suba, area]) 2391 2392 if sort_by_area: 2393 alist.sort(key=lambda x: x[1]) 2394 alist.reverse() 2395 2396 blist = [] 2397 for i, l in enumerate(alist): 2398 l[0].color(i + 1).phong() 2399 l[0].mapper.ScalarVisibilityOff() 2400 blist.append(l[0]) 2401 if i < 10: 2402 l[0].pipeline = OperationNode( 2403 f"split mesh {i}", 2404 parents=[self], 2405 comment=f"#pts {l[0].dataset.GetNumberOfPoints()}", 2406 ) 2407 return blist 2408 2409 def extract_largest_region(self) -> Self: 2410 """ 2411 Extract the largest connected part of a mesh and discard all the smaller pieces. 2412 2413 Examples: 2414 - [largestregion.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/largestregion.py) 2415 """ 2416 conn = vtki.new("PolyDataConnectivityFilter") 2417 conn.SetExtractionModeToLargestRegion() 2418 conn.ScalarConnectivityOff() 2419 conn.SetInputData(self.dataset) 2420 conn.Update() 2421 2422 m = Mesh(conn.GetOutput()) 2423 m.copy_properties_from(self) 2424 m.pipeline = OperationNode( 2425 "extract_largest_region", 2426 parents=[self], 2427 comment=f"#pts {m.dataset.GetNumberOfPoints()}", 2428 ) 2429 m.name = "MeshLargestRegion" 2430 return m 2431 2432 def boolean(self, operation: str, mesh2, method=0, tol=None) -> Self: 2433 """Volumetric union, intersection and subtraction of surfaces. 2434 2435 Use `operation` for the allowed operations `['plus', 'intersect', 'minus']`. 2436 2437 Two possible algorithms are available. 2438 Setting `method` to 0 (the default) uses the boolean operation algorithm 2439 written by Cory Quammen, Chris Weigle, and Russ Taylor (https://doi.org/10.54294/216g01); 2440 setting `method` to 1 will use the "loop" boolean algorithm 2441 written by Adam Updegrove (https://doi.org/10.1016/j.advengsoft.2016.01.015). 2442 2443 Use `tol` to specify the absolute tolerance used to determine 2444 when the distance between two points is considered to be zero (defaults to 1e-6). 2445 2446 Example: 2447 - [boolean.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/boolean.py) 2448 2449  2450 """ 2451 if method == 0: 2452 bf = vtki.new("BooleanOperationPolyDataFilter") 2453 elif method == 1: 2454 bf = vtki.new("LoopBooleanPolyDataFilter") 2455 else: 2456 raise ValueError(f"Unknown method={method}") 2457 2458 poly1 = self.compute_normals().dataset 2459 poly2 = mesh2.compute_normals().dataset 2460 2461 if operation.lower() in ("plus", "+"): 2462 bf.SetOperationToUnion() 2463 elif operation.lower() == "intersect": 2464 bf.SetOperationToIntersection() 2465 elif operation.lower() in ("minus", "-"): 2466 bf.SetOperationToDifference() 2467 2468 if tol: 2469 bf.SetTolerance(tol) 2470 2471 bf.SetInputData(0, poly1) 2472 bf.SetInputData(1, poly2) 2473 bf.Update() 2474 2475 msh = Mesh(bf.GetOutput(), c=None) 2476 msh.flat() 2477 2478 msh.pipeline = OperationNode( 2479 "boolean " + operation, 2480 parents=[self, mesh2], 2481 shape="cylinder", 2482 comment=f"#pts {msh.dataset.GetNumberOfPoints()}", 2483 ) 2484 msh.name = self.name + operation + mesh2.name 2485 return msh 2486 2487 def intersect_with(self, mesh2, tol=1e-06) -> Self: 2488 """ 2489 Intersect this Mesh with the input surface to return a set of lines. 2490 2491 Examples: 2492 - [surf_intersect.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/surf_intersect.py) 2493 2494  2495 """ 2496 bf = vtki.new("IntersectionPolyDataFilter") 2497 bf.SetGlobalWarningDisplay(0) 2498 bf.SetTolerance(tol) 2499 bf.SetInputData(0, self.dataset) 2500 bf.SetInputData(1, mesh2.dataset) 2501 bf.Update() 2502 msh = Mesh(bf.GetOutput(), c="k", alpha=1).lighting("off") 2503 msh.properties.SetLineWidth(3) 2504 msh.pipeline = OperationNode( 2505 "intersect_with", parents=[self, mesh2], comment=f"#pts {msh.npoints}" 2506 ) 2507 msh.name = "SurfaceIntersection" 2508 return msh 2509 2510 def intersect_with_line(self, p0, p1=None, return_ids=False, tol=0) -> Union[np.ndarray, Tuple[np.ndarray, np.ndarray]]: 2511 """ 2512 Return the list of points intersecting the mesh 2513 along the segment defined by two points `p0` and `p1`. 2514 2515 Use `return_ids` to return the cell ids along with point coords 2516 2517 Example: 2518 ```python 2519 from vedo import * 2520 s = Spring() 2521 pts = s.intersect_with_line([0,0,0], [1,0.1,0]) 2522 ln = Line([0,0,0], [1,0.1,0], c='blue') 2523 ps = Points(pts, r=10, c='r') 2524 show(s, ln, ps, bg='white').close() 2525 ``` 2526  2527 """ 2528 if isinstance(p0, Points): 2529 p0, p1 = p0.vertices 2530 2531 if not self.line_locator: 2532 self.line_locator = vtki.new("OBBTree") 2533 self.line_locator.SetDataSet(self.dataset) 2534 if not tol: 2535 tol = mag(np.asarray(p1) - np.asarray(p0)) / 10000 2536 self.line_locator.SetTolerance(tol) 2537 self.line_locator.BuildLocator() 2538 2539 vpts = vtki.vtkPoints() 2540 idlist = vtki.vtkIdList() 2541 self.line_locator.IntersectWithLine(p0, p1, vpts, idlist) 2542 pts = [] 2543 for i in range(vpts.GetNumberOfPoints()): 2544 intersection: MutableSequence[float] = [0, 0, 0] 2545 vpts.GetPoint(i, intersection) 2546 pts.append(intersection) 2547 pts2 = np.array(pts) 2548 2549 if return_ids: 2550 pts_ids = [] 2551 for i in range(idlist.GetNumberOfIds()): 2552 cid = idlist.GetId(i) 2553 pts_ids.append(cid) 2554 return (pts2, np.array(pts_ids).astype(np.uint32)) 2555 2556 return pts2 2557 2558 def intersect_with_plane(self, origin=(0, 0, 0), normal=(1, 0, 0)) -> Self: 2559 """ 2560 Intersect this Mesh with a plane to return a set of lines. 2561 2562 Example: 2563 ```python 2564 from vedo import * 2565 sph = Sphere() 2566 mi = sph.clone().intersect_with_plane().join() 2567 print(mi.lines) 2568 show(sph, mi, axes=1).close() 2569 ``` 2570  2571 """ 2572 plane = vtki.new("Plane") 2573 plane.SetOrigin(origin) 2574 plane.SetNormal(normal) 2575 2576 cutter = vtki.new("PolyDataPlaneCutter") 2577 cutter.SetInputData(self.dataset) 2578 cutter.SetPlane(plane) 2579 cutter.InterpolateAttributesOn() 2580 cutter.ComputeNormalsOff() 2581 cutter.Update() 2582 2583 msh = Mesh(cutter.GetOutput()) 2584 msh.c('k').lw(3).lighting("off") 2585 msh.pipeline = OperationNode( 2586 "intersect_with_plan", 2587 parents=[self], 2588 comment=f"#pts {msh.dataset.GetNumberOfPoints()}", 2589 ) 2590 msh.name = "PlaneIntersection" 2591 return msh 2592 2593 def cut_closed_surface(self, origins, normals, invert=False, return_assembly=False) -> Union[Self, "vedo.Assembly"]: 2594 """ 2595 Cut/clip a closed surface mesh with a collection of planes. 2596 This will produce a new closed surface by creating new polygonal 2597 faces where the input surface hits the planes. 2598 2599 The orientation of the polygons that form the surface is important. 2600 Polygons have a front face and a back face, and it's the back face that defines 2601 the interior or "solid" region of the closed surface. 2602 When a plane cuts through a "solid" region, a new cut face is generated, 2603 but not when a clipping plane cuts through a hole or "empty" region. 2604 This distinction is crucial when dealing with complex surfaces. 2605 Note that if a simple surface has its back faces pointing outwards, 2606 then that surface defines a hole in a potentially infinite solid. 2607 2608 Non-manifold surfaces should not be used with this method. 2609 2610 Arguments: 2611 origins : (list) 2612 list of plane origins 2613 normals : (list) 2614 list of plane normals 2615 invert : (bool) 2616 invert the clipping. 2617 return_assembly : (bool) 2618 return the cap and the clipped surfaces as a `vedo.Assembly`. 2619 2620 Example: 2621 ```python 2622 from vedo import * 2623 s = Sphere(res=50).linewidth(1) 2624 origins = [[-0.7, 0, 0], [0, -0.6, 0]] 2625 normals = [[-1, 0, 0], [0, -1, 0]] 2626 s.cut_closed_surface(origins, normals) 2627 show(s, axes=1).close() 2628 ``` 2629 """ 2630 planes = vtki.new("PlaneCollection") 2631 for p, s in zip(origins, normals): 2632 plane = vtki.vtkPlane() 2633 plane.SetOrigin(vedo.utils.make3d(p)) 2634 plane.SetNormal(vedo.utils.make3d(s)) 2635 planes.AddItem(plane) 2636 clipper = vtki.new("ClipClosedSurface") 2637 clipper.SetInputData(self.dataset) 2638 clipper.SetClippingPlanes(planes) 2639 clipper.PassPointDataOn() 2640 clipper.GenerateFacesOn() 2641 clipper.SetScalarModeToLabels() 2642 clipper.TriangulationErrorDisplayOn() 2643 clipper.SetInsideOut(not invert) 2644 2645 if return_assembly: 2646 clipper.GenerateClipFaceOutputOn() 2647 clipper.Update() 2648 parts = [] 2649 for i in range(clipper.GetNumberOfOutputPorts()): 2650 msh = Mesh(clipper.GetOutput(i)) 2651 msh.copy_properties_from(self) 2652 msh.name = "CutClosedSurface" 2653 msh.pipeline = OperationNode( 2654 "cut_closed_surface", 2655 parents=[self], 2656 comment=f"#pts {msh.dataset.GetNumberOfPoints()}", 2657 ) 2658 parts.append(msh) 2659 asse = vedo.Assembly(parts) 2660 asse.name = "CutClosedSurface" 2661 return asse 2662 2663 else: 2664 clipper.GenerateClipFaceOutputOff() 2665 clipper.Update() 2666 self._update(clipper.GetOutput()) 2667 self.flat() 2668 self.name = "CutClosedSurface" 2669 self.pipeline = OperationNode( 2670 "cut_closed_surface", 2671 parents=[self], 2672 comment=f"#pts {self.dataset.GetNumberOfPoints()}", 2673 ) 2674 return self 2675 2676 def collide_with(self, mesh2, tol=0, return_bool=False) -> Union[Self, bool]: 2677 """ 2678 Collide this Mesh with the input surface. 2679 Information is stored in `ContactCells1` and `ContactCells2`. 2680 """ 2681 ipdf = vtki.new("CollisionDetectionFilter") 2682 # ipdf.SetGlobalWarningDisplay(0) 2683 2684 transform0 = vtki.vtkTransform() 2685 transform1 = vtki.vtkTransform() 2686 2687 # ipdf.SetBoxTolerance(tol) 2688 ipdf.SetCellTolerance(tol) 2689 ipdf.SetInputData(0, self.dataset) 2690 ipdf.SetInputData(1, mesh2.dataset) 2691 ipdf.SetTransform(0, transform0) 2692 ipdf.SetTransform(1, transform1) 2693 if return_bool: 2694 ipdf.SetCollisionModeToFirstContact() 2695 else: 2696 ipdf.SetCollisionModeToAllContacts() 2697 ipdf.Update() 2698 2699 if return_bool: 2700 return bool(ipdf.GetNumberOfContacts()) 2701 2702 msh = Mesh(ipdf.GetContactsOutput(), "k", 1).lighting("off") 2703 msh.metadata["ContactCells1"] = vtk2numpy( 2704 ipdf.GetOutput(0).GetFieldData().GetArray("ContactCells") 2705 ) 2706 msh.metadata["ContactCells2"] = vtk2numpy( 2707 ipdf.GetOutput(1).GetFieldData().GetArray("ContactCells") 2708 ) 2709 msh.properties.SetLineWidth(3) 2710 2711 msh.pipeline = OperationNode( 2712 "collide_with", 2713 parents=[self, mesh2], 2714 comment=f"#pts {msh.dataset.GetNumberOfPoints()}", 2715 ) 2716 msh.name = "SurfaceCollision" 2717 return msh 2718 2719 def geodesic(self, start, end) -> Self: 2720 """ 2721 Dijkstra algorithm to compute the geodesic line. 2722 Takes as input a polygonal mesh and performs a single source shortest path calculation. 2723 2724 The output mesh contains the array "VertexIDs" that contains the ordered list of vertices 2725 traversed to get from the start vertex to the end vertex. 2726 2727 Arguments: 2728 start : (int, list) 2729 start vertex index or close point `[x,y,z]` 2730 end : (int, list) 2731 end vertex index or close point `[x,y,z]` 2732 2733 Examples: 2734 - [geodesic_curve.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/geodesic_curve.py) 2735 2736  2737 """ 2738 if is_sequence(start): 2739 cc = self.vertices 2740 pa = Points(cc) 2741 start = pa.closest_point(start, return_point_id=True) 2742 end = pa.closest_point(end, return_point_id=True) 2743 2744 dijkstra = vtki.new("DijkstraGraphGeodesicPath") 2745 dijkstra.SetInputData(self.dataset) 2746 dijkstra.SetStartVertex(end) # inverted in vtk 2747 dijkstra.SetEndVertex(start) 2748 dijkstra.Update() 2749 2750 weights = vtki.vtkDoubleArray() 2751 dijkstra.GetCumulativeWeights(weights) 2752 2753 idlist = dijkstra.GetIdList() 2754 ids = [idlist.GetId(i) for i in range(idlist.GetNumberOfIds())] 2755 2756 length = weights.GetMaxId() + 1 2757 arr = np.zeros(length) 2758 for i in range(length): 2759 arr[i] = weights.GetTuple(i)[0] 2760 2761 poly = dijkstra.GetOutput() 2762 2763 vdata = numpy2vtk(arr) 2764 vdata.SetName("CumulativeWeights") 2765 poly.GetPointData().AddArray(vdata) 2766 2767 vdata2 = numpy2vtk(ids, dtype=np.uint) 2768 vdata2.SetName("VertexIDs") 2769 poly.GetPointData().AddArray(vdata2) 2770 poly.GetPointData().Modified() 2771 2772 dmesh = Mesh(poly).copy_properties_from(self) 2773 dmesh.lw(3).alpha(1).lighting("off") 2774 dmesh.name = "GeodesicLine" 2775 2776 dmesh.pipeline = OperationNode( 2777 "GeodesicLine", 2778 parents=[self], 2779 comment=f"#steps {poly.GetNumberOfPoints()}", 2780 ) 2781 return dmesh 2782 2783 ##################################################################### 2784 ### Stuff returning a Volume object 2785 ##################################################################### 2786 def binarize( 2787 self, 2788 values=(255, 0), 2789 spacing=None, 2790 dims=None, 2791 origin=None, 2792 ) -> "vedo.Volume": 2793 """ 2794 Convert a `Mesh` into a `Volume` where 2795 the interior voxels value is set to `values[0]` (255 by default), while 2796 the exterior voxels value is set to `values[1]` (0 by default). 2797 2798 Arguments: 2799 values : (list) 2800 background and foreground values. 2801 spacing : (list) 2802 voxel spacing in x, y and z. 2803 dims : (list) 2804 dimensions (nr. of voxels) of the output volume. 2805 origin : (list) 2806 position in space of the (0,0,0) voxel. 2807 2808 Examples: 2809 - [mesh2volume.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/mesh2volume.py) 2810 2811  2812 """ 2813 assert len(values) == 2, "values must be a list of 2 values" 2814 fg_value, bg_value = values 2815 2816 bounds = self.bounds() 2817 if spacing is None: # compute spacing 2818 spacing = [0, 0, 0] 2819 diagonal = np.sqrt( 2820 (bounds[1] - bounds[0]) ** 2 2821 + (bounds[3] - bounds[2]) ** 2 2822 + (bounds[5] - bounds[4]) ** 2 2823 ) 2824 spacing[0] = spacing[1] = spacing[2] = diagonal / 250.0 2825 2826 if dims is None: # compute dimensions 2827 dim = [0, 0, 0] 2828 for i in [0, 1, 2]: 2829 dim[i] = int(np.ceil((bounds[i*2+1] - bounds[i*2]) / spacing[i])) 2830 else: 2831 dim = dims 2832 2833 white_img = vtki.vtkImageData() 2834 white_img.SetDimensions(dim) 2835 white_img.SetSpacing(spacing) 2836 white_img.SetExtent(0, dim[0]-1, 0, dim[1]-1, 0, dim[2]-1) 2837 2838 if origin is None: 2839 origin = [0, 0, 0] 2840 origin[0] = bounds[0] + spacing[0] 2841 origin[1] = bounds[2] + spacing[1] 2842 origin[2] = bounds[4] + spacing[2] 2843 white_img.SetOrigin(origin) 2844 2845 # if direction_matrix is not None: 2846 # white_img.SetDirectionMatrix(direction_matrix) 2847 2848 white_img.AllocateScalars(vtki.VTK_UNSIGNED_CHAR, 1) 2849 2850 # fill the image with foreground voxels: 2851 white_img.GetPointData().GetScalars().Fill(fg_value) 2852 2853 # polygonal data --> image stencil: 2854 pol2stenc = vtki.new("PolyDataToImageStencil") 2855 pol2stenc.SetInputData(self.dataset) 2856 pol2stenc.SetOutputOrigin(white_img.GetOrigin()) 2857 pol2stenc.SetOutputSpacing(white_img.GetSpacing()) 2858 pol2stenc.SetOutputWholeExtent(white_img.GetExtent()) 2859 pol2stenc.Update() 2860 2861 # cut the corresponding white image and set the background: 2862 imgstenc = vtki.new("ImageStencil") 2863 imgstenc.SetInputData(white_img) 2864 imgstenc.SetStencilConnection(pol2stenc.GetOutputPort()) 2865 # imgstenc.SetReverseStencil(True) 2866 imgstenc.SetBackgroundValue(bg_value) 2867 imgstenc.Update() 2868 2869 vol = vedo.Volume(imgstenc.GetOutput()) 2870 vol.name = "BinarizedVolume" 2871 vol.pipeline = OperationNode( 2872 "binarize", 2873 parents=[self], 2874 comment=f"dims={tuple(vol.dimensions())}", 2875 c="#e9c46a:#0096c7", 2876 ) 2877 return vol 2878 2879 def signed_distance(self, bounds=None, dims=(20, 20, 20), invert=False, maxradius=None) -> "vedo.Volume": 2880 """ 2881 Compute the `Volume` object whose voxels contains 2882 the signed distance from the mesh. 2883 2884 Arguments: 2885 bounds : (list) 2886 bounds of the output volume 2887 dims : (list) 2888 dimensions (nr. of voxels) of the output volume 2889 invert : (bool) 2890 flip the sign 2891 2892 Examples: 2893 - [volume_from_mesh.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/volume_from_mesh.py) 2894 """ 2895 if maxradius is not None: 2896 vedo.logger.warning( 2897 "in signedDistance(maxradius=...) is ignored. (Only valid for pointclouds)." 2898 ) 2899 if bounds is None: 2900 bounds = self.bounds() 2901 sx = (bounds[1] - bounds[0]) / dims[0] 2902 sy = (bounds[3] - bounds[2]) / dims[1] 2903 sz = (bounds[5] - bounds[4]) / dims[2] 2904 2905 img = vtki.vtkImageData() 2906 img.SetDimensions(dims) 2907 img.SetSpacing(sx, sy, sz) 2908 img.SetOrigin(bounds[0], bounds[2], bounds[4]) 2909 img.AllocateScalars(vtki.VTK_FLOAT, 1) 2910 2911 imp = vtki.new("ImplicitPolyDataDistance") 2912 imp.SetInput(self.dataset) 2913 b2 = bounds[2] 2914 b4 = bounds[4] 2915 d0, d1, d2 = dims 2916 2917 for i in range(d0): 2918 x = i * sx + bounds[0] 2919 for j in range(d1): 2920 y = j * sy + b2 2921 for k in range(d2): 2922 v = imp.EvaluateFunction((x, y, k * sz + b4)) 2923 if invert: 2924 v = -v 2925 img.SetScalarComponentFromFloat(i, j, k, 0, v) 2926 2927 vol = vedo.Volume(img) 2928 vol.name = "SignedVolume" 2929 2930 vol.pipeline = OperationNode( 2931 "signed_distance", 2932 parents=[self], 2933 comment=f"dims={tuple(vol.dimensions())}", 2934 c="#e9c46a:#0096c7", 2935 ) 2936 return vol 2937 2938 def tetralize( 2939 self, 2940 side=0.02, 2941 nmax=300_000, 2942 gap=None, 2943 subsample=False, 2944 uniform=True, 2945 seed=0, 2946 debug=False, 2947 ) -> "vedo.TetMesh": 2948 """ 2949 Tetralize a closed polygonal mesh. Return a `TetMesh`. 2950 2951 Arguments: 2952 side : (float) 2953 desired side of the single tetras as fraction of the bounding box diagonal. 2954 Typical values are in the range (0.01 - 0.03) 2955 nmax : (int) 2956 maximum random numbers to be sampled in the bounding box 2957 gap : (float) 2958 keep this minimum distance from the surface, 2959 if None an automatic choice is made. 2960 subsample : (bool) 2961 subsample input surface, the geometry might be affected 2962 (the number of original faces reduceed), but higher tet quality might be obtained. 2963 uniform : (bool) 2964 generate tets more uniformly packed in the interior of the mesh 2965 seed : (int) 2966 random number generator seed 2967 debug : (bool) 2968 show an intermediate plot with sampled points 2969 2970 Examples: 2971 - [tetralize_surface.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/tetralize_surface.py) 2972 2973  2974 """ 2975 surf = self.clone().clean().compute_normals() 2976 d = surf.diagonal_size() 2977 if gap is None: 2978 gap = side * d * np.sqrt(2 / 3) 2979 n = int(min((1 / side) ** 3, nmax)) 2980 2981 # fill the space w/ points 2982 x0, x1, y0, y1, z0, z1 = surf.bounds() 2983 2984 if uniform: 2985 pts = vedo.utils.pack_spheres([x0, x1, y0, y1, z0, z1], side * d * 1.42) 2986 pts += np.random.randn(len(pts), 3) * side * d * 1.42 / 100 # some small jitter 2987 else: 2988 disp = np.array([x0 + x1, y0 + y1, z0 + z1]) / 2 2989 np.random.seed(seed) 2990 pts = (np.random.rand(n, 3) - 0.5) * np.array([x1 - x0, y1 - y0, z1 - z0]) + disp 2991 2992 normals = surf.celldata["Normals"] 2993 cc = surf.cell_centers 2994 subpts = cc - normals * gap * 1.05 2995 pts = pts.tolist() + subpts.tolist() 2996 2997 if debug: 2998 print(".. tetralize(): subsampling and cleaning") 2999 3000 fillpts = surf.inside_points(pts) 3001 fillpts.subsample(side) 3002 3003 if gap: 3004 fillpts.distance_to(surf) 3005 fillpts.threshold("Distance", above=gap) 3006 3007 if subsample: 3008 surf.subsample(side) 3009 3010 merged_fs = vedo.merge(fillpts, surf) 3011 tmesh = merged_fs.generate_delaunay3d() 3012 tcenters = tmesh.cell_centers 3013 3014 ids = surf.inside_points(tcenters, return_ids=True) 3015 ins = np.zeros(tmesh.ncells) 3016 ins[ids] = 1 3017 3018 if debug: 3019 # vedo.pyplot.histogram(fillpts.pointdata["Distance"], xtitle=f"gap={gap}").show().close() 3020 edges = self.edges 3021 points = self.vertices 3022 elen = mag(points[edges][:, 0, :] - points[edges][:, 1, :]) 3023 histo = vedo.pyplot.histogram(elen, xtitle="edge length", xlim=(0, 3 * side * d)) 3024 print(".. edges min, max", elen.min(), elen.max()) 3025 fillpts.cmap("bone") 3026 vedo.show( 3027 [ 3028 [ 3029 f"This is a debug plot.\n\nGenerated points: {n}\ngap: {gap}", 3030 surf.wireframe().alpha(0.2), 3031 vedo.addons.Axes(surf), 3032 fillpts, 3033 Points(subpts).c("r4").ps(3), 3034 ], 3035 [f"Edges mean length: {np.mean(elen)}\n\nPress q to continue", histo], 3036 ], 3037 N=2, 3038 sharecam=False, 3039 new=True, 3040 ).close() 3041 print(".. thresholding") 3042 3043 tmesh.celldata["inside"] = ins.astype(np.uint8) 3044 tmesh.threshold("inside", above=0.9) 3045 tmesh.celldata.remove("inside") 3046 3047 if debug: 3048 print(f".. tetralize() completed, ntets = {tmesh.ncells}") 3049 3050 tmesh.pipeline = OperationNode( 3051 "tetralize", 3052 parents=[self], 3053 comment=f"#tets = {tmesh.ncells}", 3054 c="#e9c46a:#9e2a2b", 3055 ) 3056 return tmesh
29class Mesh(MeshVisual, Points): 30 """ 31 Build an instance of object `Mesh` derived from `vedo.PointCloud`. 32 """ 33 34 def __init__(self, inputobj=None, c="gold", alpha=1): 35 """ 36 Initialize a ``Mesh`` object. 37 38 Arguments: 39 inputobj : (str, vtkPolyData, vtkActor, vedo.Mesh) 40 If inputobj is `None` an empty mesh is created. 41 If inputobj is a `str` then it is interpreted as the name of a file to load as mesh. 42 If inputobj is an `vtkPolyData` or `vtkActor` or `vedo.Mesh` 43 then a shallow copy of it is created. 44 If inputobj is a `vedo.Mesh` then a shallow copy of it is created. 45 46 Examples: 47 - [buildmesh.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/buildmesh.py) 48 (and many others!) 49 50  51 """ 52 # print("INIT MESH", super()) 53 super().__init__() 54 55 self.name = "Mesh" 56 57 if inputobj is None: 58 # self.dataset = vtki.vtkPolyData() 59 pass 60 61 elif isinstance(inputobj, str): 62 self.dataset = vedo.file_io.load(inputobj).dataset 63 self.filename = inputobj 64 65 elif isinstance(inputobj, vtki.vtkPolyData): 66 # self.dataset.DeepCopy(inputobj) # NO 67 self.dataset = inputobj 68 if self.dataset.GetNumberOfCells() == 0: 69 carr = vtki.vtkCellArray() 70 for i in range(inputobj.GetNumberOfPoints()): 71 carr.InsertNextCell(1) 72 carr.InsertCellPoint(i) 73 self.dataset.SetVerts(carr) 74 75 elif isinstance(inputobj, Mesh): 76 self.dataset = inputobj.dataset 77 78 elif is_sequence(inputobj): 79 ninp = len(inputobj) 80 if ninp == 4: # assume input is [vertices, faces, lines, strips] 81 self.dataset = buildPolyData(inputobj[0], inputobj[1], inputobj[2], inputobj[3]) 82 elif ninp == 3: # assume input is [vertices, faces, lines] 83 self.dataset = buildPolyData(inputobj[0], inputobj[1], inputobj[2]) 84 elif ninp == 2: # assume input is [vertices, faces] 85 self.dataset = buildPolyData(inputobj[0], inputobj[1]) 86 elif ninp == 1: # assume input is [vertices] 87 self.dataset = buildPolyData(inputobj[0]) 88 else: 89 vedo.logger.error("input must be a list of max 4 elements.") 90 raise ValueError() 91 92 elif isinstance(inputobj, vtki.vtkActor): 93 self.dataset.DeepCopy(inputobj.GetMapper().GetInput()) 94 v = inputobj.GetMapper().GetScalarVisibility() 95 self.mapper.SetScalarVisibility(v) 96 pr = vtki.vtkProperty() 97 pr.DeepCopy(inputobj.GetProperty()) 98 self.actor.SetProperty(pr) 99 self.properties = pr 100 101 elif isinstance(inputobj, (vtki.vtkStructuredGrid, vtki.vtkRectilinearGrid)): 102 gf = vtki.new("GeometryFilter") 103 gf.SetInputData(inputobj) 104 gf.Update() 105 self.dataset = gf.GetOutput() 106 107 elif "meshlab" in str(type(inputobj)): 108 self.dataset = vedo.utils.meshlab2vedo(inputobj).dataset 109 110 elif "meshlib" in str(type(inputobj)): 111 import meshlib.mrmeshnumpy as mrmeshnumpy 112 self.dataset = buildPolyData( 113 mrmeshnumpy.getNumpyVerts(inputobj), 114 mrmeshnumpy.getNumpyFaces(inputobj.topology), 115 ) 116 117 elif "trimesh" in str(type(inputobj)): 118 self.dataset = vedo.utils.trimesh2vedo(inputobj).dataset 119 120 elif "meshio" in str(type(inputobj)): 121 # self.dataset = vedo.utils.meshio2vedo(inputobj) ##TODO 122 if len(inputobj.cells) > 0: 123 mcells = [] 124 for cellblock in inputobj.cells: 125 if cellblock.type in ("triangle", "quad"): 126 mcells += cellblock.data.tolist() 127 self.dataset = buildPolyData(inputobj.points, mcells) 128 else: 129 self.dataset = buildPolyData(inputobj.points, None) 130 # add arrays: 131 try: 132 if len(inputobj.point_data) > 0: 133 for k in inputobj.point_data.keys(): 134 vdata = numpy2vtk(inputobj.point_data[k]) 135 vdata.SetName(str(k)) 136 self.dataset.GetPointData().AddArray(vdata) 137 except AssertionError: 138 print("Could not add meshio point data, skip.") 139 140 else: 141 try: 142 gf = vtki.new("GeometryFilter") 143 gf.SetInputData(inputobj) 144 gf.Update() 145 self.dataset = gf.GetOutput() 146 except: 147 vedo.logger.error(f"cannot build mesh from type {type(inputobj)}") 148 raise RuntimeError() 149 150 self.mapper.SetInputData(self.dataset) 151 self.actor.SetMapper(self.mapper) 152 153 self.properties.SetInterpolationToPhong() 154 self.properties.SetColor(get_color(c)) 155 156 if alpha is not None: 157 self.properties.SetOpacity(alpha) 158 159 self.mapper.SetInterpolateScalarsBeforeMapping( 160 vedo.settings.interpolate_scalars_before_mapping 161 ) 162 163 if vedo.settings.use_polygon_offset: 164 self.mapper.SetResolveCoincidentTopologyToPolygonOffset() 165 pof = vedo.settings.polygon_offset_factor 166 pou = vedo.settings.polygon_offset_units 167 self.mapper.SetResolveCoincidentTopologyPolygonOffsetParameters(pof, pou) 168 169 n = self.dataset.GetNumberOfPoints() 170 self.pipeline = OperationNode(self, comment=f"#pts {n}") 171 172 def _repr_html_(self): 173 """ 174 HTML representation of the Mesh object for Jupyter Notebooks. 175 176 Returns: 177 HTML text with the image and some properties. 178 """ 179 import io 180 import base64 181 from PIL import Image 182 183 library_name = "vedo.mesh.Mesh" 184 help_url = "https://vedo.embl.es/docs/vedo/mesh.html#Mesh" 185 186 arr = self.thumbnail() 187 im = Image.fromarray(arr) 188 buffered = io.BytesIO() 189 im.save(buffered, format="PNG", quality=100) 190 encoded = base64.b64encode(buffered.getvalue()).decode("utf-8") 191 url = "data:image/png;base64," + encoded 192 image = f"<img src='{url}'></img>" 193 194 bounds = "<br/>".join( 195 [ 196 precision(min_x, 4) + " ... " + precision(max_x, 4) 197 for min_x, max_x in zip(self.bounds()[::2], self.bounds()[1::2]) 198 ] 199 ) 200 average_size = "{size:.3f}".format(size=self.average_size()) 201 202 help_text = "" 203 if self.name: 204 help_text += f"<b> {self.name}:   </b>" 205 help_text += '<b><a href="' + help_url + '" target="_blank">' + library_name + "</a></b>" 206 if self.filename: 207 dots = "" 208 if len(self.filename) > 30: 209 dots = "..." 210 help_text += f"<br/><code><i>({dots}{self.filename[-30:]})</i></code>" 211 212 pdata = "" 213 if self.dataset.GetPointData().GetScalars(): 214 if self.dataset.GetPointData().GetScalars().GetName(): 215 name = self.dataset.GetPointData().GetScalars().GetName() 216 pdata = "<tr><td><b> point data array </b></td><td>" + name + "</td></tr>" 217 218 cdata = "" 219 if self.dataset.GetCellData().GetScalars(): 220 if self.dataset.GetCellData().GetScalars().GetName(): 221 name = self.dataset.GetCellData().GetScalars().GetName() 222 cdata = "<tr><td><b> cell data array </b></td><td>" + name + "</td></tr>" 223 224 allt = [ 225 "<table>", 226 "<tr>", 227 "<td>", 228 image, 229 "</td>", 230 "<td style='text-align: center; vertical-align: center;'><br/>", 231 help_text, 232 "<table>", 233 "<tr><td><b> bounds </b> <br/> (x/y/z) </td><td>" + str(bounds) + "</td></tr>", 234 "<tr><td><b> center of mass </b></td><td>" 235 + precision(self.center_of_mass(), 3) 236 + "</td></tr>", 237 "<tr><td><b> average size </b></td><td>" + str(average_size) + "</td></tr>", 238 "<tr><td><b> nr. points / faces </b></td><td>" 239 + str(self.npoints) 240 + " / " 241 + str(self.ncells) 242 + "</td></tr>", 243 pdata, 244 cdata, 245 "</table>", 246 "</table>", 247 ] 248 return "\n".join(allt) 249 250 def faces(self, ids=()): 251 """DEPRECATED. Use property `mesh.cells` instead.""" 252 vedo.printc("WARNING: use property mesh.cells instead of mesh.faces()",c='y') 253 return self.cells 254 255 @property 256 def edges(self): 257 """Return an array containing the edges connectivity.""" 258 extractEdges = vtki.new("ExtractEdges") 259 extractEdges.SetInputData(self.dataset) 260 # eed.UseAllPointsOn() 261 extractEdges.Update() 262 lpoly = extractEdges.GetOutput() 263 264 arr1d = vtk2numpy(lpoly.GetLines().GetData()) 265 # [nids1, id0 ... idn, niids2, id0 ... idm, etc]. 266 267 i = 0 268 conn = [] 269 n = len(arr1d) 270 for _ in range(n): 271 cell = [arr1d[i + k + 1] for k in range(arr1d[i])] 272 conn.append(cell) 273 i += arr1d[i] + 1 274 if i >= n: 275 break 276 return conn # cannot always make a numpy array of it! 277 278 @property 279 def cell_normals(self): 280 """ 281 Retrieve face normals as a numpy array. 282 Check out also `compute_normals(cells=True)` and `compute_normals_with_pca()`. 283 """ 284 vtknormals = self.dataset.GetCellData().GetNormals() 285 numpy_normals = vtk2numpy(vtknormals) 286 if len(numpy_normals) == 0 and len(self.cells) != 0: 287 raise ValueError("VTK failed to return any normal vectors. You may need to call `Mesh.compute_normals()` before accessing `Mesh.cell_normals`.") 288 return numpy_normals 289 290 def compute_normals(self, points=True, cells=True, feature_angle=None, consistency=True) -> Self: 291 """ 292 Compute cell and vertex normals for the mesh. 293 294 Arguments: 295 points : (bool) 296 do the computation for the vertices too 297 cells : (bool) 298 do the computation for the cells too 299 feature_angle : (float) 300 specify the angle that defines a sharp edge. 301 If the difference in angle across neighboring polygons is greater than this value, 302 the shared edge is considered "sharp" and it is split. 303 consistency : (bool) 304 turn on/off the enforcement of consistent polygon ordering. 305 306 .. warning:: 307 If `feature_angle` is set then the Mesh can be modified, and it 308 can have a different nr. of vertices from the original. 309 310 Note that the appearance of the mesh may change if the normals are computed, 311 as shading is automatically enabled when such information is present. 312 Use `mesh.flat()` to avoid smoothing effects. 313 """ 314 pdnorm = vtki.new("PolyDataNormals") 315 pdnorm.SetInputData(self.dataset) 316 pdnorm.SetComputePointNormals(points) 317 pdnorm.SetComputeCellNormals(cells) 318 pdnorm.SetConsistency(consistency) 319 pdnorm.FlipNormalsOff() 320 if feature_angle: 321 pdnorm.SetSplitting(True) 322 pdnorm.SetFeatureAngle(feature_angle) 323 else: 324 pdnorm.SetSplitting(False) 325 pdnorm.Update() 326 out = pdnorm.GetOutput() 327 self._update(out, reset_locators=False) 328 return self 329 330 def reverse(self, cells=True, normals=False) -> Self: 331 """ 332 Reverse the order of polygonal cells 333 and/or reverse the direction of point and cell normals. 334 335 Two flags are used to control these operations: 336 - `cells=True` reverses the order of the indices in the cell connectivity list. 337 If cell is a list of IDs only those cells will be reversed. 338 - `normals=True` reverses the normals by multiplying the normal vector by -1 339 (both point and cell normals, if present). 340 """ 341 poly = self.dataset 342 343 if is_sequence(cells): 344 for cell in cells: 345 poly.ReverseCell(cell) 346 poly.GetCellData().Modified() 347 return self ############## 348 349 rev = vtki.new("ReverseSense") 350 if cells: 351 rev.ReverseCellsOn() 352 else: 353 rev.ReverseCellsOff() 354 if normals: 355 rev.ReverseNormalsOn() 356 else: 357 rev.ReverseNormalsOff() 358 rev.SetInputData(poly) 359 rev.Update() 360 self._update(rev.GetOutput(), reset_locators=False) 361 self.pipeline = OperationNode("reverse", parents=[self]) 362 return self 363 364 def volume(self) -> float: 365 """ 366 Compute the volume occupied by mesh. 367 The mesh must be triangular for this to work. 368 To triangulate a mesh use `mesh.triangulate()`. 369 """ 370 mass = vtki.new("MassProperties") 371 mass.SetGlobalWarningDisplay(0) 372 mass.SetInputData(self.dataset) 373 mass.Update() 374 mass.SetGlobalWarningDisplay(1) 375 return mass.GetVolume() 376 377 def area(self) -> float: 378 """ 379 Compute the surface area of the mesh. 380 The mesh must be triangular for this to work. 381 To triangulate a mesh use `mesh.triangulate()`. 382 """ 383 mass = vtki.new("MassProperties") 384 mass.SetGlobalWarningDisplay(0) 385 mass.SetInputData(self.dataset) 386 mass.Update() 387 mass.SetGlobalWarningDisplay(1) 388 return mass.GetSurfaceArea() 389 390 def is_closed(self) -> bool: 391 """ 392 Return `True` if the mesh is watertight. 393 Note that if the mesh contains coincident points the result may be flase. 394 Use in this case `mesh.clean()` to merge coincident points. 395 """ 396 fe = vtki.new("FeatureEdges") 397 fe.BoundaryEdgesOn() 398 fe.FeatureEdgesOff() 399 fe.NonManifoldEdgesOn() 400 fe.SetInputData(self.dataset) 401 fe.Update() 402 ne = fe.GetOutput().GetNumberOfCells() 403 return not bool(ne) 404 405 def is_manifold(self) -> bool: 406 """Return `True` if the mesh is manifold.""" 407 fe = vtki.new("FeatureEdges") 408 fe.BoundaryEdgesOff() 409 fe.FeatureEdgesOff() 410 fe.NonManifoldEdgesOn() 411 fe.SetInputData(self.dataset) 412 fe.Update() 413 ne = fe.GetOutput().GetNumberOfCells() 414 return not bool(ne) 415 416 def non_manifold_faces(self, remove=True, tol="auto") -> Self: 417 """ 418 Detect and (try to) remove non-manifold faces of a triangular mesh: 419 420 - set `remove` to `False` to mark cells without removing them. 421 - set `tol=0` for zero-tolerance, the result will be manifold but with holes. 422 - set `tol>0` to cut off non-manifold faces, and try to recover the good ones. 423 - set `tol="auto"` to make an automatic choice of the tolerance. 424 """ 425 # mark original point and cell ids 426 self.add_ids() 427 toremove = self.boundaries( 428 boundary_edges=False, 429 non_manifold_edges=True, 430 cell_edge=True, 431 return_cell_ids=True, 432 ) 433 if len(toremove) == 0: # type: ignore 434 return self 435 436 points = self.vertices 437 faces = self.cells 438 centers = self.cell_centers 439 440 copy = self.clone() 441 copy.delete_cells(toremove).clean() 442 copy.compute_normals(cells=False) 443 normals = copy.vertex_normals 444 deltas, deltas_i = [], [] 445 446 for i in vedo.utils.progressbar(toremove, delay=3, title="recover faces"): 447 pids = copy.closest_point(centers[i], n=3, return_point_id=True) 448 norms = normals[pids] 449 n = np.mean(norms, axis=0) 450 dn = np.linalg.norm(n) 451 if not dn: 452 continue 453 n = n / dn 454 455 p0, p1, p2 = points[faces[i]][:3] 456 v = np.cross(p1 - p0, p2 - p0) 457 lv = np.linalg.norm(v) 458 if not lv: 459 continue 460 v = v / lv 461 462 cosa = 1 - np.dot(n, v) 463 deltas.append(cosa) 464 deltas_i.append(i) 465 466 recover = [] 467 if len(deltas) > 0: 468 mean_delta = np.mean(deltas) 469 err_delta = np.std(deltas) 470 txt = "" 471 if tol == "auto": # automatic choice 472 tol = mean_delta / 5 473 txt = f"\n Automatic tol. : {tol: .4f}" 474 for i, cosa in zip(deltas_i, deltas): 475 if cosa < tol: 476 recover.append(i) 477 478 vedo.logger.info( 479 f"\n --------- Non manifold faces ---------" 480 f"\n Average tol. : {mean_delta: .4f} +- {err_delta: .4f}{txt}" 481 f"\n Removed faces : {len(toremove)}" # type: ignore 482 f"\n Recovered faces: {len(recover)}" 483 ) 484 485 toremove = list(set(toremove) - set(recover)) # type: ignore 486 487 if not remove: 488 mark = np.zeros(self.ncells, dtype=np.uint8) 489 mark[recover] = 1 490 mark[toremove] = 2 491 self.celldata["NonManifoldCell"] = mark 492 else: 493 self.delete_cells(toremove) # type: ignore 494 495 self.pipeline = OperationNode( 496 "non_manifold_faces", 497 parents=[self], 498 comment=f"#cells {self.dataset.GetNumberOfCells()}", 499 ) 500 return self 501 502 503 def euler_characteristic(self) -> int: 504 """ 505 Compute the Euler characteristic of the mesh. 506 The Euler characteristic is a topological invariant for surfaces. 507 """ 508 return self.npoints - len(self.edges) + self.ncells 509 510 def genus(self) -> int: 511 """ 512 Compute the genus of the mesh. 513 The genus is a topological invariant for surfaces. 514 """ 515 nb = len(self.boundaries().split()) - 1 516 return (2 - self.euler_characteristic() - nb ) / 2 517 518 def to_reeb_graph(self, field_id=0): 519 """ 520 Convert the mesh into a Reeb graph. 521 The Reeb graph is a topological structure that captures the evolution 522 of the level sets of a scalar field. 523 524 Arguments: 525 field_id : (int) 526 the id of the scalar field to use. 527 528 Example: 529 ```python 530 from vedo import * 531 mesh = Mesh("https://discourse.paraview.org/uploads/short-url/qVuZ1fiRjwhE1qYtgGE2HGXybgo.stl") 532 mesh.rotate_x(10).rotate_y(15).alpha(0.5) 533 mesh.pointdata["scalars"] = mesh.vertices[:, 2] 534 535 printc("is_closed :", mesh.is_closed()) 536 printc("is_manifold:", mesh.is_manifold()) 537 printc("euler_char :", mesh.euler_characteristic()) 538 printc("genus :", mesh.genus()) 539 540 reeb = mesh.to_reeb_graph() 541 ids = reeb[0].pointdata["Vertex Ids"] 542 pts = Points(mesh.vertices[ids], r=10) 543 544 show([[mesh, pts], reeb], N=2, sharecam=False) 545 ``` 546 """ 547 rg = vtki.new("PolyDataToReebGraphFilter") 548 rg.SetInputData(self.dataset) 549 rg.SetFieldId(field_id) 550 rg.Update() 551 gr = vedo.pyplot.DirectedGraph() 552 gr.mdg = rg.GetOutput() 553 gr.build() 554 return gr 555 556 557 def shrink(self, fraction=0.85) -> Self: 558 """ 559 Shrink the triangle polydata in the representation of the input mesh. 560 561 Examples: 562 - [shrink.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/shrink.py) 563 564  565 """ 566 # Overriding base class method core.shrink() 567 shrink = vtki.new("ShrinkPolyData") 568 shrink.SetInputData(self.dataset) 569 shrink.SetShrinkFactor(fraction) 570 shrink.Update() 571 self._update(shrink.GetOutput()) 572 self.pipeline = OperationNode("shrink", parents=[self]) 573 return self 574 575 def cap(self, return_cap=False) -> Self: 576 """ 577 Generate a "cap" on a clipped mesh, or caps sharp edges. 578 579 Examples: 580 - [cut_and_cap.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/cut_and_cap.py) 581 582  583 584 See also: `join()`, `join_segments()`, `slice()`. 585 """ 586 fe = vtki.new("FeatureEdges") 587 fe.SetInputData(self.dataset) 588 fe.BoundaryEdgesOn() 589 fe.FeatureEdgesOff() 590 fe.NonManifoldEdgesOff() 591 fe.ManifoldEdgesOff() 592 fe.Update() 593 594 stripper = vtki.new("Stripper") 595 stripper.SetInputData(fe.GetOutput()) 596 stripper.JoinContiguousSegmentsOn() 597 stripper.Update() 598 599 boundary_poly = vtki.vtkPolyData() 600 boundary_poly.SetPoints(stripper.GetOutput().GetPoints()) 601 boundary_poly.SetPolys(stripper.GetOutput().GetLines()) 602 603 rev = vtki.new("ReverseSense") 604 rev.ReverseCellsOn() 605 rev.SetInputData(boundary_poly) 606 rev.Update() 607 608 tf = vtki.new("TriangleFilter") 609 tf.SetInputData(rev.GetOutput()) 610 tf.Update() 611 612 if return_cap: 613 m = Mesh(tf.GetOutput()) 614 m.pipeline = OperationNode( 615 "cap", parents=[self], comment=f"#pts {m.dataset.GetNumberOfPoints()}" 616 ) 617 m.name = "MeshCap" 618 return m 619 620 polyapp = vtki.new("AppendPolyData") 621 polyapp.AddInputData(self.dataset) 622 polyapp.AddInputData(tf.GetOutput()) 623 polyapp.Update() 624 625 self._update(polyapp.GetOutput()) 626 self.clean() 627 628 self.pipeline = OperationNode( 629 "capped", parents=[self], comment=f"#pts {self.dataset.GetNumberOfPoints()}" 630 ) 631 return self 632 633 def join(self, polys=True, reset=False) -> Self: 634 """ 635 Generate triangle strips and/or polylines from 636 input polygons, triangle strips, and lines. 637 638 Input polygons are assembled into triangle strips only if they are triangles; 639 other types of polygons are passed through to the output and not stripped. 640 Use mesh.triangulate() to triangulate non-triangular polygons prior to running 641 this filter if you need to strip all the data. 642 643 Also note that if triangle strips or polylines are present in the input 644 they are passed through and not joined nor extended. 645 If you wish to strip these use mesh.triangulate() to fragment the input 646 into triangles and lines prior to applying join(). 647 648 Arguments: 649 polys : (bool) 650 polygonal segments will be joined if they are contiguous 651 reset : (bool) 652 reset points ordering 653 654 Warning: 655 If triangle strips or polylines exist in the input data 656 they will be passed through to the output data. 657 This filter will only construct triangle strips if triangle polygons 658 are available; and will only construct polylines if lines are available. 659 660 Example: 661 ```python 662 from vedo import * 663 c1 = Cylinder(pos=(0,0,0), r=2, height=3, axis=(1,.0,0), alpha=.1).triangulate() 664 c2 = Cylinder(pos=(0,0,2), r=1, height=2, axis=(0,.3,1), alpha=.1).triangulate() 665 intersect = c1.intersect_with(c2).join(reset=True) 666 spline = Spline(intersect).c('blue').lw(5) 667 show(c1, c2, spline, intersect.labels('id'), axes=1).close() 668 ``` 669  670 """ 671 sf = vtki.new("Stripper") 672 sf.SetPassThroughCellIds(True) 673 sf.SetPassThroughPointIds(True) 674 sf.SetJoinContiguousSegments(polys) 675 sf.SetInputData(self.dataset) 676 sf.Update() 677 if reset: 678 poly = sf.GetOutput() 679 cpd = vtki.new("CleanPolyData") 680 cpd.PointMergingOn() 681 cpd.ConvertLinesToPointsOn() 682 cpd.ConvertPolysToLinesOn() 683 cpd.ConvertStripsToPolysOn() 684 cpd.SetInputData(poly) 685 cpd.Update() 686 poly = cpd.GetOutput() 687 vpts = poly.GetCell(0).GetPoints().GetData() 688 poly.GetPoints().SetData(vpts) 689 else: 690 poly = sf.GetOutput() 691 692 self._update(poly) 693 694 self.pipeline = OperationNode( 695 "join", parents=[self], comment=f"#pts {self.dataset.GetNumberOfPoints()}" 696 ) 697 return self 698 699 def join_segments(self, closed=True, tol=1e-03) -> list: 700 """ 701 Join line segments into contiguous lines. 702 Useful to call with `triangulate()` method. 703 704 Returns: 705 list of `shapes.Lines` 706 707 Example: 708 ```python 709 from vedo import * 710 msh = Torus().alpha(0.1).wireframe() 711 intersection = msh.intersect_with_plane(normal=[1,1,1]).c('purple5') 712 slices = [s.triangulate() for s in intersection.join_segments()] 713 show(msh, intersection, merge(slices), axes=1, viewup='z') 714 ``` 715  716 """ 717 vlines = [] 718 for ipiece, outline in enumerate(self.split(must_share_edge=False)): # type: ignore 719 720 outline.clean() 721 pts = outline.vertices 722 if len(pts) < 3: 723 continue 724 avesize = outline.average_size() 725 lines = outline.lines 726 # print("---lines", lines, "in piece", ipiece) 727 tol = avesize / pts.shape[0] * tol 728 729 k = 0 730 joinedpts = [pts[k]] 731 for _ in range(len(pts)): 732 pk = pts[k] 733 for j, line in enumerate(lines): 734 735 id0, id1 = line[0], line[-1] 736 p0, p1 = pts[id0], pts[id1] 737 738 if np.linalg.norm(p0 - pk) < tol: 739 n = len(line) 740 for m in range(1, n): 741 joinedpts.append(pts[line[m]]) 742 # joinedpts.append(p1) 743 k = id1 744 lines.pop(j) 745 break 746 747 elif np.linalg.norm(p1 - pk) < tol: 748 n = len(line) 749 for m in reversed(range(0, n - 1)): 750 joinedpts.append(pts[line[m]]) 751 # joinedpts.append(p0) 752 k = id0 753 lines.pop(j) 754 break 755 756 if len(joinedpts) > 1: 757 newline = vedo.shapes.Line(joinedpts, closed=closed) 758 newline.clean() 759 newline.actor.SetProperty(self.properties) 760 newline.properties = self.properties 761 newline.pipeline = OperationNode( 762 "join_segments", 763 parents=[self], 764 comment=f"#pts {newline.dataset.GetNumberOfPoints()}", 765 ) 766 vlines.append(newline) 767 768 return vlines 769 770 def join_with_strips(self, b1, closed=True) -> Self: 771 """ 772 Join booundary lines by creating a triangle strip between them. 773 774 Example: 775 ```python 776 from vedo import * 777 m1 = Cylinder(cap=False).boundaries() 778 m2 = Cylinder(cap=False).boundaries().pos(0.2,0,1) 779 strips = m1.join_with_strips(m2) 780 show(m1, m2, strips, axes=1).close() 781 ``` 782 """ 783 b0 = self.clone().join() 784 b1 = b1.clone().join() 785 786 vertices0 = b0.vertices.tolist() 787 vertices1 = b1.vertices.tolist() 788 789 lines0 = b0.lines 790 lines1 = b1.lines 791 m = len(lines0) 792 assert m == len(lines1), ( 793 "lines must have the same number of points\n" 794 f"line has {m} points in b0 and {len(lines1)} in b1" 795 ) 796 797 strips = [] 798 points: List[Any] = [] 799 800 for j in range(m): 801 802 ids0j = list(lines0[j]) 803 ids1j = list(lines1[j]) 804 805 n = len(ids0j) 806 assert n == len(ids1j), ( 807 "lines must have the same number of points\n" 808 f"line {j} has {n} points in b0 and {len(ids1j)} in b1" 809 ) 810 811 if closed: 812 ids0j.append(ids0j[0]) 813 ids1j.append(ids1j[0]) 814 vertices0.append(vertices0[ids0j[0]]) 815 vertices1.append(vertices1[ids1j[0]]) 816 n = n + 1 817 818 strip = [] # create a triangle strip 819 npt = len(points) 820 for ipt in range(n): 821 points.append(vertices0[ids0j[ipt]]) 822 points.append(vertices1[ids1j[ipt]]) 823 824 strip = list(range(npt, npt + 2*n)) 825 strips.append(strip) 826 827 return Mesh([points, [], [], strips], c="k6") 828 829 def split_polylines(self) -> Self: 830 """Split polylines into separate segments.""" 831 tf = vtki.new("TriangleFilter") 832 tf.SetPassLines(True) 833 tf.SetPassVerts(False) 834 tf.SetInputData(self.dataset) 835 tf.Update() 836 self._update(tf.GetOutput(), reset_locators=False) 837 self.lw(0).lighting("default").pickable() 838 self.pipeline = OperationNode( 839 "split_polylines", parents=[self], 840 comment=f"#lines {self.dataset.GetNumberOfLines()}" 841 ) 842 return self 843 844 def slice(self, origin=(0, 0, 0), normal=(1, 0, 0)) -> Self: 845 """ 846 Slice a mesh with a plane and fill the contour. 847 848 Example: 849 ```python 850 from vedo import * 851 msh = Mesh(dataurl+"bunny.obj").alpha(0.1).wireframe() 852 mslice = msh.slice(normal=[0,1,0.3], origin=[0,0.16,0]) 853 mslice.c('purple5') 854 show(msh, mslice, axes=1) 855 ``` 856  857 858 See also: `join()`, `join_segments()`, `cap()`, `cut_with_plane()`. 859 """ 860 intersection = self.intersect_with_plane(origin=origin, normal=normal) 861 slices = [s.triangulate() for s in intersection.join_segments()] 862 mslices = vedo.pointcloud.merge(slices) 863 if mslices: 864 mslices.name = "MeshSlice" 865 mslices.pipeline = OperationNode("slice", parents=[self], comment=f"normal = {normal}") 866 return mslices 867 868 def triangulate(self, verts=True, lines=True) -> Self: 869 """ 870 Converts mesh polygons into triangles. 871 872 If the input mesh is only made of 2D lines (no faces) the output will be a triangulation 873 that fills the internal area. The contours may be concave, and may even contain holes, 874 i.e. a contour may contain an internal contour winding in the opposite 875 direction to indicate that it is a hole. 876 877 Arguments: 878 verts : (bool) 879 if True, break input vertex cells into individual vertex cells (one point per cell). 880 If False, the input vertex cells will be ignored. 881 lines : (bool) 882 if True, break input polylines into line segments. 883 If False, input lines will be ignored and the output will have no lines. 884 """ 885 if self.dataset.GetNumberOfPolys() or self.dataset.GetNumberOfStrips(): 886 # print("Using vtkTriangleFilter") 887 tf = vtki.new("TriangleFilter") 888 tf.SetPassLines(lines) 889 tf.SetPassVerts(verts) 890 891 elif self.dataset.GetNumberOfLines(): 892 # print("Using vtkContourTriangulator") 893 tf = vtki.new("ContourTriangulator") 894 tf.TriangulationErrorDisplayOn() 895 896 else: 897 vedo.logger.debug("input in triangulate() seems to be void! Skip.") 898 return self 899 900 tf.SetInputData(self.dataset) 901 tf.Update() 902 self._update(tf.GetOutput(), reset_locators=False) 903 self.lw(0).lighting("default").pickable() 904 905 self.pipeline = OperationNode( 906 "triangulate", parents=[self], comment=f"#cells {self.dataset.GetNumberOfCells()}" 907 ) 908 return self 909 910 def compute_cell_vertex_count(self) -> Self: 911 """ 912 Add to this mesh a cell data array containing the nr of vertices that a polygonal face has. 913 """ 914 csf = vtki.new("CellSizeFilter") 915 csf.SetInputData(self.dataset) 916 csf.SetComputeArea(False) 917 csf.SetComputeVolume(False) 918 csf.SetComputeLength(False) 919 csf.SetComputeVertexCount(True) 920 csf.SetVertexCountArrayName("VertexCount") 921 csf.Update() 922 self.dataset.GetCellData().AddArray( 923 csf.GetOutput().GetCellData().GetArray("VertexCount") 924 ) 925 return self 926 927 def compute_quality(self, metric=6) -> Self: 928 """ 929 Calculate metrics of quality for the elements of a triangular mesh. 930 This method adds to the mesh a cell array named "Quality". 931 See class 932 [vtkMeshQuality](https://vtk.org/doc/nightly/html/classvtkMeshQuality.html). 933 934 Arguments: 935 metric : (int) 936 type of available estimators are: 937 - EDGE RATIO, 0 938 - ASPECT RATIO, 1 939 - RADIUS RATIO, 2 940 - ASPECT FROBENIUS, 3 941 - MED ASPECT FROBENIUS, 4 942 - MAX ASPECT FROBENIUS, 5 943 - MIN_ANGLE, 6 944 - COLLAPSE RATIO, 7 945 - MAX ANGLE, 8 946 - CONDITION, 9 947 - SCALED JACOBIAN, 10 948 - SHEAR, 11 949 - RELATIVE SIZE SQUARED, 12 950 - SHAPE, 13 951 - SHAPE AND SIZE, 14 952 - DISTORTION, 15 953 - MAX EDGE RATIO, 16 954 - SKEW, 17 955 - TAPER, 18 956 - VOLUME, 19 957 - STRETCH, 20 958 - DIAGONAL, 21 959 - DIMENSION, 22 960 - ODDY, 23 961 - SHEAR AND SIZE, 24 962 - JACOBIAN, 25 963 - WARPAGE, 26 964 - ASPECT GAMMA, 27 965 - AREA, 28 966 - ASPECT BETA, 29 967 968 Examples: 969 - [meshquality.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/meshquality.py) 970 971  972 """ 973 qf = vtki.new("MeshQuality") 974 qf.SetInputData(self.dataset) 975 qf.SetTriangleQualityMeasure(metric) 976 qf.SaveCellQualityOn() 977 qf.Update() 978 self._update(qf.GetOutput(), reset_locators=False) 979 self.mapper.SetScalarModeToUseCellData() 980 self.pipeline = OperationNode("compute_quality", parents=[self]) 981 return self 982 983 def count_vertices(self) -> np.ndarray: 984 """Count the number of vertices each cell has and return it as a numpy array""" 985 vc = vtki.new("CountVertices") 986 vc.SetInputData(self.dataset) 987 vc.SetOutputArrayName("VertexCount") 988 vc.Update() 989 varr = vc.GetOutput().GetCellData().GetArray("VertexCount") 990 return vtk2numpy(varr) 991 992 def check_validity(self, tol=0) -> np.ndarray: 993 """ 994 Return a numpy array of possible problematic faces following this convention: 995 - Valid = 0 996 - WrongNumberOfPoints = 1 997 - IntersectingEdges = 2 998 - IntersectingFaces = 4 999 - NoncontiguousEdges = 8 1000 - Nonconvex = 10 1001 - OrientedIncorrectly = 20 1002 1003 Arguments: 1004 tol : (float) 1005 value is used as an epsilon for floating point 1006 equality checks throughout the cell checking process. 1007 """ 1008 vald = vtki.new("CellValidator") 1009 if tol: 1010 vald.SetTolerance(tol) 1011 vald.SetInputData(self.dataset) 1012 vald.Update() 1013 varr = vald.GetOutput().GetCellData().GetArray("ValidityState") 1014 return vtk2numpy(varr) 1015 1016 def compute_curvature(self, method=0) -> Self: 1017 """ 1018 Add scalars to `Mesh` that contains the curvature calculated in three different ways. 1019 1020 Variable `method` can be: 1021 - 0 = gaussian 1022 - 1 = mean curvature 1023 - 2 = max curvature 1024 - 3 = min curvature 1025 1026 Example: 1027 ```python 1028 from vedo import Torus 1029 Torus().compute_curvature().add_scalarbar().show().close() 1030 ``` 1031  1032 """ 1033 curve = vtki.new("Curvatures") 1034 curve.SetInputData(self.dataset) 1035 curve.SetCurvatureType(method) 1036 curve.Update() 1037 self._update(curve.GetOutput(), reset_locators=False) 1038 self.mapper.ScalarVisibilityOn() 1039 return self 1040 1041 def compute_elevation(self, low=(0, 0, 0), high=(0, 0, 1), vrange=(0, 1)) -> Self: 1042 """ 1043 Add to `Mesh` a scalar array that contains distance along a specified direction. 1044 1045 Arguments: 1046 low : (list) 1047 one end of the line (small scalar values) 1048 high : (list) 1049 other end of the line (large scalar values) 1050 vrange : (list) 1051 set the range of the scalar 1052 1053 Example: 1054 ```python 1055 from vedo import Sphere 1056 s = Sphere().compute_elevation(low=(0,0,0), high=(1,1,1)) 1057 s.add_scalarbar().show(axes=1).close() 1058 ``` 1059  1060 """ 1061 ef = vtki.new("ElevationFilter") 1062 ef.SetInputData(self.dataset) 1063 ef.SetLowPoint(low) 1064 ef.SetHighPoint(high) 1065 ef.SetScalarRange(vrange) 1066 ef.Update() 1067 self._update(ef.GetOutput(), reset_locators=False) 1068 self.mapper.ScalarVisibilityOn() 1069 return self 1070 1071 1072 def laplacian_diffusion(self, array_name, dt, num_steps) -> Self: 1073 """ 1074 Apply a diffusion process to a scalar array defined on the points of a mesh. 1075 1076 Arguments: 1077 array_name : (str) 1078 name of the array to diffuse. 1079 dt : (float) 1080 time step. 1081 num_steps : (int) 1082 number of iterations. 1083 """ 1084 try: 1085 import scipy.sparse 1086 import scipy.sparse.linalg 1087 except ImportError: 1088 vedo.logger.error("scipy not found. Cannot run laplacian_diffusion()") 1089 return self 1090 1091 def build_laplacian(): 1092 rows = [] 1093 cols = [] 1094 data = [] 1095 n_points = points.shape[0] 1096 avg_area = np.mean(areas) * 10000 1097 # print("avg_area", avg_area) 1098 1099 for triangle in cells: 1100 for i in range(3): 1101 for j in range(i + 1, 3): 1102 u = triangle[i] 1103 v = triangle[j] 1104 rows.append(u) 1105 cols.append(v) 1106 rows.append(v) 1107 cols.append(u) 1108 data.append(-1/avg_area) 1109 data.append(-1/avg_area) 1110 1111 L = scipy.sparse.coo_matrix( 1112 (data, (rows, cols)), shape=(n_points, n_points) 1113 ).tocsc() 1114 1115 degree = -np.array(L.sum(axis=1)).flatten() # adjust the diagonal 1116 # print("degree", degree) 1117 L.setdiag(degree) 1118 return L 1119 1120 def _diffuse(u0, L, dt, num_steps): 1121 # mean_area = np.mean(areas) * 10000 1122 # print("mean_area", mean_area) 1123 mean_area = 1 1124 I = scipy.sparse.eye(L.shape[0], format="csc") 1125 A = I - (dt/mean_area) * L 1126 u = u0 1127 for _ in range(int(num_steps)): 1128 u = A.dot(u) 1129 return u 1130 1131 self.compute_cell_size() 1132 areas = self.celldata["Area"] 1133 points = self.vertices 1134 cells = self.cells 1135 u0 = self.pointdata[array_name] 1136 1137 # Simulate diffusion 1138 L = build_laplacian() 1139 u = _diffuse(u0, L, dt, num_steps) 1140 self.pointdata[array_name] = u 1141 return self 1142 1143 1144 def subdivide(self, n=1, method=0, mel=None) -> Self: 1145 """ 1146 Increase the number of vertices of a surface mesh. 1147 1148 Arguments: 1149 n : (int) 1150 number of subdivisions. 1151 method : (int) 1152 Loop(0), Linear(1), Adaptive(2), Butterfly(3), Centroid(4) 1153 mel : (float) 1154 Maximum Edge Length (applicable to Adaptive method only). 1155 """ 1156 triangles = vtki.new("TriangleFilter") 1157 triangles.SetInputData(self.dataset) 1158 triangles.Update() 1159 tri_mesh = triangles.GetOutput() 1160 if method == 0: 1161 sdf = vtki.new("LoopSubdivisionFilter") 1162 elif method == 1: 1163 sdf = vtki.new("LinearSubdivisionFilter") 1164 elif method == 2: 1165 sdf = vtki.new("AdaptiveSubdivisionFilter") 1166 if mel is None: 1167 mel = self.diagonal_size() / np.sqrt(self.dataset.GetNumberOfPoints()) / n 1168 sdf.SetMaximumEdgeLength(mel) 1169 elif method == 3: 1170 sdf = vtki.new("ButterflySubdivisionFilter") 1171 elif method == 4: 1172 sdf = vtki.new("DensifyPolyData") 1173 else: 1174 vedo.logger.error(f"in subdivide() unknown method {method}") 1175 raise RuntimeError() 1176 1177 if method != 2: 1178 sdf.SetNumberOfSubdivisions(n) 1179 1180 sdf.SetInputData(tri_mesh) 1181 sdf.Update() 1182 1183 self._update(sdf.GetOutput()) 1184 1185 self.pipeline = OperationNode( 1186 "subdivide", 1187 parents=[self], 1188 comment=f"#pts {self.dataset.GetNumberOfPoints()}", 1189 ) 1190 return self 1191 1192 1193 def decimate(self, fraction=0.5, n=None, preserve_volume=True, regularization=0.0) -> Self: 1194 """ 1195 Downsample the number of vertices in a mesh to `fraction`. 1196 1197 This filter preserves the `pointdata` of the input dataset. In previous versions 1198 of vedo, this decimation algorithm was referred to as quadric decimation. 1199 1200 Arguments: 1201 fraction : (float) 1202 the desired target of reduction. 1203 n : (int) 1204 the desired number of final points 1205 (`fraction` is recalculated based on it). 1206 preserve_volume : (bool) 1207 Decide whether to activate volume preservation which greatly 1208 reduces errors in triangle normal direction. 1209 regularization : (float) 1210 regularize the point finding algorithm so as to have better quality 1211 mesh elements at the cost of a slightly lower precision on the 1212 geometry potentially (mostly at sharp edges). 1213 Can be useful for decimating meshes that have been triangulated on noisy data. 1214 1215 Note: 1216 Setting `fraction=0.1` leaves 10% of the original number of vertices. 1217 Internally the VTK class 1218 [vtkQuadricDecimation](https://vtk.org/doc/nightly/html/classvtkQuadricDecimation.html) 1219 is used for this operation. 1220 1221 See also: `decimate_binned()` and `decimate_pro()`. 1222 """ 1223 poly = self.dataset 1224 if n: # N = desired number of points 1225 npt = poly.GetNumberOfPoints() 1226 fraction = n / npt 1227 if fraction >= 1: 1228 return self 1229 1230 decimate = vtki.new("QuadricDecimation") 1231 decimate.SetVolumePreservation(preserve_volume) 1232 # decimate.AttributeErrorMetricOn() 1233 if regularization: 1234 decimate.SetRegularize(True) 1235 decimate.SetRegularization(regularization) 1236 1237 try: 1238 decimate.MapPointDataOn() 1239 except AttributeError: 1240 pass 1241 1242 decimate.SetTargetReduction(1 - fraction) 1243 decimate.SetInputData(poly) 1244 decimate.Update() 1245 1246 self._update(decimate.GetOutput()) 1247 self.metadata["decimate_actual_fraction"] = 1 - decimate.GetActualReduction() 1248 1249 self.pipeline = OperationNode( 1250 "decimate", 1251 parents=[self], 1252 comment=f"#pts {self.dataset.GetNumberOfPoints()}", 1253 ) 1254 return self 1255 1256 def decimate_pro( 1257 self, 1258 fraction=0.5, 1259 n=None, 1260 preserve_topology=True, 1261 preserve_boundaries=True, 1262 splitting=False, 1263 splitting_angle=75, 1264 feature_angle=0, 1265 inflection_point_ratio=10, 1266 vertex_degree=0, 1267 ) -> Self: 1268 """ 1269 Downsample the number of vertices in a mesh to `fraction`. 1270 1271 This filter preserves the `pointdata` of the input dataset. 1272 1273 Arguments: 1274 fraction : (float) 1275 The desired target of reduction. 1276 Setting `fraction=0.1` leaves 10% of the original number of vertices. 1277 n : (int) 1278 the desired number of final points (`fraction` is recalculated based on it). 1279 preserve_topology : (bool) 1280 If on, mesh splitting and hole elimination will not occur. 1281 This may limit the maximum reduction that may be achieved. 1282 preserve_boundaries : (bool) 1283 Turn on/off the deletion of vertices on the boundary of a mesh. 1284 Control whether mesh boundaries are preserved during decimation. 1285 feature_angle : (float) 1286 Specify the angle that defines a feature. 1287 This angle is used to define what an edge is 1288 (i.e., if the surface normal between two adjacent triangles 1289 is >= FeatureAngle, an edge exists). 1290 splitting : (bool) 1291 Turn on/off the splitting of the mesh at corners, 1292 along edges, at non-manifold points, or anywhere else a split is required. 1293 Turning splitting off will better preserve the original topology of the mesh, 1294 but you may not obtain the requested reduction. 1295 splitting_angle : (float) 1296 Specify the angle that defines a sharp edge. 1297 This angle is used to control the splitting of the mesh. 1298 A split line exists when the surface normals between two edge connected triangles 1299 are >= `splitting_angle`. 1300 inflection_point_ratio : (float) 1301 An inflection point occurs when the ratio of reduction error between two iterations 1302 is greater than or equal to the `inflection_point_ratio` value. 1303 vertex_degree : (int) 1304 If the number of triangles connected to a vertex exceeds it then the vertex will be split. 1305 1306 Note: 1307 Setting `fraction=0.1` leaves 10% of the original number of vertices 1308 1309 See also: 1310 `decimate()` and `decimate_binned()`. 1311 """ 1312 poly = self.dataset 1313 if n: # N = desired number of points 1314 npt = poly.GetNumberOfPoints() 1315 fraction = n / npt 1316 if fraction >= 1: 1317 return self 1318 1319 decimate = vtki.new("DecimatePro") 1320 decimate.SetPreserveTopology(preserve_topology) 1321 decimate.SetBoundaryVertexDeletion(preserve_boundaries) 1322 if feature_angle: 1323 decimate.SetFeatureAngle(feature_angle) 1324 decimate.SetSplitting(splitting) 1325 decimate.SetSplitAngle(splitting_angle) 1326 decimate.SetInflectionPointRatio(inflection_point_ratio) 1327 if vertex_degree: 1328 decimate.SetDegree(vertex_degree) 1329 1330 decimate.SetTargetReduction(1 - fraction) 1331 decimate.SetInputData(poly) 1332 decimate.Update() 1333 self._update(decimate.GetOutput()) 1334 1335 self.pipeline = OperationNode( 1336 "decimate_pro", 1337 parents=[self], 1338 comment=f"#pts {self.dataset.GetNumberOfPoints()}", 1339 ) 1340 return self 1341 1342 def decimate_binned(self, divisions=(), use_clustering=False) -> Self: 1343 """ 1344 Downsample the number of vertices in a mesh. 1345 1346 This filter preserves the `celldata` of the input dataset, 1347 if `use_clustering=True` also the `pointdata` will be preserved in the result. 1348 1349 Arguments: 1350 divisions : (list) 1351 number of divisions along x, y and z axes. 1352 auto_adjust : (bool) 1353 if True, the number of divisions is automatically adjusted to 1354 create more uniform cells. 1355 use_clustering : (bool) 1356 use [vtkQuadricClustering](https://vtk.org/doc/nightly/html/classvtkQuadricClustering.html) 1357 instead of 1358 [vtkBinnedDecimation](https://vtk.org/doc/nightly/html/classvtkBinnedDecimation.html). 1359 1360 See also: `decimate()` and `decimate_pro()`. 1361 """ 1362 if use_clustering: 1363 decimate = vtki.new("QuadricClustering") 1364 decimate.CopyCellDataOn() 1365 else: 1366 decimate = vtki.new("BinnedDecimation") 1367 decimate.ProducePointDataOn() 1368 decimate.ProduceCellDataOn() 1369 1370 decimate.SetInputData(self.dataset) 1371 1372 if len(divisions) == 0: 1373 decimate.SetAutoAdjustNumberOfDivisions(1) 1374 else: 1375 decimate.SetAutoAdjustNumberOfDivisions(0) 1376 decimate.SetNumberOfDivisions(divisions) 1377 decimate.Update() 1378 1379 self._update(decimate.GetOutput()) 1380 self.metadata["decimate_binned_divisions"] = decimate.GetNumberOfDivisions() 1381 self.pipeline = OperationNode( 1382 "decimate_binned", 1383 parents=[self], 1384 comment=f"#pts {self.dataset.GetNumberOfPoints()}", 1385 ) 1386 return self 1387 1388 def generate_random_points(self, n: int, min_radius=0.0) -> "Points": 1389 """ 1390 Generate `n` uniformly distributed random points 1391 inside the polygonal mesh. 1392 1393 A new point data array is added to the output points 1394 called "OriginalCellID" which contains the index of 1395 the cell ID in which the point was generated. 1396 1397 Arguments: 1398 n : (int) 1399 number of points to generate. 1400 min_radius: (float) 1401 impose a minimum distance between points. 1402 If `min_radius` is set to 0, the points are 1403 generated uniformly at random inside the mesh. 1404 If `min_radius` is set to a positive value, 1405 the points are generated uniformly at random 1406 inside the mesh, but points closer than `min_radius` 1407 to any other point are discarded. 1408 1409 Returns a `vedo.Points` object. 1410 1411 Note: 1412 Consider using `points.probe(msh)` or 1413 `points.interpolate_data_from(msh)` 1414 to interpolate existing mesh data onto the new points. 1415 1416 Example: 1417 ```python 1418 from vedo import * 1419 msh = Mesh(dataurl + "panther.stl").lw(2) 1420 pts = msh.generate_random_points(20000, min_radius=0.5) 1421 print("Original cell ids:", pts.pointdata["OriginalCellID"]) 1422 show(pts, msh, axes=1).close() 1423 ``` 1424 """ 1425 cmesh = self.clone().clean().triangulate().compute_cell_size() 1426 triangles = cmesh.cells 1427 vertices = cmesh.vertices 1428 cumul = np.cumsum(cmesh.celldata["Area"]) 1429 1430 out_pts = [] 1431 orig_cell = [] 1432 for _ in range(n): 1433 # choose a triangle based on area 1434 random_area = np.random.random() * cumul[-1] 1435 it = np.searchsorted(cumul, random_area) 1436 A, B, C = vertices[triangles[it]] 1437 # calculate the random point in the triangle 1438 r1, r2 = np.random.random(2) 1439 if r1 + r2 > 1: 1440 r1 = 1 - r1 1441 r2 = 1 - r2 1442 out_pts.append((1 - r1 - r2) * A + r1 * B + r2 * C) 1443 orig_cell.append(it) 1444 nporig_cell = np.array(orig_cell, dtype=np.uint32) 1445 1446 vpts = Points(out_pts) 1447 vpts.pointdata["OriginalCellID"] = nporig_cell 1448 1449 if min_radius > 0: 1450 vpts.subsample(min_radius, absolute=True) 1451 1452 vpts.point_size(5).color("k1") 1453 vpts.name = "RandomPoints" 1454 vpts.pipeline = OperationNode( 1455 "generate_random_points", c="#edabab", parents=[self]) 1456 return vpts 1457 1458 def delete_cells(self, ids: List[int]) -> Self: 1459 """ 1460 Remove cells from the mesh object by their ID. 1461 Points (vertices) are not removed (you may use `clean()` to remove those). 1462 """ 1463 self.dataset.BuildLinks() 1464 for cid in ids: 1465 self.dataset.DeleteCell(cid) 1466 self.dataset.RemoveDeletedCells() 1467 self.dataset.Modified() 1468 self.mapper.Modified() 1469 self.pipeline = OperationNode( 1470 "delete_cells", 1471 parents=[self], 1472 comment=f"#cells {self.dataset.GetNumberOfCells()}", 1473 ) 1474 return self 1475 1476 def delete_cells_by_point_index(self, indices: List[int]) -> Self: 1477 """ 1478 Delete a list of vertices identified by any of their vertex index. 1479 1480 See also `delete_cells()`. 1481 1482 Examples: 1483 - [delete_mesh_pts.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/delete_mesh_pts.py) 1484 1485  1486 """ 1487 cell_ids = vtki.vtkIdList() 1488 self.dataset.BuildLinks() 1489 n = 0 1490 for i in np.unique(indices): 1491 self.dataset.GetPointCells(i, cell_ids) 1492 for j in range(cell_ids.GetNumberOfIds()): 1493 self.dataset.DeleteCell(cell_ids.GetId(j)) # flag cell 1494 n += 1 1495 1496 self.dataset.RemoveDeletedCells() 1497 self.dataset.Modified() 1498 self.pipeline = OperationNode("delete_cells_by_point_index", parents=[self]) 1499 return self 1500 1501 def collapse_edges(self, distance: float, iterations=1) -> Self: 1502 """ 1503 Collapse mesh edges so that are all above `distance`. 1504 1505 Example: 1506 ```python 1507 from vedo import * 1508 np.random.seed(2) 1509 grid1 = Grid().add_gaussian_noise(0.8).triangulate().lw(1) 1510 grid1.celldata['scalar'] = grid1.cell_centers[:,1] 1511 grid2 = grid1.clone().collapse_edges(0.1) 1512 show(grid1, grid2, N=2, axes=1) 1513 ``` 1514 """ 1515 for _ in range(iterations): 1516 medges = self.edges 1517 pts = self.vertices 1518 newpts = np.array(pts) 1519 moved = [] 1520 for e in medges: 1521 if len(e) == 2: 1522 id0, id1 = e 1523 p0, p1 = pts[id0], pts[id1] 1524 if (np.linalg.norm(p1-p0) < distance 1525 and id0 not in moved 1526 and id1 not in moved 1527 ): 1528 p = (p0 + p1) / 2 1529 newpts[id0] = p 1530 newpts[id1] = p 1531 moved += [id0, id1] 1532 self.vertices = newpts 1533 cpd = vtki.new("CleanPolyData") 1534 cpd.ConvertLinesToPointsOff() 1535 cpd.ConvertPolysToLinesOff() 1536 cpd.ConvertStripsToPolysOff() 1537 cpd.SetInputData(self.dataset) 1538 cpd.Update() 1539 self._update(cpd.GetOutput()) 1540 1541 self.pipeline = OperationNode( 1542 "collapse_edges", 1543 parents=[self], 1544 comment=f"#pts {self.dataset.GetNumberOfPoints()}", 1545 ) 1546 return self 1547 1548 def adjacency_list(self) -> List[set]: 1549 """ 1550 Computes the adjacency list for mesh edge-graph. 1551 1552 Returns: 1553 a list with i-th entry being the set if indices of vertices connected by an edge to i-th vertex 1554 """ 1555 inc = [set()] * self.nvertices 1556 for cell in self.cells: 1557 nc = len(cell) 1558 if nc > 1: 1559 for i in range(nc-1): 1560 ci = cell[i] 1561 inc[ci] = inc[ci].union({cell[i-1], cell[i+1]}) 1562 return inc 1563 1564 def graph_ball(self, index, n: int) -> set: 1565 """ 1566 Computes the ball of radius `n` in the mesh' edge-graph metric centred in vertex `index`. 1567 1568 Arguments: 1569 index : (int) 1570 index of the vertex 1571 n : (int) 1572 radius in the graph metric 1573 1574 Returns: 1575 the set of indices of the vertices which are at most `n` edges from vertex `index`. 1576 """ 1577 if n == 0: 1578 return {index} 1579 else: 1580 al = self.adjacency_list() 1581 ball = {index} 1582 i = 0 1583 while i < n and len(ball) < self.nvertices: 1584 for v in ball: 1585 ball = ball.union(al[v]) 1586 i += 1 1587 return ball 1588 1589 def smooth(self, niter=15, pass_band=0.1, edge_angle=15, feature_angle=60, boundary=False) -> Self: 1590 """ 1591 Adjust mesh point positions using the so-called "Windowed Sinc" method. 1592 1593 Arguments: 1594 niter : (int) 1595 number of iterations. 1596 pass_band : (float) 1597 set the pass_band value for the windowed sinc filter. 1598 edge_angle : (float) 1599 edge angle to control smoothing along edges (either interior or boundary). 1600 feature_angle : (float) 1601 specifies the feature angle for sharp edge identification. 1602 boundary : (bool) 1603 specify if boundary should also be smoothed or kept unmodified 1604 1605 Examples: