vedo.tetmesh
Work with tetrahedral meshes.
1#!/usr/bin/env python3 2# -*- coding: utf-8 -*- 3from weakref import ref as weak_ref_to 4 5import vedo.vtkclasses as vtk 6 7import numpy as np 8import vedo 9from vedo import utils 10from vedo.core import PointAlgorithms 11from vedo.mesh import Mesh 12from vedo.file_io import download 13from vedo.visual import MeshVisual 14from vedo.transformations import LinearTransform 15 16__docformat__ = "google" 17 18__doc__ = """ 19Work with tetrahedral meshes. 20 21 22""" 23 24__all__ = ["UnstructuredGrid", "TetMesh"] 25 26######################################################################### 27class UnstructuredGrid(PointAlgorithms, MeshVisual): 28 """Support for UnstructuredGrid objects.""" 29 30 def __init__(self, inputobj=None): 31 """ 32 Support for UnstructuredGrid objects. 33 34 Arguments: 35 inputobj : (list, vtkUnstructuredGrid, str) 36 A list in the form `[points, cells, celltypes]`, 37 or a vtkUnstructuredGrid object, or a filename 38 39 Celltypes are identified by the following 40 [convention](https://vtk.org/doc/nightly/html/vtkCellType_8h_source.html). 41 """ 42 super().__init__() 43 44 self.dataset = None 45 46 self.mapper = vtk.new("PolyDataMapper") 47 self._actor = vtk.vtkActor() 48 self._actor.retrieve_object = weak_ref_to(self) 49 self._actor.SetMapper(self.mapper) 50 self.properties = self._actor.GetProperty() 51 52 self.transform = LinearTransform() 53 54 self.name = "UnstructuredGrid" 55 self.filename = "" 56 self.info = {} 57 self.time = 0 58 self.rendered_at = set() 59 self._cmap_name = "" # remember the cmap name for self._keypress 60 61 ################### 62 inputtype = str(type(inputobj)) 63 64 if inputobj is None: 65 self.dataset = vtk.vtkUnstructuredGrid() 66 67 elif utils.is_sequence(inputobj): 68 69 pts, cells, celltypes = inputobj 70 assert len(cells) == len(celltypes) 71 72 self.dataset = vtk.vtkUnstructuredGrid() 73 74 if not utils.is_sequence(cells[0]): 75 tets = [] 76 nf = cells[0] + 1 77 for i, cl in enumerate(cells): 78 if i in (nf, 0): 79 k = i + 1 80 nf = cl + k 81 cell = [cells[j + k] for j in range(cl)] 82 tets.append(cell) 83 cells = tets 84 85 # This would fill the points and use those to define orientation 86 vpts = utils.numpy2vtk(pts, dtype=np.float32) 87 points = vtk.vtkPoints() 88 points.SetData(vpts) 89 self.dataset.SetPoints(points) 90 91 # Fill cells 92 # https://vtk.org/doc/nightly/html/vtkCellType_8h_source.html 93 for i, ct in enumerate(celltypes): 94 cell_conn = cells[i] 95 if ct == vtk.VTK_HEXAHEDRON: 96 cell = vtk.vtkHexahedron() 97 elif ct == vtk.VTK_TETRA: 98 cell = vtk.vtkTetra() 99 elif ct == vtk.VTK_VOXEL: 100 cell = vtk.vtkVoxel() 101 elif ct == vtk.VTK_WEDGE: 102 cell = vtk.vtkWedge() 103 elif ct == vtk.VTK_PYRAMID: 104 cell = vtk.vtkPyramid() 105 elif ct == vtk.VTK_HEXAGONAL_PRISM: 106 cell = vtk.vtkHexagonalPrism() 107 elif ct == vtk.VTK_PENTAGONAL_PRISM: 108 cell = vtk.vtkPentagonalPrism() 109 else: 110 print("UnstructuredGrid: cell type", ct, "not supported. Skip.") 111 continue 112 cpids = cell.GetPointIds() 113 for j, pid in enumerate(cell_conn): 114 cpids.SetId(j, pid) 115 self.dataset.InsertNextCell(ct, cpids) 116 117 elif "UnstructuredGrid" in inputtype: 118 self.dataset = inputobj 119 120 elif isinstance(inputobj, str): 121 if "https://" in inputobj: 122 inputobj = download(inputobj, verbose=False) 123 self.filename = inputobj 124 if inputobj.endswith(".vtu"): 125 reader = vtk.new("XMLUnstructuredGridReader") 126 else: 127 reader = vtk.new("UnstructuredGridReader") 128 self.filename = inputobj 129 reader.SetFileName(inputobj) 130 reader.Update() 131 self.dataset = reader.GetOutput() 132 133 else: 134 vedo.logger.error(f"cannot understand input type {inputtype}") 135 return 136 137 self.pipeline = utils.OperationNode( 138 self, comment=f"#cells {self.dataset.GetNumberOfCells()}", 139 c="#4cc9f0", 140 ) 141 142 # ------------------------------------------------------------------ 143 def __str__(self): 144 """Print a string summary of the `UnstructuredGrid` object.""" 145 module = self.__class__.__module__ 146 name = self.__class__.__name__ 147 out = vedo.printc( 148 f"{module}.{name} at ({hex(self.memory_address())})".ljust(75), 149 c="m", bold=True, invert=True, return_string=True, 150 ) 151 out += "\x1b[0m\u001b[35m" 152 153 out += "nr. of verts".ljust(14) + ": " + str(self.npoints) + "\n" 154 out += "nr. of cells".ljust(14) + ": " + str(self.ncells) + "\n" 155 156 if self.npoints: 157 out+="size".ljust(14)+ ": average=" + utils.precision(self.average_size(),6) 158 out+=", diagonal="+ utils.precision(self.diagonal_size(), 6)+ "\n" 159 out+="center of mass".ljust(14) + ": " + utils.precision(self.center_of_mass(),6)+"\n" 160 161 bnds = self.bounds() 162 bx1, bx2 = utils.precision(bnds[0], 3), utils.precision(bnds[1], 3) 163 by1, by2 = utils.precision(bnds[2], 3), utils.precision(bnds[3], 3) 164 bz1, bz2 = utils.precision(bnds[4], 3), utils.precision(bnds[5], 3) 165 out += "bounds".ljust(14) + ":" 166 out += " x=(" + bx1 + ", " + bx2 + ")," 167 out += " y=(" + by1 + ", " + by2 + ")," 168 out += " z=(" + bz1 + ", " + bz2 + ")\n" 169 170 for key in self.pointdata.keys(): 171 arr = self.pointdata[key] 172 rng = utils.precision(arr.min(), 3) + ", " + utils.precision(arr.max(), 3) 173 mark_active = "pointdata" 174 a_scalars = self.dataset.GetPointData().GetScalars() 175 a_vectors = self.dataset.GetPointData().GetVectors() 176 a_tensors = self.dataset.GetPointData().GetTensors() 177 if a_scalars and a_scalars.GetName() == key: 178 mark_active += " *" 179 elif a_vectors and a_vectors.GetName() == key: 180 mark_active += " **" 181 elif a_tensors and a_tensors.GetName() == key: 182 mark_active += " ***" 183 out += mark_active.ljust(14) + f': "{key}" ({arr.dtype}), ndim={arr.ndim}' 184 out += f", range=({rng})\n" 185 186 for key in self.celldata.keys(): 187 arr = self.celldata[key] 188 rng = utils.precision(arr.min(), 3) + ", " + utils.precision(arr.max(), 3) 189 mark_active = "celldata" 190 a_scalars = self.dataset.GetCellData().GetScalars() 191 a_vectors = self.dataset.GetCellData().GetVectors() 192 a_tensors = self.dataset.GetCellData().GetTensors() 193 if a_scalars and a_scalars.GetName() == key: 194 mark_active += " *" 195 elif a_vectors and a_vectors.GetName() == key: 196 mark_active += " **" 197 elif a_tensors and a_tensors.GetName() == key: 198 mark_active += " ***" 199 out += mark_active.ljust(14) + f': "{key}" ({arr.dtype}), ndim={arr.ndim}' 200 out += f", range=({rng})\n" 201 202 for key in self.metadata.keys(): 203 arr = self.metadata[key] 204 out+= "metadata".ljust(14) + ": " + f'"{key}" ({len(arr)} values)\n' 205 206 return out.rstrip() + "\x1b[0m" 207 208 209 def _repr_html_(self): 210 """ 211 HTML representation of the UnstructuredGrid object for Jupyter Notebooks. 212 213 Returns: 214 HTML text with the image and some properties. 215 """ 216 import io 217 import base64 218 from PIL import Image 219 220 library_name = "vedo.tetmesh.UnstructuredGrid" 221 help_url = "https://vedo.embl.es/docs/vedo/tetmesh.html" 222 223 # self.mapper.SetInputData(self.dataset) 224 arr = self.thumbnail() 225 im = Image.fromarray(arr) 226 buffered = io.BytesIO() 227 im.save(buffered, format="PNG", quality=100) 228 encoded = base64.b64encode(buffered.getvalue()).decode("utf-8") 229 url = "data:image/png;base64," + encoded 230 image = f"<img src='{url}'></img>" 231 232 bounds = "<br/>".join( 233 [ 234 utils.precision(min_x,4) + " ... " + utils.precision(max_x,4) 235 for min_x, max_x in zip(self.bounds()[::2], self.bounds()[1::2]) 236 ] 237 ) 238 239 help_text = "" 240 if self.name: 241 help_text += f"<b> {self.name}:   </b>" 242 help_text += '<b><a href="' + help_url + '" target="_blank">' + library_name + "</a></b>" 243 if self.filename: 244 dots = "" 245 if len(self.filename) > 30: 246 dots = "..." 247 help_text += f"<br/><code><i>({dots}{self.filename[-30:]})</i></code>" 248 249 pdata = "" 250 if self.dataset.GetPointData().GetScalars(): 251 if self.dataset.GetPointData().GetScalars().GetName(): 252 name = self.dataset.GetPointData().GetScalars().GetName() 253 pdata = "<tr><td><b> point data array </b></td><td>" + name + "</td></tr>" 254 255 cdata = "" 256 if self.dataset.GetCellData().GetScalars(): 257 if self.dataset.GetCellData().GetScalars().GetName(): 258 name = self.dataset.GetCellData().GetScalars().GetName() 259 cdata = "<tr><td><b> cell data array </b></td><td>" + name + "</td></tr>" 260 261 pts = self.vertices 262 cm = np.mean(pts, axis=0) 263 264 all = [ 265 "<table>", 266 "<tr>", 267 "<td>", image, "</td>", 268 "<td style='text-align: center; vertical-align: center;'><br/>", help_text, 269 "<table>", 270 "<tr><td><b> bounds </b> <br/> (x/y/z) </td><td>" + str(bounds) + "</td></tr>", 271 "<tr><td><b> center of mass </b></td><td>" + utils.precision(cm,3) + "</td></tr>", 272 # "<tr><td><b> average size </b></td><td>" + str(average_size) + "</td></tr>", 273 "<tr><td><b> nr. points / cells </b></td><td>" 274 + str(self.npoints) + " / " + str(self.ncells) + "</td></tr>", 275 pdata, 276 cdata, 277 "</table>", 278 "</table>", 279 ] 280 return "\n".join(all) 281 282 @property 283 def actor(self): 284 """Return the `vtkActor` of the object.""" 285 # print("building actor") 286 gf = vtk.new("GeometryFilter") 287 gf.SetInputData(self.dataset) 288 gf.Update() 289 out = gf.GetOutput() 290 self.mapper.SetInputData(out) 291 self.mapper.Modified() 292 return self._actor 293 294 @actor.setter 295 def actor(self, _): 296 pass 297 298 def _update(self, data, reset_locators=False): 299 self.dataset = data 300 # self.mapper.SetInputData(data) 301 # self.mapper.Modified() 302 ## self.actor.Modified() 303 return self 304 305 def copy(self, deep=True): 306 """Return a copy of the object. Alias of `clone()`.""" 307 return self.clone(deep=deep) 308 309 def clone(self, deep=True): 310 """Clone the UnstructuredGrid object to yield an exact copy.""" 311 ug = vtk.vtkUnstructuredGrid() 312 if deep: 313 ug.DeepCopy(self.dataset) 314 else: 315 ug.ShallowCopy(self.dataset) 316 if isinstance(self, vedo.UnstructuredGrid): 317 cloned = vedo.UnstructuredGrid(ug) 318 else: 319 cloned = vedo.TetMesh(ug) 320 321 cloned.copy_properties_from(self) 322 323 cloned.pipeline = utils.OperationNode( 324 "clone", parents=[self], shape='diamond', c='#bbe1ed', 325 ) 326 return cloned 327 328 def bounds(self): 329 """ 330 Get the object bounds. 331 Returns a list in format `[xmin,xmax, ymin,ymax, zmin,zmax]`. 332 """ 333 # OVERRIDE CommonAlgorithms.bounds() which is too slow 334 return self.dataset.GetBounds() 335 336 def threshold(self, name=None, above=None, below=None, on="cells"): 337 """ 338 Threshold the tetrahedral mesh by a cell scalar value. 339 Reduce to only tets which satisfy the threshold limits. 340 341 - if `above = below` will only select tets with that specific value. 342 - if `above > below` selection range is flipped. 343 344 Set keyword "on" to either "cells" or "points". 345 """ 346 th = vtk.new("Threshold") 347 th.SetInputData(self.dataset) 348 349 if name is None: 350 if self.celldata.keys(): 351 name = self.celldata.keys()[0] 352 th.SetInputArrayToProcess(0, 0, 0, 1, name) 353 elif self.pointdata.keys(): 354 name = self.pointdata.keys()[0] 355 th.SetInputArrayToProcess(0, 0, 0, 0, name) 356 if name is None: 357 vedo.logger.warning("cannot find active array. Skip.") 358 return self 359 else: 360 if on.startswith("c"): 361 th.SetInputArrayToProcess(0, 0, 0, 1, name) 362 else: 363 th.SetInputArrayToProcess(0, 0, 0, 0, name) 364 365 if above is not None: 366 th.SetLowerThreshold(above) 367 368 if below is not None: 369 th.SetUpperThreshold(below) 370 371 th.Update() 372 return self._update(th.GetOutput()) 373 374 def isosurface(self, value=None, flying_edges=True): 375 """ 376 Return an `Mesh` isosurface extracted from the `Volume` object. 377 378 Set `value` as single float or list of values to draw the isosurface(s). 379 Use flying_edges for faster results (but sometimes can interfere with `smooth()`). 380 381 Examples: 382 - [isosurfaces.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/isosurfaces.py) 383 384  385 """ 386 scrange = self.dataset.GetScalarRange() 387 388 if flying_edges: 389 cf = vtk.new("FlyingEdges3D") 390 cf.InterpolateAttributesOn() 391 else: 392 cf = vtk.new("ContourFilter") 393 cf.UseScalarTreeOn() 394 395 cf.SetInputData(self.dataset) 396 cf.ComputeNormalsOn() 397 398 if utils.is_sequence(value): 399 cf.SetNumberOfContours(len(value)) 400 for i, t in enumerate(value): 401 cf.SetValue(i, t) 402 else: 403 if value is None: 404 value = (2 * scrange[0] + scrange[1]) / 3.0 405 # print("automatic isosurface value =", value) 406 cf.SetValue(0, value) 407 408 cf.Update() 409 poly = cf.GetOutput() 410 411 out = vedo.mesh.Mesh(poly, c=None).phong() 412 out.mapper.SetScalarRange(scrange[0], scrange[1]) 413 414 out.pipeline = utils.OperationNode( 415 "isosurface", 416 parents=[self], 417 comment=f"#pts {out.dataset.GetNumberOfPoints()}", 418 c="#4cc9f0:#e9c46a", 419 ) 420 return out 421 422 def tomesh(self, fill=True, shrink=1.0): 423 """ 424 Build a polygonal `Mesh` from the current object. 425 426 If `fill=True`, the interior faces of all the cells are created. 427 (setting a `shrink` value slightly smaller than the default 1.0 428 can avoid flickering due to internal adjacent faces). 429 430 If `fill=False`, only the boundary faces will be generated. 431 """ 432 gf = vtk.new("GeometryFilter") 433 if fill: 434 sf = vtk.new("ShrinkFilter") 435 sf.SetInputData(self.dataset) 436 sf.SetShrinkFactor(shrink) 437 sf.Update() 438 gf.SetInputData(sf.GetOutput()) 439 gf.Update() 440 poly = gf.GetOutput() 441 else: 442 gf.SetInputData(self.dataset) 443 gf.Update() 444 poly = gf.GetOutput() 445 446 msh = vedo.mesh.Mesh(poly) 447 msh.copy_properties_from(self) 448 449 msh.pipeline = utils.OperationNode( 450 "tomesh", parents=[self], comment=f"fill={fill}", c="#9e2a2b:#e9c46a" 451 ) 452 return msh 453 454 def extract_cell_by_type(self, ctype): 455 """Extract a specific cell type and return a new `UnstructuredGrid`.""" 456 uarr = self.dataset.GetCellTypesArray() 457 ctarrtyp = np.where(utils.vtk2numpy(uarr) == ctype)[0] 458 uarrtyp = utils.numpy2vtk(ctarrtyp, deep=False, dtype="id") 459 selection_node = vtk.new("SelectionNode") 460 selection_node.SetFieldType(vtk.get_class("SelectionNode").CELL) 461 selection_node.SetContentType(vtk.get_class("SelectionNode").INDICES) 462 selection_node.SetSelectionList(uarrtyp) 463 selection = vtk.new("Selection") 464 selection.AddNode(selection_node) 465 es = vtk.new("ExtractSelection") 466 es.SetInputData(0, self.dataset) 467 es.SetInputData(1, selection) 468 es.Update() 469 470 ug = UnstructuredGrid(es.GetOutput()) 471 472 ug.pipeline = utils.OperationNode( 473 "extract_cell_type", comment=f"type {ctype}", 474 c="#edabab", parents=[self], 475 ) 476 return ug 477 478 def extract_cells_by_id(self, idlist, use_point_ids=False): 479 """Return a new `UnstructuredGrid` composed of the specified subset of indices.""" 480 selection_node = vtk.new("SelectionNode") 481 if use_point_ids: 482 selection_node.SetFieldType(vtk.get_class("SelectionNode").POINT) 483 contcells = vtk.get_class("SelectionNode").CONTAINING_CELLS() 484 selection_node.GetProperties().Set(contcells, 1) 485 else: 486 selection_node.SetFieldType(vtk.get_class("SelectionNode").CELL) 487 selection_node.SetContentType(vtk.get_class("SelectionNode").INDICES) 488 vidlist = utils.numpy2vtk(idlist, dtype="id") 489 selection_node.SetSelectionList(vidlist) 490 selection = vtk.new("Selection") 491 selection.AddNode(selection_node) 492 es = vtk.new("ExtractSelection") 493 es.SetInputData(0, self) 494 es.SetInputData(1, selection) 495 es.Update() 496 497 ug = vedo.tetmesh.UnstructuredGrid(es.GetOutput()) 498 pr = vtk.vtkProperty() 499 pr.DeepCopy(self.properties) 500 ug.SetProperty(pr) 501 ug.properties = pr 502 503 ug.mapper.SetLookupTable(utils.ctf2lut(self)) 504 ug.pipeline = utils.OperationNode( 505 "extract_cells_by_id", 506 parents=[self], 507 comment=f"#cells {self.dataset.GetNumberOfCells()}", 508 c="#9e2a2b", 509 ) 510 return ug 511 512 def find_cell(self, p): 513 """Locate the cell that contains a point and return the cell ID.""" 514 cell = vtk.vtkTetra() 515 cell_id = vtk.mutable(0) 516 tol2 = vtk.mutable(0) 517 sub_id = vtk.mutable(0) 518 pcoords = [0, 0, 0] 519 weights = [0, 0, 0] 520 cid = self.dataset.FindCell( 521 p, cell, cell_id, tol2, sub_id, pcoords, weights) 522 return cid 523 524 def clean(self): 525 """ 526 Cleanup unused points and empty cells 527 """ 528 cl = vtk.new("StaticCleanUnstructuredGrid") 529 cl.SetInputData(self.dataset) 530 cl.RemoveUnusedPointsOn() 531 cl.ProduceMergeMapOff() 532 cl.AveragePointDataOff() 533 cl.Update() 534 535 self._update(cl.GetOutput()) 536 self.pipeline = utils.OperationNode( 537 "clean", 538 parents=[self], 539 comment=f"#cells {self.dataset.GetNumberOfCells()}", 540 c="#9e2a2b", 541 ) 542 return self 543 544 def extract_cells_on_plane(self, origin, normal): 545 """ 546 Extract cells that are lying of the specified surface. 547 """ 548 bf = vtk.new("3DLinearGridCrinkleExtractor") 549 bf.SetInputData(self.dataset) 550 bf.CopyPointDataOn() 551 bf.CopyCellDataOn() 552 bf.RemoveUnusedPointsOff() 553 554 plane = vtk.new("Plane") 555 plane.SetOrigin(origin) 556 plane.SetNormal(normal) 557 bf.SetImplicitFunction(plane) 558 bf.Update() 559 560 self._update(bf.GetOutput(), reset_locators=False) 561 self.pipeline = utils.OperationNode( 562 "extract_cells_on_plane", 563 parents=[self], 564 comment=f"#cells {self.dataset.GetNumberOfCells()}", 565 c="#9e2a2b", 566 ) 567 return self 568 569 def extract_cells_on_sphere(self, center, radius): 570 """ 571 Extract cells that are lying of the specified surface. 572 """ 573 bf = vtk.new("3DLinearGridCrinkleExtractor") 574 bf.SetInputData(self.dataset) 575 bf.CopyPointDataOn() 576 bf.CopyCellDataOn() 577 bf.RemoveUnusedPointsOff() 578 579 sph = vtk.new("Sphere") 580 sph.SetRadius(radius) 581 sph.SetCenter(center) 582 bf.SetImplicitFunction(sph) 583 bf.Update() 584 585 self._update(bf.GetOutput()) 586 self.pipeline = utils.OperationNode( 587 "extract_cells_on_sphere", 588 parents=[self], 589 comment=f"#cells {self.dataset.GetNumberOfCells()}", 590 c="#9e2a2b", 591 ) 592 return self 593 594 def extract_cells_on_cylinder(self, center, axis, radius): 595 """ 596 Extract cells that are lying of the specified surface. 597 """ 598 bf = vtk.new("3DLinearGridCrinkleExtractor") 599 bf.SetInputData(self.dataset) 600 bf.CopyPointDataOn() 601 bf.CopyCellDataOn() 602 bf.RemoveUnusedPointsOff() 603 604 cyl = vtk.new("Cylinder") 605 cyl.SetRadius(radius) 606 cyl.SetCenter(center) 607 cyl.SetAxis(axis) 608 bf.SetImplicitFunction(cyl) 609 bf.Update() 610 611 self.pipeline = utils.OperationNode( 612 "extract_cells_on_cylinder", 613 parents=[self], 614 comment=f"#cells {self.dataset.GetNumberOfCells()}", 615 c="#9e2a2b", 616 ) 617 self._update(bf.GetOutput()) 618 return self 619 620 def cut_with_plane(self, origin=(0, 0, 0), normal="x"): 621 """ 622 Cut the object with the plane defined by a point and a normal. 623 624 Arguments: 625 origin : (list) 626 the cutting plane goes through this point 627 normal : (list, str) 628 normal vector to the cutting plane 629 """ 630 # if isinstance(self, vedo.Volume): 631 # raise RuntimeError("cut_with_plane() is not applicable to Volume objects.") 632 633 strn = str(normal) 634 if strn == "x": normal = (1, 0, 0) 635 elif strn == "y": normal = (0, 1, 0) 636 elif strn == "z": normal = (0, 0, 1) 637 elif strn == "-x": normal = (-1, 0, 0) 638 elif strn == "-y": normal = (0, -1, 0) 639 elif strn == "-z": normal = (0, 0, -1) 640 plane = vtk.new("Plane") 641 plane.SetOrigin(origin) 642 plane.SetNormal(normal) 643 clipper = vtk.new("ClipDataSet") 644 clipper.SetInputData(self.dataset) 645 clipper.SetClipFunction(plane) 646 clipper.GenerateClipScalarsOff() 647 clipper.GenerateClippedOutputOff() 648 clipper.SetValue(0) 649 clipper.Update() 650 cout = clipper.GetOutput() 651 652 if isinstance(cout, vtk.vtkUnstructuredGrid): 653 ug = vedo.UnstructuredGrid(cout) 654 if isinstance(self, vedo.UnstructuredGrid): 655 self._update(cout) 656 self.pipeline = utils.OperationNode("cut_with_plane", parents=[self], c="#9e2a2b") 657 return self 658 ug.pipeline = utils.OperationNode("cut_with_plane", parents=[self], c="#9e2a2b") 659 return ug 660 661 else: 662 self._update(cout) 663 self.pipeline = utils.OperationNode("cut_with_plane", parents=[self], c="#9e2a2b") 664 return self 665 666 def cut_with_box(self, box): 667 """ 668 Cut the grid with the specified bounding box. 669 670 Parameter box has format [xmin, xmax, ymin, ymax, zmin, zmax]. 671 If an object is passed, its bounding box are used. 672 673 This method always returns a TetMesh object. 674 675 Example: 676 ```python 677 from vedo import * 678 tetmesh = TetMesh(dataurl+'limb_ugrid.vtk') 679 tetmesh.color('rainbow') 680 cu = Cube(side=500).x(500) # any Mesh works 681 tetmesh.cut_with_box(cu).show(axes=1) 682 ``` 683 684  685 """ 686 bc = vtk.new("BoxClipDataSet") 687 bc.SetInputData(self.dataset) 688 try: 689 boxb = box.bounds() 690 except AttributeError: 691 boxb = box 692 693 bc.SetBoxClip(*boxb) 694 bc.Update() 695 cout = bc.GetOutput() 696 697 # output of vtkBoxClipDataSet is always tetrahedrons 698 tm = vedo.TetMesh(cout) 699 tm.pipeline = utils.OperationNode("cut_with_box", parents=[self], c="#9e2a2b") 700 return tm 701 702 703 def cut_with_mesh( 704 self, mesh, invert=False, whole_cells=False, on_boundary=False 705 ): 706 """ 707 Cut a UnstructuredGrid or TetMesh with a Mesh. 708 709 Use `invert` to return cut off part of the input object. 710 """ 711 ug = self.dataset 712 713 ippd = vtk.new("ImplicitPolyDataDistance") 714 ippd.SetInput(mesh.dataset) 715 716 if whole_cells or on_boundary: 717 clipper = vtk.new("ExtractGeometry") 718 clipper.SetInputData(ug) 719 clipper.SetImplicitFunction(ippd) 720 clipper.SetExtractInside(not invert) 721 clipper.SetExtractBoundaryCells(False) 722 if on_boundary: 723 clipper.SetExtractBoundaryCells(True) 724 clipper.SetExtractOnlyBoundaryCells(True) 725 else: 726 signed_dists = vtk.vtkFloatArray() 727 signed_dists.SetNumberOfComponents(1) 728 signed_dists.SetName("SignedDistance") 729 for pointId in range(ug.GetNumberOfPoints()): 730 p = ug.GetPoint(pointId) 731 signed_dist = ippd.EvaluateFunction(p) 732 signed_dists.InsertNextValue(signed_dist) 733 ug.GetPointData().AddArray(signed_dists) 734 ug.GetPointData().SetActiveScalars("SignedDistance") # NEEDED 735 clipper = vtk.new("ClipDataSet") 736 clipper.SetInputData(ug) 737 clipper.SetInsideOut(not invert) 738 clipper.SetValue(0.0) 739 740 clipper.Update() 741 742 out = vedo.UnstructuredGrid(clipper.GetOutput()) 743 out.pipeline = utils.OperationNode("cut_with_mesh", parents=[self], c="#9e2a2b") 744 return out 745 746########################################################################## 747class TetMesh(UnstructuredGrid): 748 """The class describing tetrahedral meshes.""" 749 750 def __init__(self, inputobj=None): 751 """ 752 Arguments: 753 inputobj : (vtkDataSet, list, str) 754 list of points and tet indices, or filename 755 """ 756 super().__init__() 757 758 self.dataset = None 759 760 self.mapper = vtk.new("PolyDataMapper") 761 self._actor = vtk.vtkActor() 762 self._actor.retrieve_object = weak_ref_to(self) 763 self._actor.SetMapper(self.mapper) 764 self.properties = self._actor.GetProperty() 765 766 self.transform = LinearTransform() 767 768 self.name = "TetMesh" 769 self.filename = "" 770 771 772 # inputtype = str(type(inputobj)) 773 # print('TetMesh inputtype', inputtype) 774 775 ################### 776 if inputobj is None: 777 self.dataset = vtk.vtkUnstructuredGrid() 778 779 elif isinstance(inputobj, vtk.vtkUnstructuredGrid): 780 self.dataset = inputobj 781 782 elif isinstance(inputobj, UnstructuredGrid): 783 self.dataset = inputobj.dataset 784 785 elif isinstance(inputobj, vtk.vtkRectilinearGrid): 786 r2t = vtk.new("RectilinearGridToTetrahedra") 787 r2t.SetInputData(inputobj) 788 r2t.RememberVoxelIdOn() 789 r2t.SetTetraPerCellTo6() 790 r2t.Update() 791 self.dataset = r2t.GetOutput() 792 793 elif isinstance(inputobj, vtk.vtkDataSet): 794 r2t = vtk.new("DataSetTriangleFilter") 795 r2t.SetInputData(inputobj) 796 r2t.TetrahedraOnlyOn() 797 r2t.Update() 798 self.dataset = r2t.GetOutput() 799 800 elif isinstance(inputobj, str): 801 if "https://" in inputobj: 802 inputobj = download(inputobj, verbose=False) 803 if inputobj.endswith(".vtu"): 804 reader = vtk.new("XMLUnstructuredGridReader") 805 else: 806 reader = vtk.new("UnstructuredGridReader") 807 self.filename = inputobj 808 reader.SetFileName(inputobj) 809 reader.Update() 810 ug = reader.GetOutput() 811 tt = vtk.new("DataSetTriangleFilter") 812 tt.SetInputData(ug) 813 tt.SetTetrahedraOnly(True) 814 tt.Update() 815 self.dataset = tt.GetOutput() 816 817 elif utils.is_sequence(inputobj): 818 self.dataset = vtk.vtkUnstructuredGrid() 819 820 points, cells = inputobj 821 if len(points) == 0: 822 return 823 if not utils.is_sequence(points[0]): 824 return 825 if len(cells) == 0: 826 return 827 828 if not utils.is_sequence(cells[0]): 829 tets = [] 830 nf = cells[0] + 1 831 for i, cl in enumerate(cells): 832 if i in (nf, 0): 833 k = i + 1 834 nf = cl + k 835 cell = [cells[j + k] for j in range(cl)] 836 tets.append(cell) 837 cells = tets 838 839 source_points = vtk.vtkPoints() 840 varr = utils.numpy2vtk(points, dtype=np.float32) 841 source_points.SetData(varr) 842 self.dataset.SetPoints(source_points) 843 844 source_tets = vtk.vtkCellArray() 845 for f in cells: 846 ele = vtk.vtkTetra() 847 pid = ele.GetPointIds() 848 for i, fi in enumerate(f): 849 pid.SetId(i, fi) 850 source_tets.InsertNextCell(ele) 851 self.dataset.SetCells(vtk.VTK_TETRA, source_tets) 852 853 else: 854 vedo.logger.error(f"cannot understand input type {inputtype}") 855 return 856 857 self.pipeline = utils.OperationNode( 858 self, comment=f"#tets {self.dataset.GetNumberOfCells()}", 859 c="#9e2a2b", 860 ) 861 862 ################################################################## 863 def __str__(self): 864 """Print a string summary of the `TetMesh` object.""" 865 module = self.__class__.__module__ 866 name = self.__class__.__name__ 867 out = vedo.printc( 868 f"{module}.{name} at ({hex(self.memory_address())})".ljust(75), 869 c="c", bold=True, invert=True, return_string=True, 870 ) 871 out += "\x1b[0m\u001b[36m" 872 873 out += "nr. of verts".ljust(14) + ": " + str(self.npoints) + "\n" 874 out += "nr. of tetras".ljust(14)+ ": " + str(self.ncells) + "\n" 875 876 if self.npoints: 877 out+="size".ljust(14)+ ": average=" + utils.precision(self.average_size(),6) 878 out+=", diagonal="+ utils.precision(self.diagonal_size(), 6)+ "\n" 879 out+="center of mass".ljust(14) + ": " + utils.precision(self.center_of_mass(),6)+"\n" 880 881 bnds = self.bounds() 882 bx1, bx2 = utils.precision(bnds[0], 3), utils.precision(bnds[1], 3) 883 by1, by2 = utils.precision(bnds[2], 3), utils.precision(bnds[3], 3) 884 bz1, bz2 = utils.precision(bnds[4], 3), utils.precision(bnds[5], 3) 885 out += "bounds".ljust(14) + ":" 886 out += " x=(" + bx1 + ", " + bx2 + ")," 887 out += " y=(" + by1 + ", " + by2 + ")," 888 out += " z=(" + bz1 + ", " + bz2 + ")\n" 889 890 for key in self.pointdata.keys(): 891 arr = self.pointdata[key] 892 rng = utils.precision(arr.min(), 3) + ", " + utils.precision(arr.max(), 3) 893 mark_active = "pointdata" 894 a_scalars = self.dataset.GetPointData().GetScalars() 895 a_vectors = self.dataset.GetPointData().GetVectors() 896 a_tensors = self.dataset.GetPointData().GetTensors() 897 if a_scalars and a_scalars.GetName() == key: 898 mark_active += " *" 899 elif a_vectors and a_vectors.GetName() == key: 900 mark_active += " **" 901 elif a_tensors and a_tensors.GetName() == key: 902 mark_active += " ***" 903 out += mark_active.ljust(14) + f': "{key}" ({arr.dtype}), ndim={arr.ndim}' 904 out += f", range=({rng})\n" 905 906 for key in self.celldata.keys(): 907 arr = self.celldata[key] 908 rng = utils.precision(arr.min(), 3) + ", " + utils.precision(arr.max(), 3) 909 mark_active = "celldata" 910 a_scalars = self.dataset.GetCellData().GetScalars() 911 a_vectors = self.dataset.GetCellData().GetVectors() 912 a_tensors = self.dataset.GetCellData().GetTensors() 913 if a_scalars and a_scalars.GetName() == key: 914 mark_active += " *" 915 elif a_vectors and a_vectors.GetName() == key: 916 mark_active += " **" 917 elif a_tensors and a_tensors.GetName() == key: 918 mark_active += " ***" 919 out += mark_active.ljust(14) + f': "{key}" ({arr.dtype}), ndim={arr.ndim}' 920 out += f", range=({rng})\n" 921 922 for key in self.metadata.keys(): 923 arr = self.metadata[key] 924 out+= "metadata".ljust(14) + ": " + f'"{key}" ({len(arr)} values)\n' 925 926 return out.rstrip() + "\x1b[0m" 927 928 929 def _repr_html_(self): 930 """ 931 HTML representation of the TetMesh object for Jupyter Notebooks. 932 933 Returns: 934 HTML text with the image and some properties. 935 """ 936 import io 937 import base64 938 from PIL import Image 939 940 library_name = "vedo.tetmesh.TetMesh" 941 help_url = "https://vedo.embl.es/docs/vedo/tetmesh.html" 942 943 arr = self.thumbnail() 944 im = Image.fromarray(arr) 945 buffered = io.BytesIO() 946 im.save(buffered, format="PNG", quality=100) 947 encoded = base64.b64encode(buffered.getvalue()).decode("utf-8") 948 url = "data:image/png;base64," + encoded 949 image = f"<img src='{url}'></img>" 950 951 bounds = "<br/>".join( 952 [ 953 utils.precision(min_x,4) + " ... " + utils.precision(max_x,4) 954 for min_x, max_x in zip(self.bounds()[::2], self.bounds()[1::2]) 955 ] 956 ) 957 958 help_text = "" 959 if self.name: 960 help_text += f"<b> {self.name}:   </b>" 961 help_text += '<b><a href="' + help_url + '" target="_blank">' + library_name + "</a></b>" 962 if self.filename: 963 dots = "" 964 if len(self.filename) > 30: 965 dots = "..." 966 help_text += f"<br/><code><i>({dots}{self.filename[-30:]})</i></code>" 967 968 pdata = "" 969 if self.dataset.GetPointData().GetScalars(): 970 if self.dataset.GetPointData().GetScalars().GetName(): 971 name = self.dataset.GetPointData().GetScalars().GetName() 972 pdata = "<tr><td><b> point data array </b></td><td>" + name + "</td></tr>" 973 974 cdata = "" 975 if self.dataset.GetCellData().GetScalars(): 976 if self.dataset.GetCellData().GetScalars().GetName(): 977 name = self.dataset.GetCellData().GetScalars().GetName() 978 cdata = "<tr><td><b> cell data array </b></td><td>" + name + "</td></tr>" 979 980 pts = self.vertices 981 cm = np.mean(pts, axis=0) 982 983 allt = [ 984 "<table>", 985 "<tr>", 986 "<td>", image, "</td>", 987 "<td style='text-align: center; vertical-align: center;'><br/>", help_text, 988 "<table>", 989 "<tr><td><b> bounds </b> <br/> (x/y/z) </td><td>" + str(bounds) + "</td></tr>", 990 "<tr><td><b> center of mass </b></td><td>" + utils.precision(cm,3) + "</td></tr>", 991 "<tr><td><b> nr. points / tets </b></td><td>" 992 + str(self.npoints) + " / " + str(self.ncells) + "</td></tr>", 993 pdata, 994 cdata, 995 "</table>", 996 "</table>", 997 ] 998 return "\n".join(allt) 999 1000 def compute_quality(self, metric=7): 1001 """ 1002 Calculate functions of quality for the elements of a triangular mesh. 1003 This method adds to the mesh a cell array named "Quality". 1004 1005 Arguments: 1006 metric : (int) 1007 type of estimators: 1008 - EDGE RATIO, 0 1009 - ASPECT RATIO, 1 1010 - RADIUS RATIO, 2 1011 - ASPECT FROBENIUS, 3 1012 - MIN_ANGLE, 4 1013 - COLLAPSE RATIO, 5 1014 - ASPECT GAMMA, 6 1015 - VOLUME, 7 1016 - ... 1017 1018 See class [vtkMeshQuality](https://vtk.org/doc/nightly/html/classvtkMeshQuality.html) 1019 for an explanation of the meaning of each metric.. 1020 """ 1021 qf = vtk.new("MeshQuality") 1022 qf.SetInputData(self.dataset) 1023 qf.SetTetQualityMeasure(metric) 1024 qf.SaveCellQualityOn() 1025 qf.Update() 1026 self._update(qf.GetOutput()) 1027 return utils.vtk2numpy(qf.GetOutput().GetCellData().GetArray("Quality")) 1028 1029 def check_validity(self, tol=0): 1030 """ 1031 Return an array of possible problematic tets following this convention: 1032 ```python 1033 Valid = 0 1034 WrongNumberOfPoints = 01 1035 IntersectingEdges = 02 1036 IntersectingFaces = 04 1037 NoncontiguousEdges = 08 1038 Nonconvex = 10 1039 OrientedIncorrectly = 20 1040 ``` 1041 1042 Arguments: 1043 tol : (float) 1044 This value is used as an epsilon for floating point 1045 equality checks throughout the cell checking process. 1046 """ 1047 vald = vtk.new("CellValidator") 1048 if tol: 1049 vald.SetTolerance(tol) 1050 vald.SetInputData(self.dataset) 1051 vald.Update() 1052 varr = vald.GetOutput().GetCellData().GetArray("ValidityState") 1053 return utils.vtk2numpy(varr) 1054 1055 def decimate(self, scalars_name, fraction=0.5, n=0): 1056 """ 1057 Downsample the number of tets in a TetMesh to a specified fraction. 1058 Either `fraction` or `n` must be set. 1059 1060 Arguments: 1061 fraction : (float) 1062 the desired final fraction of the total. 1063 n : (int) 1064 the desired number of final tets 1065 1066 .. note:: setting `fraction=0.1` leaves 10% of the original nr of tets. 1067 """ 1068 decimate = vtk.new("UnstructuredGridQuadricDecimation") 1069 decimate.SetInputData(self.dataset) 1070 decimate.SetScalarsName(scalars_name) 1071 1072 if n: # n = desired number of points 1073 decimate.SetNumberOfTetsOutput(n) 1074 else: 1075 decimate.SetTargetReduction(1 - fraction) 1076 decimate.Update() 1077 self._update(decimate.GetOutput()) 1078 self.pipeline = utils.OperationNode( 1079 "decimate", comment=f"array: {scalars_name}", 1080 c="#edabab", parents=[self], 1081 ) 1082 return self 1083 1084 def subdvide(self): 1085 """ 1086 Increase the number of tets of a `TetMesh`. 1087 Subdivide one tetrahedron into twelve for every tetra. 1088 """ 1089 sd = vtk.new("SubdivideTetra") 1090 sd.SetInputData(self.dataset) 1091 sd.Update() 1092 self._update(sd.GetOutput()) 1093 self.pipeline = utils.OperationNode( 1094 "subdvide", c="#edabab", parents=[self], 1095 ) 1096 return self 1097 1098 def isosurface(self, value=True): 1099 """ 1100 Return a `vedo.Mesh` isosurface. 1101 1102 Set `value` to a single value or list of values to compute the isosurface(s). 1103 """ 1104 if not self.dataset.GetPointData().GetScalars(): 1105 self.map_cells_to_points() 1106 scrange = self.dataset.GetPointData().GetScalars().GetRange() 1107 cf = vtk.new("ContourFilter") # vtk.new("ContourGrid") 1108 cf.SetInputData(self.dataset) 1109 1110 if utils.is_sequence(value): 1111 cf.SetNumberOfContours(len(value)) 1112 for i, t in enumerate(value): 1113 cf.SetValue(i, t) 1114 cf.Update() 1115 else: 1116 if value is True: 1117 value = (2 * scrange[0] + scrange[1]) / 3.0 1118 cf.SetValue(0, value) 1119 cf.Update() 1120 1121 clp = vtk.new("CleanPolyData") 1122 clp.SetInputData(cf.GetOutput()) 1123 clp.Update() 1124 msh = Mesh(clp.GetOutput(), c=None).phong() 1125 msh.copy_properties_from(self) 1126 msh.pipeline = utils.OperationNode("isosurface", c="#edabab", parents=[self]) 1127 return msh 1128 1129 def slice(self, origin=(0, 0, 0), normal=(1, 0, 0)): 1130 """ 1131 Return a 2D slice of the mesh by a plane passing through origin and assigned normal. 1132 """ 1133 strn = str(normal) 1134 if strn == "x": normal = (1, 0, 0) 1135 elif strn == "y": normal = (0, 1, 0) 1136 elif strn == "z": normal = (0, 0, 1) 1137 elif strn == "-x": normal = (-1, 0, 0) 1138 elif strn == "-y": normal = (0, -1, 0) 1139 elif strn == "-z": normal = (0, 0, -1) 1140 plane = vtk.new("Plane") 1141 plane.SetOrigin(origin) 1142 plane.SetNormal(normal) 1143 1144 cc = vtk.new("Cutter") 1145 cc.SetInputData(self.dataset) 1146 cc.SetCutFunction(plane) 1147 cc.Update() 1148 msh = Mesh(cc.GetOutput()).flat().lighting("ambient") 1149 msh.copy_properties_from(self) 1150 msh.pipeline = utils.OperationNode("slice", c="#edabab", parents=[self]) 1151 return msh
28class UnstructuredGrid(PointAlgorithms, MeshVisual): 29 """Support for UnstructuredGrid objects.""" 30 31 def __init__(self, inputobj=None): 32 """ 33 Support for UnstructuredGrid objects. 34 35 Arguments: 36 inputobj : (list, vtkUnstructuredGrid, str) 37 A list in the form `[points, cells, celltypes]`, 38 or a vtkUnstructuredGrid object, or a filename 39 40 Celltypes are identified by the following 41 [convention](https://vtk.org/doc/nightly/html/vtkCellType_8h_source.html). 42 """ 43 super().__init__() 44 45 self.dataset = None 46 47 self.mapper = vtk.new("PolyDataMapper") 48 self._actor = vtk.vtkActor() 49 self._actor.retrieve_object = weak_ref_to(self) 50 self._actor.SetMapper(self.mapper) 51 self.properties = self._actor.GetProperty() 52 53 self.transform = LinearTransform() 54 55 self.name = "UnstructuredGrid" 56 self.filename = "" 57 self.info = {} 58 self.time = 0 59 self.rendered_at = set() 60 self._cmap_name = "" # remember the cmap name for self._keypress 61 62 ################### 63 inputtype = str(type(inputobj)) 64 65 if inputobj is None: 66 self.dataset = vtk.vtkUnstructuredGrid() 67 68 elif utils.is_sequence(inputobj): 69 70 pts, cells, celltypes = inputobj 71 assert len(cells) == len(celltypes) 72 73 self.dataset = vtk.vtkUnstructuredGrid() 74 75 if not utils.is_sequence(cells[0]): 76 tets = [] 77 nf = cells[0] + 1 78 for i, cl in enumerate(cells): 79 if i in (nf, 0): 80 k = i + 1 81 nf = cl + k 82 cell = [cells[j + k] for j in range(cl)] 83 tets.append(cell) 84 cells = tets 85 86 # This would fill the points and use those to define orientation 87 vpts = utils.numpy2vtk(pts, dtype=np.float32) 88 points = vtk.vtkPoints() 89 points.SetData(vpts) 90 self.dataset.SetPoints(points) 91 92 # Fill cells 93 # https://vtk.org/doc/nightly/html/vtkCellType_8h_source.html 94 for i, ct in enumerate(celltypes): 95 cell_conn = cells[i] 96 if ct == vtk.VTK_HEXAHEDRON: 97 cell = vtk.vtkHexahedron() 98 elif ct == vtk.VTK_TETRA: 99 cell = vtk.vtkTetra() 100 elif ct == vtk.VTK_VOXEL: 101 cell = vtk.vtkVoxel() 102 elif ct == vtk.VTK_WEDGE: 103 cell = vtk.vtkWedge() 104 elif ct == vtk.VTK_PYRAMID: 105 cell = vtk.vtkPyramid() 106 elif ct == vtk.VTK_HEXAGONAL_PRISM: 107 cell = vtk.vtkHexagonalPrism() 108 elif ct == vtk.VTK_PENTAGONAL_PRISM: 109 cell = vtk.vtkPentagonalPrism() 110 else: 111 print("UnstructuredGrid: cell type", ct, "not supported. Skip.") 112 continue 113 cpids = cell.GetPointIds() 114 for j, pid in enumerate(cell_conn): 115 cpids.SetId(j, pid) 116 self.dataset.InsertNextCell(ct, cpids) 117 118 elif "UnstructuredGrid" in inputtype: 119 self.dataset = inputobj 120 121 elif isinstance(inputobj, str): 122 if "https://" in inputobj: 123 inputobj = download(inputobj, verbose=False) 124 self.filename = inputobj 125 if inputobj.endswith(".vtu"): 126 reader = vtk.new("XMLUnstructuredGridReader") 127 else: 128 reader = vtk.new("UnstructuredGridReader") 129 self.filename = inputobj 130 reader.SetFileName(inputobj) 131 reader.Update() 132 self.dataset = reader.GetOutput() 133 134 else: 135 vedo.logger.error(f"cannot understand input type {inputtype}") 136 return 137 138 self.pipeline = utils.OperationNode( 139 self, comment=f"#cells {self.dataset.GetNumberOfCells()}", 140 c="#4cc9f0", 141 ) 142 143 # ------------------------------------------------------------------ 144 def __str__(self): 145 """Print a string summary of the `UnstructuredGrid` object.""" 146 module = self.__class__.__module__ 147 name = self.__class__.__name__ 148 out = vedo.printc( 149 f"{module}.{name} at ({hex(self.memory_address())})".ljust(75), 150 c="m", bold=True, invert=True, return_string=True, 151 ) 152 out += "\x1b[0m\u001b[35m" 153 154 out += "nr. of verts".ljust(14) + ": " + str(self.npoints) + "\n" 155 out += "nr. of cells".ljust(14) + ": " + str(self.ncells) + "\n" 156 157 if self.npoints: 158 out+="size".ljust(14)+ ": average=" + utils.precision(self.average_size(),6) 159 out+=", diagonal="+ utils.precision(self.diagonal_size(), 6)+ "\n" 160 out+="center of mass".ljust(14) + ": " + utils.precision(self.center_of_mass(),6)+"\n" 161 162 bnds = self.bounds() 163 bx1, bx2 = utils.precision(bnds[0], 3), utils.precision(bnds[1], 3) 164 by1, by2 = utils.precision(bnds[2], 3), utils.precision(bnds[3], 3) 165 bz1, bz2 = utils.precision(bnds[4], 3), utils.precision(bnds[5], 3) 166 out += "bounds".ljust(14) + ":" 167 out += " x=(" + bx1 + ", " + bx2 + ")," 168 out += " y=(" + by1 + ", " + by2 + ")," 169 out += " z=(" + bz1 + ", " + bz2 + ")\n" 170 171 for key in self.pointdata.keys(): 172 arr = self.pointdata[key] 173 rng = utils.precision(arr.min(), 3) + ", " + utils.precision(arr.max(), 3) 174 mark_active = "pointdata" 175 a_scalars = self.dataset.GetPointData().GetScalars() 176 a_vectors = self.dataset.GetPointData().GetVectors() 177 a_tensors = self.dataset.GetPointData().GetTensors() 178 if a_scalars and a_scalars.GetName() == key: 179 mark_active += " *" 180 elif a_vectors and a_vectors.GetName() == key: 181 mark_active += " **" 182 elif a_tensors and a_tensors.GetName() == key: 183 mark_active += " ***" 184 out += mark_active.ljust(14) + f': "{key}" ({arr.dtype}), ndim={arr.ndim}' 185 out += f", range=({rng})\n" 186 187 for key in self.celldata.keys(): 188 arr = self.celldata[key] 189 rng = utils.precision(arr.min(), 3) + ", " + utils.precision(arr.max(), 3) 190 mark_active = "celldata" 191 a_scalars = self.dataset.GetCellData().GetScalars() 192 a_vectors = self.dataset.GetCellData().GetVectors() 193 a_tensors = self.dataset.GetCellData().GetTensors() 194 if a_scalars and a_scalars.GetName() == key: 195 mark_active += " *" 196 elif a_vectors and a_vectors.GetName() == key: 197 mark_active += " **" 198 elif a_tensors and a_tensors.GetName() == key: 199 mark_active += " ***" 200 out += mark_active.ljust(14) + f': "{key}" ({arr.dtype}), ndim={arr.ndim}' 201 out += f", range=({rng})\n" 202 203 for key in self.metadata.keys(): 204 arr = self.metadata[key] 205 out+= "metadata".ljust(14) + ": " + f'"{key}" ({len(arr)} values)\n' 206 207 return out.rstrip() + "\x1b[0m" 208 209 210 def _repr_html_(self): 211 """ 212 HTML representation of the UnstructuredGrid object for Jupyter Notebooks. 213 214 Returns: 215 HTML text with the image and some properties. 216 """ 217 import io 218 import base64 219 from PIL import Image 220 221 library_name = "vedo.tetmesh.UnstructuredGrid" 222 help_url = "https://vedo.embl.es/docs/vedo/tetmesh.html" 223 224 # self.mapper.SetInputData(self.dataset) 225 arr = self.thumbnail() 226 im = Image.fromarray(arr) 227 buffered = io.BytesIO() 228 im.save(buffered, format="PNG", quality=100) 229 encoded = base64.b64encode(buffered.getvalue()).decode("utf-8") 230 url = "data:image/png;base64," + encoded 231 image = f"<img src='{url}'></img>" 232 233 bounds = "<br/>".join( 234 [ 235 utils.precision(min_x,4) + " ... " + utils.precision(max_x,4) 236 for min_x, max_x in zip(self.bounds()[::2], self.bounds()[1::2]) 237 ] 238 ) 239 240 help_text = "" 241 if self.name: 242 help_text += f"<b> {self.name}:   </b>" 243 help_text += '<b><a href="' + help_url + '" target="_blank">' + library_name + "</a></b>" 244 if self.filename: 245 dots = "" 246 if len(self.filename) > 30: 247 dots = "..." 248 help_text += f"<br/><code><i>({dots}{self.filename[-30:]})</i></code>" 249 250 pdata = "" 251 if self.dataset.GetPointData().GetScalars(): 252 if self.dataset.GetPointData().GetScalars().GetName(): 253 name = self.dataset.GetPointData().GetScalars().GetName() 254 pdata = "<tr><td><b> point data array </b></td><td>" + name + "</td></tr>" 255 256 cdata = "" 257 if self.dataset.GetCellData().GetScalars(): 258 if self.dataset.GetCellData().GetScalars().GetName(): 259 name = self.dataset.GetCellData().GetScalars().GetName() 260 cdata = "<tr><td><b> cell data array </b></td><td>" + name + "</td></tr>" 261 262 pts = self.vertices 263 cm = np.mean(pts, axis=0) 264 265 all = [ 266 "<table>", 267 "<tr>", 268 "<td>", image, "</td>", 269 "<td style='text-align: center; vertical-align: center;'><br/>", help_text, 270 "<table>", 271 "<tr><td><b> bounds </b> <br/> (x/y/z) </td><td>" + str(bounds) + "</td></tr>", 272 "<tr><td><b> center of mass </b></td><td>" + utils.precision(cm,3) + "</td></tr>", 273 # "<tr><td><b> average size </b></td><td>" + str(average_size) + "</td></tr>", 274 "<tr><td><b> nr. points / cells </b></td><td>" 275 + str(self.npoints) + " / " + str(self.ncells) + "</td></tr>", 276 pdata, 277 cdata, 278 "</table>", 279 "</table>", 280 ] 281 return "\n".join(all) 282 283 @property 284 def actor(self): 285 """Return the `vtkActor` of the object.""" 286 # print("building actor") 287 gf = vtk.new("GeometryFilter") 288 gf.SetInputData(self.dataset) 289 gf.Update() 290 out = gf.GetOutput() 291 self.mapper.SetInputData(out) 292 self.mapper.Modified() 293 return self._actor 294 295 @actor.setter 296 def actor(self, _): 297 pass 298 299 def _update(self, data, reset_locators=False): 300 self.dataset = data 301 # self.mapper.SetInputData(data) 302 # self.mapper.Modified() 303 ## self.actor.Modified() 304 return self 305 306 def copy(self, deep=True): 307 """Return a copy of the object. Alias of `clone()`.""" 308 return self.clone(deep=deep) 309 310 def clone(self, deep=True): 311 """Clone the UnstructuredGrid object to yield an exact copy.""" 312 ug = vtk.vtkUnstructuredGrid() 313 if deep: 314 ug.DeepCopy(self.dataset) 315 else: 316 ug.ShallowCopy(self.dataset) 317 if isinstance(self, vedo.UnstructuredGrid): 318 cloned = vedo.UnstructuredGrid(ug) 319 else: 320 cloned = vedo.TetMesh(ug) 321 322 cloned.copy_properties_from(self) 323 324 cloned.pipeline = utils.OperationNode( 325 "clone", parents=[self], shape='diamond', c='#bbe1ed', 326 ) 327 return cloned 328 329 def bounds(self): 330 """ 331 Get the object bounds. 332 Returns a list in format `[xmin,xmax, ymin,ymax, zmin,zmax]`. 333 """ 334 # OVERRIDE CommonAlgorithms.bounds() which is too slow 335 return self.dataset.GetBounds() 336 337 def threshold(self, name=None, above=None, below=None, on="cells"): 338 """ 339 Threshold the tetrahedral mesh by a cell scalar value. 340 Reduce to only tets which satisfy the threshold limits. 341 342 - if `above = below` will only select tets with that specific value. 343 - if `above > below` selection range is flipped. 344 345 Set keyword "on" to either "cells" or "points". 346 """ 347 th = vtk.new("Threshold") 348 th.SetInputData(self.dataset) 349 350 if name is None: 351 if self.celldata.keys(): 352 name = self.celldata.keys()[0] 353 th.SetInputArrayToProcess(0, 0, 0, 1, name) 354 elif self.pointdata.keys(): 355 name = self.pointdata.keys()[0] 356 th.SetInputArrayToProcess(0, 0, 0, 0, name) 357 if name is None: 358 vedo.logger.warning("cannot find active array. Skip.") 359 return self 360 else: 361 if on.startswith("c"): 362 th.SetInputArrayToProcess(0, 0, 0, 1, name) 363 else: 364 th.SetInputArrayToProcess(0, 0, 0, 0, name) 365 366 if above is not None: 367 th.SetLowerThreshold(above) 368 369 if below is not None: 370 th.SetUpperThreshold(below) 371 372 th.Update() 373 return self._update(th.GetOutput()) 374 375 def isosurface(self, value=None, flying_edges=True): 376 """ 377 Return an `Mesh` isosurface extracted from the `Volume` object. 378 379 Set `value` as single float or list of values to draw the isosurface(s). 380 Use flying_edges for faster results (but sometimes can interfere with `smooth()`). 381 382 Examples: 383 - [isosurfaces.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/isosurfaces.py) 384 385  386 """ 387 scrange = self.dataset.GetScalarRange() 388 389 if flying_edges: 390 cf = vtk.new("FlyingEdges3D") 391 cf.InterpolateAttributesOn() 392 else: 393 cf = vtk.new("ContourFilter") 394 cf.UseScalarTreeOn() 395 396 cf.SetInputData(self.dataset) 397 cf.ComputeNormalsOn() 398 399 if utils.is_sequence(value): 400 cf.SetNumberOfContours(len(value)) 401 for i, t in enumerate(value): 402 cf.SetValue(i, t) 403 else: 404 if value is None: 405 value = (2 * scrange[0] + scrange[1]) / 3.0 406 # print("automatic isosurface value =", value) 407 cf.SetValue(0, value) 408 409 cf.Update() 410 poly = cf.GetOutput() 411 412 out = vedo.mesh.Mesh(poly, c=None).phong() 413 out.mapper.SetScalarRange(scrange[0], scrange[1]) 414 415 out.pipeline = utils.OperationNode( 416 "isosurface", 417 parents=[self], 418 comment=f"#pts {out.dataset.GetNumberOfPoints()}", 419 c="#4cc9f0:#e9c46a", 420 ) 421 return out 422 423 def tomesh(self, fill=True, shrink=1.0): 424 """ 425 Build a polygonal `Mesh` from the current object. 426 427 If `fill=True`, the interior faces of all the cells are created. 428 (setting a `shrink` value slightly smaller than the default 1.0 429 can avoid flickering due to internal adjacent faces). 430 431 If `fill=False`, only the boundary faces will be generated. 432 """ 433 gf = vtk.new("GeometryFilter") 434 if fill: 435 sf = vtk.new("ShrinkFilter") 436 sf.SetInputData(self.dataset) 437 sf.SetShrinkFactor(shrink) 438 sf.Update() 439 gf.SetInputData(sf.GetOutput()) 440 gf.Update() 441 poly = gf.GetOutput() 442 else: 443 gf.SetInputData(self.dataset) 444 gf.Update() 445 poly = gf.GetOutput() 446 447 msh = vedo.mesh.Mesh(poly) 448 msh.copy_properties_from(self) 449 450 msh.pipeline = utils.OperationNode( 451 "tomesh", parents=[self], comment=f"fill={fill}", c="#9e2a2b:#e9c46a" 452 ) 453 return msh 454 455 def extract_cell_by_type(self, ctype): 456 """Extract a specific cell type and return a new `UnstructuredGrid`.""" 457 uarr = self.dataset.GetCellTypesArray() 458 ctarrtyp = np.where(utils.vtk2numpy(uarr) == ctype)[0] 459 uarrtyp = utils.numpy2vtk(ctarrtyp, deep=False, dtype="id") 460 selection_node = vtk.new("SelectionNode") 461 selection_node.SetFieldType(vtk.get_class("SelectionNode").CELL) 462 selection_node.SetContentType(vtk.get_class("SelectionNode").INDICES) 463 selection_node.SetSelectionList(uarrtyp) 464 selection = vtk.new("Selection") 465 selection.AddNode(selection_node) 466 es = vtk.new("ExtractSelection") 467 es.SetInputData(0, self.dataset) 468 es.SetInputData(1, selection) 469 es.Update() 470 471 ug = UnstructuredGrid(es.GetOutput()) 472 473 ug.pipeline = utils.OperationNode( 474 "extract_cell_type", comment=f"type {ctype}", 475 c="#edabab", parents=[self], 476 ) 477 return ug 478 479 def extract_cells_by_id(self, idlist, use_point_ids=False): 480 """Return a new `UnstructuredGrid` composed of the specified subset of indices.""" 481 selection_node = vtk.new("SelectionNode") 482 if use_point_ids: 483 selection_node.SetFieldType(vtk.get_class("SelectionNode").POINT) 484 contcells = vtk.get_class("SelectionNode").CONTAINING_CELLS() 485 selection_node.GetProperties().Set(contcells, 1) 486 else: 487 selection_node.SetFieldType(vtk.get_class("SelectionNode").CELL) 488 selection_node.SetContentType(vtk.get_class("SelectionNode").INDICES) 489 vidlist = utils.numpy2vtk(idlist, dtype="id") 490 selection_node.SetSelectionList(vidlist) 491 selection = vtk.new("Selection") 492 selection.AddNode(selection_node) 493 es = vtk.new("ExtractSelection") 494 es.SetInputData(0, self) 495 es.SetInputData(1, selection) 496 es.Update() 497 498 ug = vedo.tetmesh.UnstructuredGrid(es.GetOutput()) 499 pr = vtk.vtkProperty() 500 pr.DeepCopy(self.properties) 501 ug.SetProperty(pr) 502 ug.properties = pr 503 504 ug.mapper.SetLookupTable(utils.ctf2lut(self)) 505 ug.pipeline = utils.OperationNode( 506 "extract_cells_by_id", 507 parents=[self], 508 comment=f"#cells {self.dataset.GetNumberOfCells()}", 509 c="#9e2a2b", 510 ) 511 return ug 512 513 def find_cell(self, p): 514 """Locate the cell that contains a point and return the cell ID.""" 515 cell = vtk.vtkTetra() 516 cell_id = vtk.mutable(0) 517 tol2 = vtk.mutable(0) 518 sub_id = vtk.mutable(0) 519 pcoords = [0, 0, 0] 520 weights = [0, 0, 0] 521 cid = self.dataset.FindCell( 522 p, cell, cell_id, tol2, sub_id, pcoords, weights) 523 return cid 524 525 def clean(self): 526 """ 527 Cleanup unused points and empty cells 528 """ 529 cl = vtk.new("StaticCleanUnstructuredGrid") 530 cl.SetInputData(self.dataset) 531 cl.RemoveUnusedPointsOn() 532 cl.ProduceMergeMapOff() 533 cl.AveragePointDataOff() 534 cl.Update() 535 536 self._update(cl.GetOutput()) 537 self.pipeline = utils.OperationNode( 538 "clean", 539 parents=[self], 540 comment=f"#cells {self.dataset.GetNumberOfCells()}", 541 c="#9e2a2b", 542 ) 543 return self 544 545 def extract_cells_on_plane(self, origin, normal): 546 """ 547 Extract cells that are lying of the specified surface. 548 """ 549 bf = vtk.new("3DLinearGridCrinkleExtractor") 550 bf.SetInputData(self.dataset) 551 bf.CopyPointDataOn() 552 bf.CopyCellDataOn() 553 bf.RemoveUnusedPointsOff() 554 555 plane = vtk.new("Plane") 556 plane.SetOrigin(origin) 557 plane.SetNormal(normal) 558 bf.SetImplicitFunction(plane) 559 bf.Update() 560 561 self._update(bf.GetOutput(), reset_locators=False) 562 self.pipeline = utils.OperationNode( 563 "extract_cells_on_plane", 564 parents=[self], 565 comment=f"#cells {self.dataset.GetNumberOfCells()}", 566 c="#9e2a2b", 567 ) 568 return self 569 570 def extract_cells_on_sphere(self, center, radius): 571 """ 572 Extract cells that are lying of the specified surface. 573 """ 574 bf = vtk.new("3DLinearGridCrinkleExtractor") 575 bf.SetInputData(self.dataset) 576 bf.CopyPointDataOn() 577 bf.CopyCellDataOn() 578 bf.RemoveUnusedPointsOff() 579 580 sph = vtk.new("Sphere") 581 sph.SetRadius(radius) 582 sph.SetCenter(center) 583 bf.SetImplicitFunction(sph) 584 bf.Update() 585 586 self._update(bf.GetOutput()) 587 self.pipeline = utils.OperationNode( 588 "extract_cells_on_sphere", 589 parents=[self], 590 comment=f"#cells {self.dataset.GetNumberOfCells()}", 591 c="#9e2a2b", 592 ) 593 return self 594 595 def extract_cells_on_cylinder(self, center, axis, radius): 596 """ 597 Extract cells that are lying of the specified surface. 598 """ 599 bf = vtk.new("3DLinearGridCrinkleExtractor") 600 bf.SetInputData(self.dataset) 601 bf.CopyPointDataOn() 602 bf.CopyCellDataOn() 603 bf.RemoveUnusedPointsOff() 604 605 cyl = vtk.new("Cylinder") 606 cyl.SetRadius(radius) 607 cyl.SetCenter(center) 608 cyl.SetAxis(axis) 609 bf.SetImplicitFunction(cyl) 610 bf.Update() 611 612 self.pipeline = utils.OperationNode( 613 "extract_cells_on_cylinder", 614 parents=[self], 615 comment=f"#cells {self.dataset.GetNumberOfCells()}", 616 c="#9e2a2b", 617 ) 618 self._update(bf.GetOutput()) 619 return self 620 621 def cut_with_plane(self, origin=(0, 0, 0), normal="x"): 622 """ 623 Cut the object with the plane defined by a point and a normal. 624 625 Arguments: 626 origin : (list) 627 the cutting plane goes through this point 628 normal : (list, str) 629 normal vector to the cutting plane 630 """ 631 # if isinstance(self, vedo.Volume): 632 # raise RuntimeError("cut_with_plane() is not applicable to Volume objects.") 633 634 strn = str(normal) 635 if strn == "x": normal = (1, 0, 0) 636 elif strn == "y": normal = (0, 1, 0) 637 elif strn == "z": normal = (0, 0, 1) 638 elif strn == "-x": normal = (-1, 0, 0) 639 elif strn == "-y": normal = (0, -1, 0) 640 elif strn == "-z": normal = (0, 0, -1) 641 plane = vtk.new("Plane") 642 plane.SetOrigin(origin) 643 plane.SetNormal(normal) 644 clipper = vtk.new("ClipDataSet") 645 clipper.SetInputData(self.dataset) 646 clipper.SetClipFunction(plane) 647 clipper.GenerateClipScalarsOff() 648 clipper.GenerateClippedOutputOff() 649 clipper.SetValue(0) 650 clipper.Update() 651 cout = clipper.GetOutput() 652 653 if isinstance(cout, vtk.vtkUnstructuredGrid): 654 ug = vedo.UnstructuredGrid(cout) 655 if isinstance(self, vedo.UnstructuredGrid): 656 self._update(cout) 657 self.pipeline = utils.OperationNode("cut_with_plane", parents=[self], c="#9e2a2b") 658 return self 659 ug.pipeline = utils.OperationNode("cut_with_plane", parents=[self], c="#9e2a2b") 660 return ug 661 662 else: 663 self._update(cout) 664 self.pipeline = utils.OperationNode("cut_with_plane", parents=[self], c="#9e2a2b") 665 return self 666 667 def cut_with_box(self, box): 668 """ 669 Cut the grid with the specified bounding box. 670 671 Parameter box has format [xmin, xmax, ymin, ymax, zmin, zmax]. 672 If an object is passed, its bounding box are used. 673 674 This method always returns a TetMesh object. 675 676 Example: 677 ```python 678 from vedo import * 679 tetmesh = TetMesh(dataurl+'limb_ugrid.vtk') 680 tetmesh.color('rainbow') 681 cu = Cube(side=500).x(500) # any Mesh works 682 tetmesh.cut_with_box(cu).show(axes=1) 683 ``` 684 685  686 """ 687 bc = vtk.new("BoxClipDataSet") 688 bc.SetInputData(self.dataset) 689 try: 690 boxb = box.bounds() 691 except AttributeError: 692 boxb = box 693 694 bc.SetBoxClip(*boxb) 695 bc.Update() 696 cout = bc.GetOutput() 697 698 # output of vtkBoxClipDataSet is always tetrahedrons 699 tm = vedo.TetMesh(cout) 700 tm.pipeline = utils.OperationNode("cut_with_box", parents=[self], c="#9e2a2b") 701 return tm 702 703 704 def cut_with_mesh( 705 self, mesh, invert=False, whole_cells=False, on_boundary=False 706 ): 707 """ 708 Cut a UnstructuredGrid or TetMesh with a Mesh. 709 710 Use `invert` to return cut off part of the input object. 711 """ 712 ug = self.dataset 713 714 ippd = vtk.new("ImplicitPolyDataDistance") 715 ippd.SetInput(mesh.dataset) 716 717 if whole_cells or on_boundary: 718 clipper = vtk.new("ExtractGeometry") 719 clipper.SetInputData(ug) 720 clipper.SetImplicitFunction(ippd) 721 clipper.SetExtractInside(not invert) 722 clipper.SetExtractBoundaryCells(False) 723 if on_boundary: 724 clipper.SetExtractBoundaryCells(True) 725 clipper.SetExtractOnlyBoundaryCells(True) 726 else: 727 signed_dists = vtk.vtkFloatArray() 728 signed_dists.SetNumberOfComponents(1) 729 signed_dists.SetName("SignedDistance") 730 for pointId in range(ug.GetNumberOfPoints()): 731 p = ug.GetPoint(pointId) 732 signed_dist = ippd.EvaluateFunction(p) 733 signed_dists.InsertNextValue(signed_dist) 734 ug.GetPointData().AddArray(signed_dists) 735 ug.GetPointData().SetActiveScalars("SignedDistance") # NEEDED 736 clipper = vtk.new("ClipDataSet") 737 clipper.SetInputData(ug) 738 clipper.SetInsideOut(not invert) 739 clipper.SetValue(0.0) 740 741 clipper.Update() 742 743 out = vedo.UnstructuredGrid(clipper.GetOutput()) 744 out.pipeline = utils.OperationNode("cut_with_mesh", parents=[self], c="#9e2a2b") 745 return out
Support for UnstructuredGrid objects.
31 def __init__(self, inputobj=None): 32 """ 33 Support for UnstructuredGrid objects. 34 35 Arguments: 36 inputobj : (list, vtkUnstructuredGrid, str) 37 A list in the form `[points, cells, celltypes]`, 38 or a vtkUnstructuredGrid object, or a filename 39 40 Celltypes are identified by the following 41 [convention](https://vtk.org/doc/nightly/html/vtkCellType_8h_source.html). 42 """ 43 super().__init__() 44 45 self.dataset = None 46 47 self.mapper = vtk.new("PolyDataMapper") 48 self._actor = vtk.vtkActor() 49 self._actor.retrieve_object = weak_ref_to(self) 50 self._actor.SetMapper(self.mapper) 51 self.properties = self._actor.GetProperty() 52 53 self.transform = LinearTransform() 54 55 self.name = "UnstructuredGrid" 56 self.filename = "" 57 self.info = {} 58 self.time = 0 59 self.rendered_at = set() 60 self._cmap_name = "" # remember the cmap name for self._keypress 61 62 ################### 63 inputtype = str(type(inputobj)) 64 65 if inputobj is None: 66 self.dataset = vtk.vtkUnstructuredGrid() 67 68 elif utils.is_sequence(inputobj): 69 70 pts, cells, celltypes = inputobj 71 assert len(cells) == len(celltypes) 72 73 self.dataset = vtk.vtkUnstructuredGrid() 74 75 if not utils.is_sequence(cells[0]): 76 tets = [] 77 nf = cells[0] + 1 78 for i, cl in enumerate(cells): 79 if i in (nf, 0): 80 k = i + 1 81 nf = cl + k 82 cell = [cells[j + k] for j in range(cl)] 83 tets.append(cell) 84 cells = tets 85 86 # This would fill the points and use those to define orientation 87 vpts = utils.numpy2vtk(pts, dtype=np.float32) 88 points = vtk.vtkPoints() 89 points.SetData(vpts) 90 self.dataset.SetPoints(points) 91 92 # Fill cells 93 # https://vtk.org/doc/nightly/html/vtkCellType_8h_source.html 94 for i, ct in enumerate(celltypes): 95 cell_conn = cells[i] 96 if ct == vtk.VTK_HEXAHEDRON: 97 cell = vtk.vtkHexahedron() 98 elif ct == vtk.VTK_TETRA: 99 cell = vtk.vtkTetra() 100 elif ct == vtk.VTK_VOXEL: 101 cell = vtk.vtkVoxel() 102 elif ct == vtk.VTK_WEDGE: 103 cell = vtk.vtkWedge() 104 elif ct == vtk.VTK_PYRAMID: 105 cell = vtk.vtkPyramid() 106 elif ct == vtk.VTK_HEXAGONAL_PRISM: 107 cell = vtk.vtkHexagonalPrism() 108 elif ct == vtk.VTK_PENTAGONAL_PRISM: 109 cell = vtk.vtkPentagonalPrism() 110 else: 111 print("UnstructuredGrid: cell type", ct, "not supported. Skip.") 112 continue 113 cpids = cell.GetPointIds() 114 for j, pid in enumerate(cell_conn): 115 cpids.SetId(j, pid) 116 self.dataset.InsertNextCell(ct, cpids) 117 118 elif "UnstructuredGrid" in inputtype: 119 self.dataset = inputobj 120 121 elif isinstance(inputobj, str): 122 if "https://" in inputobj: 123 inputobj = download(inputobj, verbose=False) 124 self.filename = inputobj 125 if inputobj.endswith(".vtu"): 126 reader = vtk.new("XMLUnstructuredGridReader") 127 else: 128 reader = vtk.new("UnstructuredGridReader") 129 self.filename = inputobj 130 reader.SetFileName(inputobj) 131 reader.Update() 132 self.dataset = reader.GetOutput() 133 134 else: 135 vedo.logger.error(f"cannot understand input type {inputtype}") 136 return 137 138 self.pipeline = utils.OperationNode( 139 self, comment=f"#cells {self.dataset.GetNumberOfCells()}", 140 c="#4cc9f0", 141 )
Support for UnstructuredGrid objects.
Arguments:
- inputobj : (list, vtkUnstructuredGrid, str)
A list in the form
[points, cells, celltypes]
, or a vtkUnstructuredGrid object, or a filename
Celltypes are identified by the following convention.
306 def copy(self, deep=True): 307 """Return a copy of the object. Alias of `clone()`.""" 308 return self.clone(deep=deep)
Return a copy of the object. Alias of clone()
.
310 def clone(self, deep=True): 311 """Clone the UnstructuredGrid object to yield an exact copy.""" 312 ug = vtk.vtkUnstructuredGrid() 313 if deep: 314 ug.DeepCopy(self.dataset) 315 else: 316 ug.ShallowCopy(self.dataset) 317 if isinstance(self, vedo.UnstructuredGrid): 318 cloned = vedo.UnstructuredGrid(ug) 319 else: 320 cloned = vedo.TetMesh(ug) 321 322 cloned.copy_properties_from(self) 323 324 cloned.pipeline = utils.OperationNode( 325 "clone", parents=[self], shape='diamond', c='#bbe1ed', 326 ) 327 return cloned
Clone the UnstructuredGrid object to yield an exact copy.
329 def bounds(self): 330 """ 331 Get the object bounds. 332 Returns a list in format `[xmin,xmax, ymin,ymax, zmin,zmax]`. 333 """ 334 # OVERRIDE CommonAlgorithms.bounds() which is too slow 335 return self.dataset.GetBounds()
Get the object bounds.
Returns a list in format [xmin,xmax, ymin,ymax, zmin,zmax]
.
337 def threshold(self, name=None, above=None, below=None, on="cells"): 338 """ 339 Threshold the tetrahedral mesh by a cell scalar value. 340 Reduce to only tets which satisfy the threshold limits. 341 342 - if `above = below` will only select tets with that specific value. 343 - if `above > below` selection range is flipped. 344 345 Set keyword "on" to either "cells" or "points". 346 """ 347 th = vtk.new("Threshold") 348 th.SetInputData(self.dataset) 349 350 if name is None: 351 if self.celldata.keys(): 352 name = self.celldata.keys()[0] 353 th.SetInputArrayToProcess(0, 0, 0, 1, name) 354 elif self.pointdata.keys(): 355 name = self.pointdata.keys()[0] 356 th.SetInputArrayToProcess(0, 0, 0, 0, name) 357 if name is None: 358 vedo.logger.warning("cannot find active array. Skip.") 359 return self 360 else: 361 if on.startswith("c"): 362 th.SetInputArrayToProcess(0, 0, 0, 1, name) 363 else: 364 th.SetInputArrayToProcess(0, 0, 0, 0, name) 365 366 if above is not None: 367 th.SetLowerThreshold(above) 368 369 if below is not None: 370 th.SetUpperThreshold(below) 371 372 th.Update() 373 return self._update(th.GetOutput())
Threshold the tetrahedral mesh by a cell scalar value. Reduce to only tets which satisfy the threshold limits.
- if
above = below
will only select tets with that specific value. - if
above > below
selection range is flipped.
Set keyword "on" to either "cells" or "points".
375 def isosurface(self, value=None, flying_edges=True): 376 """ 377 Return an `Mesh` isosurface extracted from the `Volume` object. 378 379 Set `value` as single float or list of values to draw the isosurface(s). 380 Use flying_edges for faster results (but sometimes can interfere with `smooth()`). 381 382 Examples: 383 - [isosurfaces.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/isosurfaces.py) 384 385  386 """ 387 scrange = self.dataset.GetScalarRange() 388 389 if flying_edges: 390 cf = vtk.new("FlyingEdges3D") 391 cf.InterpolateAttributesOn() 392 else: 393 cf = vtk.new("ContourFilter") 394 cf.UseScalarTreeOn() 395 396 cf.SetInputData(self.dataset) 397 cf.ComputeNormalsOn() 398 399 if utils.is_sequence(value): 400 cf.SetNumberOfContours(len(value)) 401 for i, t in enumerate(value): 402 cf.SetValue(i, t) 403 else: 404 if value is None: 405 value = (2 * scrange[0] + scrange[1]) / 3.0 406 # print("automatic isosurface value =", value) 407 cf.SetValue(0, value) 408 409 cf.Update() 410 poly = cf.GetOutput() 411 412 out = vedo.mesh.Mesh(poly, c=None).phong() 413 out.mapper.SetScalarRange(scrange[0], scrange[1]) 414 415 out.pipeline = utils.OperationNode( 416 "isosurface", 417 parents=[self], 418 comment=f"#pts {out.dataset.GetNumberOfPoints()}", 419 c="#4cc9f0:#e9c46a", 420 ) 421 return out
Return an Mesh
isosurface extracted from the Volume
object.
Set value
as single float or list of values to draw the isosurface(s).
Use flying_edges for faster results (but sometimes can interfere with smooth()
).
Examples:
423 def tomesh(self, fill=True, shrink=1.0): 424 """ 425 Build a polygonal `Mesh` from the current object. 426 427 If `fill=True`, the interior faces of all the cells are created. 428 (setting a `shrink` value slightly smaller than the default 1.0 429 can avoid flickering due to internal adjacent faces). 430 431 If `fill=False`, only the boundary faces will be generated. 432 """ 433 gf = vtk.new("GeometryFilter") 434 if fill: 435 sf = vtk.new("ShrinkFilter") 436 sf.SetInputData(self.dataset) 437 sf.SetShrinkFactor(shrink) 438 sf.Update() 439 gf.SetInputData(sf.GetOutput()) 440 gf.Update() 441 poly = gf.GetOutput() 442 else: 443 gf.SetInputData(self.dataset) 444 gf.Update() 445 poly = gf.GetOutput() 446 447 msh = vedo.mesh.Mesh(poly) 448 msh.copy_properties_from(self) 449 450 msh.pipeline = utils.OperationNode( 451 "tomesh", parents=[self], comment=f"fill={fill}", c="#9e2a2b:#e9c46a" 452 ) 453 return msh
Build a polygonal Mesh
from the current object.
If fill=True
, the interior faces of all the cells are created.
(setting a shrink
value slightly smaller than the default 1.0
can avoid flickering due to internal adjacent faces).
If fill=False
, only the boundary faces will be generated.
455 def extract_cell_by_type(self, ctype): 456 """Extract a specific cell type and return a new `UnstructuredGrid`.""" 457 uarr = self.dataset.GetCellTypesArray() 458 ctarrtyp = np.where(utils.vtk2numpy(uarr) == ctype)[0] 459 uarrtyp = utils.numpy2vtk(ctarrtyp, deep=False, dtype="id") 460 selection_node = vtk.new("SelectionNode") 461 selection_node.SetFieldType(vtk.get_class("SelectionNode").CELL) 462 selection_node.SetContentType(vtk.get_class("SelectionNode").INDICES) 463 selection_node.SetSelectionList(uarrtyp) 464 selection = vtk.new("Selection") 465 selection.AddNode(selection_node) 466 es = vtk.new("ExtractSelection") 467 es.SetInputData(0, self.dataset) 468 es.SetInputData(1, selection) 469 es.Update() 470 471 ug = UnstructuredGrid(es.GetOutput()) 472 473 ug.pipeline = utils.OperationNode( 474 "extract_cell_type", comment=f"type {ctype}", 475 c="#edabab", parents=[self], 476 ) 477 return ug
Extract a specific cell type and return a new UnstructuredGrid
.
479 def extract_cells_by_id(self, idlist, use_point_ids=False): 480 """Return a new `UnstructuredGrid` composed of the specified subset of indices.""" 481 selection_node = vtk.new("SelectionNode") 482 if use_point_ids: 483 selection_node.SetFieldType(vtk.get_class("SelectionNode").POINT) 484 contcells = vtk.get_class("SelectionNode").CONTAINING_CELLS() 485 selection_node.GetProperties().Set(contcells, 1) 486 else: 487 selection_node.SetFieldType(vtk.get_class("SelectionNode").CELL) 488 selection_node.SetContentType(vtk.get_class("SelectionNode").INDICES) 489 vidlist = utils.numpy2vtk(idlist, dtype="id") 490 selection_node.SetSelectionList(vidlist) 491 selection = vtk.new("Selection") 492 selection.AddNode(selection_node) 493 es = vtk.new("ExtractSelection") 494 es.SetInputData(0, self) 495 es.SetInputData(1, selection) 496 es.Update() 497 498 ug = vedo.tetmesh.UnstructuredGrid(es.GetOutput()) 499 pr = vtk.vtkProperty() 500 pr.DeepCopy(self.properties) 501 ug.SetProperty(pr) 502 ug.properties = pr 503 504 ug.mapper.SetLookupTable(utils.ctf2lut(self)) 505 ug.pipeline = utils.OperationNode( 506 "extract_cells_by_id", 507 parents=[self], 508 comment=f"#cells {self.dataset.GetNumberOfCells()}", 509 c="#9e2a2b", 510 ) 511 return ug
Return a new UnstructuredGrid
composed of the specified subset of indices.
513 def find_cell(self, p): 514 """Locate the cell that contains a point and return the cell ID.""" 515 cell = vtk.vtkTetra() 516 cell_id = vtk.mutable(0) 517 tol2 = vtk.mutable(0) 518 sub_id = vtk.mutable(0) 519 pcoords = [0, 0, 0] 520 weights = [0, 0, 0] 521 cid = self.dataset.FindCell( 522 p, cell, cell_id, tol2, sub_id, pcoords, weights) 523 return cid
Locate the cell that contains a point and return the cell ID.
525 def clean(self): 526 """ 527 Cleanup unused points and empty cells 528 """ 529 cl = vtk.new("StaticCleanUnstructuredGrid") 530 cl.SetInputData(self.dataset) 531 cl.RemoveUnusedPointsOn() 532 cl.ProduceMergeMapOff() 533 cl.AveragePointDataOff() 534 cl.Update() 535 536 self._update(cl.GetOutput()) 537 self.pipeline = utils.OperationNode( 538 "clean", 539 parents=[self], 540 comment=f"#cells {self.dataset.GetNumberOfCells()}", 541 c="#9e2a2b", 542 ) 543 return self
Cleanup unused points and empty cells
545 def extract_cells_on_plane(self, origin, normal): 546 """ 547 Extract cells that are lying of the specified surface. 548 """ 549 bf = vtk.new("3DLinearGridCrinkleExtractor") 550 bf.SetInputData(self.dataset) 551 bf.CopyPointDataOn() 552 bf.CopyCellDataOn() 553 bf.RemoveUnusedPointsOff() 554 555 plane = vtk.new("Plane") 556 plane.SetOrigin(origin) 557 plane.SetNormal(normal) 558 bf.SetImplicitFunction(plane) 559 bf.Update() 560 561 self._update(bf.GetOutput(), reset_locators=False) 562 self.pipeline = utils.OperationNode( 563 "extract_cells_on_plane", 564 parents=[self], 565 comment=f"#cells {self.dataset.GetNumberOfCells()}", 566 c="#9e2a2b", 567 ) 568 return self
Extract cells that are lying of the specified surface.
570 def extract_cells_on_sphere(self, center, radius): 571 """ 572 Extract cells that are lying of the specified surface. 573 """ 574 bf = vtk.new("3DLinearGridCrinkleExtractor") 575 bf.SetInputData(self.dataset) 576 bf.CopyPointDataOn() 577 bf.CopyCellDataOn() 578 bf.RemoveUnusedPointsOff() 579 580 sph = vtk.new("Sphere") 581 sph.SetRadius(radius) 582 sph.SetCenter(center) 583 bf.SetImplicitFunction(sph) 584 bf.Update() 585 586 self._update(bf.GetOutput()) 587 self.pipeline = utils.OperationNode( 588 "extract_cells_on_sphere", 589 parents=[self], 590 comment=f"#cells {self.dataset.GetNumberOfCells()}", 591 c="#9e2a2b", 592 ) 593 return self
Extract cells that are lying of the specified surface.
595 def extract_cells_on_cylinder(self, center, axis, radius): 596 """ 597 Extract cells that are lying of the specified surface. 598 """ 599 bf = vtk.new("3DLinearGridCrinkleExtractor") 600 bf.SetInputData(self.dataset) 601 bf.CopyPointDataOn() 602 bf.CopyCellDataOn() 603 bf.RemoveUnusedPointsOff() 604 605 cyl = vtk.new("Cylinder") 606 cyl.SetRadius(radius) 607 cyl.SetCenter(center) 608 cyl.SetAxis(axis) 609 bf.SetImplicitFunction(cyl) 610 bf.Update() 611 612 self.pipeline = utils.OperationNode( 613 "extract_cells_on_cylinder", 614 parents=[self], 615 comment=f"#cells {self.dataset.GetNumberOfCells()}", 616 c="#9e2a2b", 617 ) 618 self._update(bf.GetOutput()) 619 return self
Extract cells that are lying of the specified surface.
621 def cut_with_plane(self, origin=(0, 0, 0), normal="x"): 622 """ 623 Cut the object with the plane defined by a point and a normal. 624 625 Arguments: 626 origin : (list) 627 the cutting plane goes through this point 628 normal : (list, str) 629 normal vector to the cutting plane 630 """ 631 # if isinstance(self, vedo.Volume): 632 # raise RuntimeError("cut_with_plane() is not applicable to Volume objects.") 633 634 strn = str(normal) 635 if strn == "x": normal = (1, 0, 0) 636 elif strn == "y": normal = (0, 1, 0) 637 elif strn == "z": normal = (0, 0, 1) 638 elif strn == "-x": normal = (-1, 0, 0) 639 elif strn == "-y": normal = (0, -1, 0) 640 elif strn == "-z": normal = (0, 0, -1) 641 plane = vtk.new("Plane") 642 plane.SetOrigin(origin) 643 plane.SetNormal(normal) 644 clipper = vtk.new("ClipDataSet") 645 clipper.SetInputData(self.dataset) 646 clipper.SetClipFunction(plane) 647 clipper.GenerateClipScalarsOff() 648 clipper.GenerateClippedOutputOff() 649 clipper.SetValue(0) 650 clipper.Update() 651 cout = clipper.GetOutput() 652 653 if isinstance(cout, vtk.vtkUnstructuredGrid): 654 ug = vedo.UnstructuredGrid(cout) 655 if isinstance(self, vedo.UnstructuredGrid): 656 self._update(cout) 657 self.pipeline = utils.OperationNode("cut_with_plane", parents=[self], c="#9e2a2b") 658 return self 659 ug.pipeline = utils.OperationNode("cut_with_plane", parents=[self], c="#9e2a2b") 660 return ug 661 662 else: 663 self._update(cout) 664 self.pipeline = utils.OperationNode("cut_with_plane", parents=[self], c="#9e2a2b") 665 return self
Cut the object with the plane defined by a point and a normal.
Arguments:
- origin : (list) the cutting plane goes through this point
- normal : (list, str) normal vector to the cutting plane
667 def cut_with_box(self, box): 668 """ 669 Cut the grid with the specified bounding box. 670 671 Parameter box has format [xmin, xmax, ymin, ymax, zmin, zmax]. 672 If an object is passed, its bounding box are used. 673 674 This method always returns a TetMesh object. 675 676 Example: 677 ```python 678 from vedo import * 679 tetmesh = TetMesh(dataurl+'limb_ugrid.vtk') 680 tetmesh.color('rainbow') 681 cu = Cube(side=500).x(500) # any Mesh works 682 tetmesh.cut_with_box(cu).show(axes=1) 683 ``` 684 685  686 """ 687 bc = vtk.new("BoxClipDataSet") 688 bc.SetInputData(self.dataset) 689 try: 690 boxb = box.bounds() 691 except AttributeError: 692 boxb = box 693 694 bc.SetBoxClip(*boxb) 695 bc.Update() 696 cout = bc.GetOutput() 697 698 # output of vtkBoxClipDataSet is always tetrahedrons 699 tm = vedo.TetMesh(cout) 700 tm.pipeline = utils.OperationNode("cut_with_box", parents=[self], c="#9e2a2b") 701 return tm
Cut the grid with the specified bounding box.
Parameter box has format [xmin, xmax, ymin, ymax, zmin, zmax]. If an object is passed, its bounding box are used.
This method always returns a TetMesh object.
Example:
from vedo import *
tetmesh = TetMesh(dataurl+'limb_ugrid.vtk')
tetmesh.color('rainbow')
cu = Cube(side=500).x(500) # any Mesh works
tetmesh.cut_with_box(cu).show(axes=1)
704 def cut_with_mesh( 705 self, mesh, invert=False, whole_cells=False, on_boundary=False 706 ): 707 """ 708 Cut a UnstructuredGrid or TetMesh with a Mesh. 709 710 Use `invert` to return cut off part of the input object. 711 """ 712 ug = self.dataset 713 714 ippd = vtk.new("ImplicitPolyDataDistance") 715 ippd.SetInput(mesh.dataset) 716 717 if whole_cells or on_boundary: 718 clipper = vtk.new("ExtractGeometry") 719 clipper.SetInputData(ug) 720 clipper.SetImplicitFunction(ippd) 721 clipper.SetExtractInside(not invert) 722 clipper.SetExtractBoundaryCells(False) 723 if on_boundary: 724 clipper.SetExtractBoundaryCells(True) 725 clipper.SetExtractOnlyBoundaryCells(True) 726 else: 727 signed_dists = vtk.vtkFloatArray() 728 signed_dists.SetNumberOfComponents(1) 729 signed_dists.SetName("SignedDistance") 730 for pointId in range(ug.GetNumberOfPoints()): 731 p = ug.GetPoint(pointId) 732 signed_dist = ippd.EvaluateFunction(p) 733 signed_dists.InsertNextValue(signed_dist) 734 ug.GetPointData().AddArray(signed_dists) 735 ug.GetPointData().SetActiveScalars("SignedDistance") # NEEDED 736 clipper = vtk.new("ClipDataSet") 737 clipper.SetInputData(ug) 738 clipper.SetInsideOut(not invert) 739 clipper.SetValue(0.0) 740 741 clipper.Update() 742 743 out = vedo.UnstructuredGrid(clipper.GetOutput()) 744 out.pipeline = utils.OperationNode("cut_with_mesh", parents=[self], c="#9e2a2b") 745 return out
Cut a UnstructuredGrid or TetMesh with a Mesh.
Use invert
to return cut off part of the input object.
Inherited Members
- vedo.core.PointAlgorithms
- apply_transform
- apply_transform_from_actor
- pos
- shift
- x
- y
- z
- rotate
- rotate_x
- rotate_y
- rotate_z
- reorient
- scale
- vedo.core.CommonAlgorithms
- pointdata
- celldata
- metadata
- memory_address
- memory_size
- modified
- box
- xbounds
- ybounds
- zbounds
- diagonal_size
- average_size
- center_of_mass
- copy_data_from
- inputdata
- npoints
- nvertices
- ncells
- points
- cell_centers
- lines
- lines_as_flat_array
- mark_boundaries
- find_cells_in_bounds
- find_cells_along_line
- find_cells_along_plane
- delete_cells_by_point_index
- map_cells_to_points
- vertices
- coordinates
- cells_as_flat_array
- cells
- map_points_to_cells
- resample_data_from
- interpolate_data_from
- add_ids
- gradient
- divergence
- vorticity
- probe
- compute_cell_size
- integrate_arrays_over_domain
- write
- shrink
- vedo.visual.MeshVisual
- follow_camera
- wireframe
- flat
- phong
- backface_culling
- render_lines_as_tubes
- frontface_culling
- backcolor
- bc
- linewidth
- lw
- linecolor
- lc
- texture
- vedo.visual.PointsVisual
- clone2d
- copy_properties_from
- color
- c
- alpha
- lut_color_at
- opacity
- force_opaque
- force_translucent
- point_size
- ps
- render_points_as_spheres
- lighting
- point_blurring
- cellcolors
- pointcolors
- cmap
- add_trail
- update_trail
- add_shadow
- update_shadows
- labels
- labels2d
- legend
- flagpole
- flagpost
- caption
748class TetMesh(UnstructuredGrid): 749 """The class describing tetrahedral meshes.""" 750 751 def __init__(self, inputobj=None): 752 """ 753 Arguments: 754 inputobj : (vtkDataSet, list, str) 755 list of points and tet indices, or filename 756 """ 757 super().__init__() 758 759 self.dataset = None 760 761 self.mapper = vtk.new("PolyDataMapper") 762 self._actor = vtk.vtkActor() 763 self._actor.retrieve_object = weak_ref_to(self) 764 self._actor.SetMapper(self.mapper) 765 self.properties = self._actor.GetProperty() 766 767 self.transform = LinearTransform() 768 769 self.name = "TetMesh" 770 self.filename = "" 771 772 773 # inputtype = str(type(inputobj)) 774 # print('TetMesh inputtype', inputtype) 775 776 ################### 777 if inputobj is None: 778 self.dataset = vtk.vtkUnstructuredGrid() 779 780 elif isinstance(inputobj, vtk.vtkUnstructuredGrid): 781 self.dataset = inputobj 782 783 elif isinstance(inputobj, UnstructuredGrid): 784 self.dataset = inputobj.dataset 785 786 elif isinstance(inputobj, vtk.vtkRectilinearGrid): 787 r2t = vtk.new("RectilinearGridToTetrahedra") 788 r2t.SetInputData(inputobj) 789 r2t.RememberVoxelIdOn() 790 r2t.SetTetraPerCellTo6() 791 r2t.Update() 792 self.dataset = r2t.GetOutput() 793 794 elif isinstance(inputobj, vtk.vtkDataSet): 795 r2t = vtk.new("DataSetTriangleFilter") 796 r2t.SetInputData(inputobj) 797 r2t.TetrahedraOnlyOn() 798 r2t.Update() 799 self.dataset = r2t.GetOutput() 800 801 elif isinstance(inputobj, str): 802 if "https://" in inputobj: 803 inputobj = download(inputobj, verbose=False) 804 if inputobj.endswith(".vtu"): 805 reader = vtk.new("XMLUnstructuredGridReader") 806 else: 807 reader = vtk.new("UnstructuredGridReader") 808 self.filename = inputobj 809 reader.SetFileName(inputobj) 810 reader.Update() 811 ug = reader.GetOutput() 812 tt = vtk.new("DataSetTriangleFilter") 813 tt.SetInputData(ug) 814 tt.SetTetrahedraOnly(True) 815 tt.Update() 816 self.dataset = tt.GetOutput() 817 818 elif utils.is_sequence(inputobj): 819 self.dataset = vtk.vtkUnstructuredGrid() 820 821 points, cells = inputobj 822 if len(points) == 0: 823 return 824 if not utils.is_sequence(points[0]): 825 return 826 if len(cells) == 0: 827 return 828 829 if not utils.is_sequence(cells[0]): 830 tets = [] 831 nf = cells[0] + 1 832 for i, cl in enumerate(cells): 833 if i in (nf, 0): 834 k = i + 1 835 nf = cl + k 836 cell = [cells[j + k] for j in range(cl)] 837 tets.append(cell) 838 cells = tets 839 840 source_points = vtk.vtkPoints() 841 varr = utils.numpy2vtk(points, dtype=np.float32) 842 source_points.SetData(varr) 843 self.dataset.SetPoints(source_points) 844 845 source_tets = vtk.vtkCellArray() 846 for f in cells: 847 ele = vtk.vtkTetra() 848 pid = ele.GetPointIds() 849 for i, fi in enumerate(f): 850 pid.SetId(i, fi) 851 source_tets.InsertNextCell(ele) 852 self.dataset.SetCells(vtk.VTK_TETRA, source_tets) 853 854 else: 855 vedo.logger.error(f"cannot understand input type {inputtype}") 856 return 857 858 self.pipeline = utils.OperationNode( 859 self, comment=f"#tets {self.dataset.GetNumberOfCells()}", 860 c="#9e2a2b", 861 ) 862 863 ################################################################## 864 def __str__(self): 865 """Print a string summary of the `TetMesh` object.""" 866 module = self.__class__.__module__ 867 name = self.__class__.__name__ 868 out = vedo.printc( 869 f"{module}.{name} at ({hex(self.memory_address())})".ljust(75), 870 c="c", bold=True, invert=True, return_string=True, 871 ) 872 out += "\x1b[0m\u001b[36m" 873 874 out += "nr. of verts".ljust(14) + ": " + str(self.npoints) + "\n" 875 out += "nr. of tetras".ljust(14)+ ": " + str(self.ncells) + "\n" 876 877 if self.npoints: 878 out+="size".ljust(14)+ ": average=" + utils.precision(self.average_size(),6) 879 out+=", diagonal="+ utils.precision(self.diagonal_size(), 6)+ "\n" 880 out+="center of mass".ljust(14) + ": " + utils.precision(self.center_of_mass(),6)+"\n" 881 882 bnds = self.bounds() 883 bx1, bx2 = utils.precision(bnds[0], 3), utils.precision(bnds[1], 3) 884 by1, by2 = utils.precision(bnds[2], 3), utils.precision(bnds[3], 3) 885 bz1, bz2 = utils.precision(bnds[4], 3), utils.precision(bnds[5], 3) 886 out += "bounds".ljust(14) + ":" 887 out += " x=(" + bx1 + ", " + bx2 + ")," 888 out += " y=(" + by1 + ", " + by2 + ")," 889 out += " z=(" + bz1 + ", " + bz2 + ")\n" 890 891 for key in self.pointdata.keys(): 892 arr = self.pointdata[key] 893 rng = utils.precision(arr.min(), 3) + ", " + utils.precision(arr.max(), 3) 894 mark_active = "pointdata" 895 a_scalars = self.dataset.GetPointData().GetScalars() 896 a_vectors = self.dataset.GetPointData().GetVectors() 897 a_tensors = self.dataset.GetPointData().GetTensors() 898 if a_scalars and a_scalars.GetName() == key: 899 mark_active += " *" 900 elif a_vectors and a_vectors.GetName() == key: 901 mark_active += " **" 902 elif a_tensors and a_tensors.GetName() == key: 903 mark_active += " ***" 904 out += mark_active.ljust(14) + f': "{key}" ({arr.dtype}), ndim={arr.ndim}' 905 out += f", range=({rng})\n" 906 907 for key in self.celldata.keys(): 908 arr = self.celldata[key] 909 rng = utils.precision(arr.min(), 3) + ", " + utils.precision(arr.max(), 3) 910 mark_active = "celldata" 911 a_scalars = self.dataset.GetCellData().GetScalars() 912 a_vectors = self.dataset.GetCellData().GetVectors() 913 a_tensors = self.dataset.GetCellData().GetTensors() 914 if a_scalars and a_scalars.GetName() == key: 915 mark_active += " *" 916 elif a_vectors and a_vectors.GetName() == key: 917 mark_active += " **" 918 elif a_tensors and a_tensors.GetName() == key: 919 mark_active += " ***" 920 out += mark_active.ljust(14) + f': "{key}" ({arr.dtype}), ndim={arr.ndim}' 921 out += f", range=({rng})\n" 922 923 for key in self.metadata.keys(): 924 arr = self.metadata[key] 925 out+= "metadata".ljust(14) + ": " + f'"{key}" ({len(arr)} values)\n' 926 927 return out.rstrip() + "\x1b[0m" 928 929 930 def _repr_html_(self): 931 """ 932 HTML representation of the TetMesh object for Jupyter Notebooks. 933 934 Returns: 935 HTML text with the image and some properties. 936 """ 937 import io 938 import base64 939 from PIL import Image 940 941 library_name = "vedo.tetmesh.TetMesh" 942 help_url = "https://vedo.embl.es/docs/vedo/tetmesh.html" 943 944 arr = self.thumbnail() 945 im = Image.fromarray(arr) 946 buffered = io.BytesIO() 947 im.save(buffered, format="PNG", quality=100) 948 encoded = base64.b64encode(buffered.getvalue()).decode("utf-8") 949 url = "data:image/png;base64," + encoded 950 image = f"<img src='{url}'></img>" 951 952 bounds = "<br/>".join( 953 [ 954 utils.precision(min_x,4) + " ... " + utils.precision(max_x,4) 955 for min_x, max_x in zip(self.bounds()[::2], self.bounds()[1::2]) 956 ] 957 ) 958 959 help_text = "" 960 if self.name: 961 help_text += f"<b> {self.name}:   </b>" 962 help_text += '<b><a href="' + help_url + '" target="_blank">' + library_name + "</a></b>" 963 if self.filename: 964 dots = "" 965 if len(self.filename) > 30: 966 dots = "..." 967 help_text += f"<br/><code><i>({dots}{self.filename[-30:]})</i></code>" 968 969 pdata = "" 970 if self.dataset.GetPointData().GetScalars(): 971 if self.dataset.GetPointData().GetScalars().GetName(): 972 name = self.dataset.GetPointData().GetScalars().GetName() 973 pdata = "<tr><td><b> point data array </b></td><td>" + name + "</td></tr>" 974 975 cdata = "" 976 if self.dataset.GetCellData().GetScalars(): 977 if self.dataset.GetCellData().GetScalars().GetName(): 978 name = self.dataset.GetCellData().GetScalars().GetName() 979 cdata = "<tr><td><b> cell data array </b></td><td>" + name + "</td></tr>" 980 981 pts = self.vertices 982 cm = np.mean(pts, axis=0) 983 984 allt = [ 985 "<table>", 986 "<tr>", 987 "<td>", image, "</td>", 988 "<td style='text-align: center; vertical-align: center;'><br/>", help_text, 989 "<table>", 990 "<tr><td><b> bounds </b> <br/> (x/y/z) </td><td>" + str(bounds) + "</td></tr>", 991 "<tr><td><b> center of mass </b></td><td>" + utils.precision(cm,3) + "</td></tr>", 992 "<tr><td><b> nr. points / tets </b></td><td>" 993 + str(self.npoints) + " / " + str(self.ncells) + "</td></tr>", 994 pdata, 995 cdata, 996 "</table>", 997 "</table>", 998 ] 999 return "\n".join(allt) 1000 1001 def compute_quality(self, metric=7): 1002 """ 1003 Calculate functions of quality for the elements of a triangular mesh. 1004 This method adds to the mesh a cell array named "Quality". 1005 1006 Arguments: 1007 metric : (int) 1008 type of estimators: 1009 - EDGE RATIO, 0 1010 - ASPECT RATIO, 1 1011 - RADIUS RATIO, 2 1012 - ASPECT FROBENIUS, 3 1013 - MIN_ANGLE, 4 1014 - COLLAPSE RATIO, 5 1015 - ASPECT GAMMA, 6 1016 - VOLUME, 7 1017 - ... 1018 1019 See class [vtkMeshQuality](https://vtk.org/doc/nightly/html/classvtkMeshQuality.html) 1020 for an explanation of the meaning of each metric.. 1021 """ 1022 qf = vtk.new("MeshQuality") 1023 qf.SetInputData(self.dataset) 1024 qf.SetTetQualityMeasure(metric) 1025 qf.SaveCellQualityOn() 1026 qf.Update() 1027 self._update(qf.GetOutput()) 1028 return utils.vtk2numpy(qf.GetOutput().GetCellData().GetArray("Quality")) 1029 1030 def check_validity(self, tol=0): 1031 """ 1032 Return an array of possible problematic tets following this convention: 1033 ```python 1034 Valid = 0 1035 WrongNumberOfPoints = 01 1036 IntersectingEdges = 02 1037 IntersectingFaces = 04 1038 NoncontiguousEdges = 08 1039 Nonconvex = 10 1040 OrientedIncorrectly = 20 1041 ``` 1042 1043 Arguments: 1044 tol : (float) 1045 This value is used as an epsilon for floating point 1046 equality checks throughout the cell checking process. 1047 """ 1048 vald = vtk.new("CellValidator") 1049 if tol: 1050 vald.SetTolerance(tol) 1051 vald.SetInputData(self.dataset) 1052 vald.Update() 1053 varr = vald.GetOutput().GetCellData().GetArray("ValidityState") 1054 return utils.vtk2numpy(varr) 1055 1056 def decimate(self, scalars_name, fraction=0.5, n=0): 1057 """ 1058 Downsample the number of tets in a TetMesh to a specified fraction. 1059 Either `fraction` or `n` must be set. 1060 1061 Arguments: 1062 fraction : (float) 1063 the desired final fraction of the total. 1064 n : (int) 1065 the desired number of final tets 1066 1067 .. note:: setting `fraction=0.1` leaves 10% of the original nr of tets. 1068 """ 1069 decimate = vtk.new("UnstructuredGridQuadricDecimation") 1070 decimate.SetInputData(self.dataset) 1071 decimate.SetScalarsName(scalars_name) 1072 1073 if n: # n = desired number of points 1074 decimate.SetNumberOfTetsOutput(n) 1075 else: 1076 decimate.SetTargetReduction(1 - fraction) 1077 decimate.Update() 1078 self._update(decimate.GetOutput()) 1079 self.pipeline = utils.OperationNode( 1080 "decimate", comment=f"array: {scalars_name}", 1081 c="#edabab", parents=[self], 1082 ) 1083 return self 1084 1085 def subdvide(self): 1086 """ 1087 Increase the number of tets of a `TetMesh`. 1088 Subdivide one tetrahedron into twelve for every tetra. 1089 """ 1090 sd = vtk.new("SubdivideTetra") 1091 sd.SetInputData(self.dataset) 1092 sd.Update() 1093 self._update(sd.GetOutput()) 1094 self.pipeline = utils.OperationNode( 1095 "subdvide", c="#edabab", parents=[self], 1096 ) 1097 return self 1098 1099 def isosurface(self, value=True): 1100 """ 1101 Return a `vedo.Mesh` isosurface. 1102 1103 Set `value` to a single value or list of values to compute the isosurface(s). 1104 """ 1105 if not self.dataset.GetPointData().GetScalars(): 1106 self.map_cells_to_points() 1107 scrange = self.dataset.GetPointData().GetScalars().GetRange() 1108 cf = vtk.new("ContourFilter") # vtk.new("ContourGrid") 1109 cf.SetInputData(self.dataset) 1110 1111 if utils.is_sequence(value): 1112 cf.SetNumberOfContours(len(value)) 1113 for i, t in enumerate(value): 1114 cf.SetValue(i, t) 1115 cf.Update() 1116 else: 1117 if value is True: 1118 value = (2 * scrange[0] + scrange[1]) / 3.0 1119 cf.SetValue(0, value) 1120 cf.Update() 1121 1122 clp = vtk.new("CleanPolyData") 1123 clp.SetInputData(cf.GetOutput()) 1124 clp.Update() 1125 msh = Mesh(clp.GetOutput(), c=None).phong() 1126 msh.copy_properties_from(self) 1127 msh.pipeline = utils.OperationNode("isosurface", c="#edabab", parents=[self]) 1128 return msh 1129 1130 def slice(self, origin=(0, 0, 0), normal=(1, 0, 0)): 1131 """ 1132 Return a 2D slice of the mesh by a plane passing through origin and assigned normal. 1133 """ 1134 strn = str(normal) 1135 if strn == "x": normal = (1, 0, 0) 1136 elif strn == "y": normal = (0, 1, 0) 1137 elif strn == "z": normal = (0, 0, 1) 1138 elif strn == "-x": normal = (-1, 0, 0) 1139 elif strn == "-y": normal = (0, -1, 0) 1140 elif strn == "-z": normal = (0, 0, -1) 1141 plane = vtk.new("Plane") 1142 plane.SetOrigin(origin) 1143 plane.SetNormal(normal) 1144 1145 cc = vtk.new("Cutter") 1146 cc.SetInputData(self.dataset) 1147 cc.SetCutFunction(plane) 1148 cc.Update() 1149 msh = Mesh(cc.GetOutput()).flat().lighting("ambient") 1150 msh.copy_properties_from(self) 1151 msh.pipeline = utils.OperationNode("slice", c="#edabab", parents=[self]) 1152 return msh
The class describing tetrahedral meshes.
751 def __init__(self, inputobj=None): 752 """ 753 Arguments: 754 inputobj : (vtkDataSet, list, str) 755 list of points and tet indices, or filename 756 """ 757 super().__init__() 758 759 self.dataset = None 760 761 self.mapper = vtk.new("PolyDataMapper") 762 self._actor = vtk.vtkActor() 763 self._actor.retrieve_object = weak_ref_to(self) 764 self._actor.SetMapper(self.mapper) 765 self.properties = self._actor.GetProperty() 766 767 self.transform = LinearTransform() 768 769 self.name = "TetMesh" 770 self.filename = "" 771 772 773 # inputtype = str(type(inputobj)) 774 # print('TetMesh inputtype', inputtype) 775 776 ################### 777 if inputobj is None: 778 self.dataset = vtk.vtkUnstructuredGrid() 779 780 elif isinstance(inputobj, vtk.vtkUnstructuredGrid): 781 self.dataset = inputobj 782 783 elif isinstance(inputobj, UnstructuredGrid): 784 self.dataset = inputobj.dataset 785 786 elif isinstance(inputobj, vtk.vtkRectilinearGrid): 787 r2t = vtk.new("RectilinearGridToTetrahedra") 788 r2t.SetInputData(inputobj) 789 r2t.RememberVoxelIdOn() 790 r2t.SetTetraPerCellTo6() 791 r2t.Update() 792 self.dataset = r2t.GetOutput() 793 794 elif isinstance(inputobj, vtk.vtkDataSet): 795 r2t = vtk.new("DataSetTriangleFilter") 796 r2t.SetInputData(inputobj) 797 r2t.TetrahedraOnlyOn() 798 r2t.Update() 799 self.dataset = r2t.GetOutput() 800 801 elif isinstance(inputobj, str): 802 if "https://" in inputobj: 803 inputobj = download(inputobj, verbose=False) 804 if inputobj.endswith(".vtu"): 805 reader = vtk.new("XMLUnstructuredGridReader") 806 else: 807 reader = vtk.new("UnstructuredGridReader") 808 self.filename = inputobj 809 reader.SetFileName(inputobj) 810 reader.Update() 811 ug = reader.GetOutput() 812 tt = vtk.new("DataSetTriangleFilter") 813 tt.SetInputData(ug) 814 tt.SetTetrahedraOnly(True) 815 tt.Update() 816 self.dataset = tt.GetOutput() 817 818 elif utils.is_sequence(inputobj): 819 self.dataset = vtk.vtkUnstructuredGrid() 820 821 points, cells = inputobj 822 if len(points) == 0: 823 return 824 if not utils.is_sequence(points[0]): 825 return 826 if len(cells) == 0: 827 return 828 829 if not utils.is_sequence(cells[0]): 830 tets = [] 831 nf = cells[0] + 1 832 for i, cl in enumerate(cells): 833 if i in (nf, 0): 834 k = i + 1 835 nf = cl + k 836 cell = [cells[j + k] for j in range(cl)] 837 tets.append(cell) 838 cells = tets 839 840 source_points = vtk.vtkPoints() 841 varr = utils.numpy2vtk(points, dtype=np.float32) 842 source_points.SetData(varr) 843 self.dataset.SetPoints(source_points) 844 845 source_tets = vtk.vtkCellArray() 846 for f in cells: 847 ele = vtk.vtkTetra() 848 pid = ele.GetPointIds() 849 for i, fi in enumerate(f): 850 pid.SetId(i, fi) 851 source_tets.InsertNextCell(ele) 852 self.dataset.SetCells(vtk.VTK_TETRA, source_tets) 853 854 else: 855 vedo.logger.error(f"cannot understand input type {inputtype}") 856 return 857 858 self.pipeline = utils.OperationNode( 859 self, comment=f"#tets {self.dataset.GetNumberOfCells()}", 860 c="#9e2a2b", 861 )
Arguments:
- inputobj : (vtkDataSet, list, str) list of points and tet indices, or filename
1001 def compute_quality(self, metric=7): 1002 """ 1003 Calculate functions of quality for the elements of a triangular mesh. 1004 This method adds to the mesh a cell array named "Quality". 1005 1006 Arguments: 1007 metric : (int) 1008 type of estimators: 1009 - EDGE RATIO, 0 1010 - ASPECT RATIO, 1 1011 - RADIUS RATIO, 2 1012 - ASPECT FROBENIUS, 3 1013 - MIN_ANGLE, 4 1014 - COLLAPSE RATIO, 5 1015 - ASPECT GAMMA, 6 1016 - VOLUME, 7 1017 - ... 1018 1019 See class [vtkMeshQuality](https://vtk.org/doc/nightly/html/classvtkMeshQuality.html) 1020 for an explanation of the meaning of each metric.. 1021 """ 1022 qf = vtk.new("MeshQuality") 1023 qf.SetInputData(self.dataset) 1024 qf.SetTetQualityMeasure(metric) 1025 qf.SaveCellQualityOn() 1026 qf.Update() 1027 self._update(qf.GetOutput()) 1028 return utils.vtk2numpy(qf.GetOutput().GetCellData().GetArray("Quality"))
Calculate functions of quality for the elements of a triangular mesh. This method adds to the mesh a cell array named "Quality".
Arguments:
- metric : (int) type of estimators: - EDGE RATIO, 0 - ASPECT RATIO, 1 - RADIUS RATIO, 2 - ASPECT FROBENIUS, 3 - MIN_ANGLE, 4 - COLLAPSE RATIO, 5 - ASPECT GAMMA, 6 - VOLUME, 7 - ...
See class vtkMeshQuality for an explanation of the meaning of each metric..
1030 def check_validity(self, tol=0): 1031 """ 1032 Return an array of possible problematic tets following this convention: 1033 ```python 1034 Valid = 0 1035 WrongNumberOfPoints = 01 1036 IntersectingEdges = 02 1037 IntersectingFaces = 04 1038 NoncontiguousEdges = 08 1039 Nonconvex = 10 1040 OrientedIncorrectly = 20 1041 ``` 1042 1043 Arguments: 1044 tol : (float) 1045 This value is used as an epsilon for floating point 1046 equality checks throughout the cell checking process. 1047 """ 1048 vald = vtk.new("CellValidator") 1049 if tol: 1050 vald.SetTolerance(tol) 1051 vald.SetInputData(self.dataset) 1052 vald.Update() 1053 varr = vald.GetOutput().GetCellData().GetArray("ValidityState") 1054 return utils.vtk2numpy(varr)
Return an array of possible problematic tets following this convention:
Valid = 0
WrongNumberOfPoints = 01
IntersectingEdges = 02
IntersectingFaces = 04
NoncontiguousEdges = 08
Nonconvex = 10
OrientedIncorrectly = 20
Arguments:
- tol : (float) This value is used as an epsilon for floating point equality checks throughout the cell checking process.
1056 def decimate(self, scalars_name, fraction=0.5, n=0): 1057 """ 1058 Downsample the number of tets in a TetMesh to a specified fraction. 1059 Either `fraction` or `n` must be set. 1060 1061 Arguments: 1062 fraction : (float) 1063 the desired final fraction of the total. 1064 n : (int) 1065 the desired number of final tets 1066 1067 .. note:: setting `fraction=0.1` leaves 10% of the original nr of tets. 1068 """ 1069 decimate = vtk.new("UnstructuredGridQuadricDecimation") 1070 decimate.SetInputData(self.dataset) 1071 decimate.SetScalarsName(scalars_name) 1072 1073 if n: # n = desired number of points 1074 decimate.SetNumberOfTetsOutput(n) 1075 else: 1076 decimate.SetTargetReduction(1 - fraction) 1077 decimate.Update() 1078 self._update(decimate.GetOutput()) 1079 self.pipeline = utils.OperationNode( 1080 "decimate", comment=f"array: {scalars_name}", 1081 c="#edabab", parents=[self], 1082 ) 1083 return self
Downsample the number of tets in a TetMesh to a specified fraction.
Either fraction
or n
must be set.
Arguments:
- fraction : (float) the desired final fraction of the total.
- n : (int) the desired number of final tets
setting fraction=0.1
leaves 10% of the original nr of tets.
1085 def subdvide(self): 1086 """ 1087 Increase the number of tets of a `TetMesh`. 1088 Subdivide one tetrahedron into twelve for every tetra. 1089 """ 1090 sd = vtk.new("SubdivideTetra") 1091 sd.SetInputData(self.dataset) 1092 sd.Update() 1093 self._update(sd.GetOutput()) 1094 self.pipeline = utils.OperationNode( 1095 "subdvide", c="#edabab", parents=[self], 1096 ) 1097 return self
Increase the number of tets of a TetMesh
.
Subdivide one tetrahedron into twelve for every tetra.
1099 def isosurface(self, value=True): 1100 """ 1101 Return a `vedo.Mesh` isosurface. 1102 1103 Set `value` to a single value or list of values to compute the isosurface(s). 1104 """ 1105 if not self.dataset.GetPointData().GetScalars(): 1106 self.map_cells_to_points() 1107 scrange = self.dataset.GetPointData().GetScalars().GetRange() 1108 cf = vtk.new("ContourFilter") # vtk.new("ContourGrid") 1109 cf.SetInputData(self.dataset) 1110 1111 if utils.is_sequence(value): 1112 cf.SetNumberOfContours(len(value)) 1113 for i, t in enumerate(value): 1114 cf.SetValue(i, t) 1115 cf.Update() 1116 else: 1117 if value is True: 1118 value = (2 * scrange[0] + scrange[1]) / 3.0 1119 cf.SetValue(0, value) 1120 cf.Update() 1121 1122 clp = vtk.new("CleanPolyData") 1123 clp.SetInputData(cf.GetOutput()) 1124 clp.Update() 1125 msh = Mesh(clp.GetOutput(), c=None).phong() 1126 msh.copy_properties_from(self) 1127 msh.pipeline = utils.OperationNode("isosurface", c="#edabab", parents=[self]) 1128 return msh
Return a vedo.Mesh
isosurface.
Set value
to a single value or list of values to compute the isosurface(s).
1130 def slice(self, origin=(0, 0, 0), normal=(1, 0, 0)): 1131 """ 1132 Return a 2D slice of the mesh by a plane passing through origin and assigned normal. 1133 """ 1134 strn = str(normal) 1135 if strn == "x": normal = (1, 0, 0) 1136 elif strn == "y": normal = (0, 1, 0) 1137 elif strn == "z": normal = (0, 0, 1) 1138 elif strn == "-x": normal = (-1, 0, 0) 1139 elif strn == "-y": normal = (0, -1, 0) 1140 elif strn == "-z": normal = (0, 0, -1) 1141 plane = vtk.new("Plane") 1142 plane.SetOrigin(origin) 1143 plane.SetNormal(normal) 1144 1145 cc = vtk.new("Cutter") 1146 cc.SetInputData(self.dataset) 1147 cc.SetCutFunction(plane) 1148 cc.Update() 1149 msh = Mesh(cc.GetOutput()).flat().lighting("ambient") 1150 msh.copy_properties_from(self) 1151 msh.pipeline = utils.OperationNode("slice", c="#edabab", parents=[self]) 1152 return msh
Return a 2D slice of the mesh by a plane passing through origin and assigned normal.
Inherited Members
- UnstructuredGrid
- actor
- copy
- clone
- bounds
- threshold
- tomesh
- extract_cell_by_type
- extract_cells_by_id
- find_cell
- clean
- extract_cells_on_plane
- extract_cells_on_sphere
- extract_cells_on_cylinder
- cut_with_plane
- cut_with_box
- cut_with_mesh
- vedo.core.PointAlgorithms
- apply_transform
- apply_transform_from_actor
- pos
- shift
- x
- y
- z
- rotate
- rotate_x
- rotate_y
- rotate_z
- reorient
- scale
- vedo.core.CommonAlgorithms
- pointdata
- celldata
- metadata
- memory_address
- memory_size
- modified
- box
- xbounds
- ybounds
- zbounds
- diagonal_size
- average_size
- center_of_mass
- copy_data_from
- inputdata
- npoints
- nvertices
- ncells
- points
- cell_centers
- lines
- lines_as_flat_array
- mark_boundaries
- find_cells_in_bounds
- find_cells_along_line
- find_cells_along_plane
- delete_cells_by_point_index
- map_cells_to_points
- vertices
- coordinates
- cells_as_flat_array
- cells
- map_points_to_cells
- resample_data_from
- interpolate_data_from
- add_ids
- gradient
- divergence
- vorticity
- probe
- compute_cell_size
- integrate_arrays_over_domain
- write
- shrink
- vedo.visual.MeshVisual
- follow_camera
- wireframe
- flat
- phong
- backface_culling
- render_lines_as_tubes
- frontface_culling
- backcolor
- bc
- linewidth
- lw
- linecolor
- lc
- texture
- vedo.visual.PointsVisual
- clone2d
- copy_properties_from
- color
- c
- alpha
- lut_color_at
- opacity
- force_opaque
- force_translucent
- point_size
- ps
- render_points_as_spheres
- lighting
- point_blurring
- cellcolors
- pointcolors
- cmap
- add_trail
- update_trail
- add_shadow
- update_shadows
- labels
- labels2d
- legend
- flagpole
- flagpost
- caption