vedo.grids
Work with tetrahedral meshes.
1#!/usr/bin/env python3 2# -*- coding: utf-8 -*- 3import os 4import time 5from weakref import ref as weak_ref_to 6from typing import Any 7from typing_extensions import Self 8import numpy as np 9 10import vedo.vtkclasses as vtki # a wrapper for lazy imports 11 12import vedo 13from vedo import utils 14from vedo.core import PointAlgorithms 15from vedo.mesh import Mesh 16from vedo.file_io import download 17from vedo.visual import MeshVisual 18from vedo.transformations import LinearTransform 19 20__docformat__ = "google" 21 22__doc__ = """ 23Work with tetrahedral meshes. 24 25 26""" 27 28__all__ = [ 29 "UnstructuredGrid", 30 "TetMesh", 31 "RectilinearGrid", 32 "StructuredGrid", 33] 34 35 36######################################################################### 37class UnstructuredGrid(PointAlgorithms, MeshVisual): 38 """Support for UnstructuredGrid objects.""" 39 40 def __init__(self, inputobj=None): 41 """ 42 Support for UnstructuredGrid objects. 43 44 Arguments: 45 inputobj : (list, vtkUnstructuredGrid, str) 46 A list in the form `[points, cells, celltypes]`, 47 or a vtkUnstructuredGrid object, or a filename 48 49 Celltypes are identified by the following 50 [convention](https://vtk.org/doc/nightly/html/vtkCellType_8h_source.html). 51 """ 52 super().__init__() 53 54 self.dataset = None 55 56 self.mapper = vtki.new("PolyDataMapper") 57 self._actor = vtki.vtkActor() 58 self._actor.retrieve_object = weak_ref_to(self) 59 self._actor.SetMapper(self.mapper) 60 self.properties = self._actor.GetProperty() 61 62 self.transform = LinearTransform() 63 self.point_locator = None 64 self.cell_locator = None 65 self.line_locator = None 66 67 self.name = "UnstructuredGrid" 68 self.filename = "" 69 self.file_size = "" 70 71 self.info = {} 72 self.time = time.time() 73 self.rendered_at = set() 74 75 ###################################### 76 inputtype = str(type(inputobj)) 77 78 if inputobj is None: 79 self.dataset = vtki.vtkUnstructuredGrid() 80 81 elif utils.is_sequence(inputobj): 82 83 pts, cells, celltypes = inputobj 84 assert len(cells) == len(celltypes) 85 86 self.dataset = vtki.vtkUnstructuredGrid() 87 88 if not utils.is_sequence(cells[0]): 89 tets = [] 90 nf = cells[0] + 1 91 for i, cl in enumerate(cells): 92 if i in (nf, 0): 93 k = i + 1 94 nf = cl + k 95 cell = [cells[j + k] for j in range(cl)] 96 tets.append(cell) 97 cells = tets 98 99 # This would fill the points and use those to define orientation 100 vpts = utils.numpy2vtk(pts, dtype=np.float32) 101 points = vtki.vtkPoints() 102 points.SetData(vpts) 103 self.dataset.SetPoints(points) 104 105 # Fill cells 106 # https://vtk.org/doc/nightly/html/vtkCellType_8h_source.html 107 for i, ct in enumerate(celltypes): 108 if ct == vtki.cell_types["VERTEX"]: 109 cell = vtki.vtkVertex() 110 elif ct == vtki.cell_types["POLY_VERTEX"]: 111 cell = vtki.vtkPolyVertex() 112 elif ct == vtki.cell_types["TETRA"]: 113 cell = vtki.vtkTetra() 114 elif ct == vtki.cell_types["WEDGE"]: 115 cell = vtki.vtkWedge() 116 elif ct == vtki.cell_types["LINE"]: 117 cell = vtki.vtkLine() 118 elif ct == vtki.cell_types["POLY_LINE"]: 119 cell = vtki.vtkPolyLine() 120 elif ct == vtki.cell_types["TRIANGLE"]: 121 cell = vtki.vtkTriangle() 122 elif ct == vtki.cell_types["TRIANGLE_STRIP"]: 123 cell = vtki.vtkTriangleStrip() 124 elif ct == vtki.cell_types["POLYGON"]: 125 cell = vtki.vtkPolygon() 126 elif ct == vtki.cell_types["PIXEL"]: 127 cell = vtki.vtkPixel() 128 elif ct == vtki.cell_types["QUAD"]: 129 cell = vtki.vtkQuad() 130 elif ct == vtki.cell_types["VOXEL"]: 131 cell = vtki.vtkVoxel() 132 elif ct == vtki.cell_types["PYRAMID"]: 133 cell = vtki.vtkPyramid() 134 elif ct == vtki.cell_types["HEXAHEDRON"]: 135 cell = vtki.vtkHexahedron() 136 elif ct == vtki.cell_types["HEXAGONAL_PRISM"]: 137 cell = vtki.vtkHexagonalPrism() 138 elif ct == vtki.cell_types["PENTAGONAL_PRISM"]: 139 cell = vtki.vtkPentagonalPrism() 140 elif ct == vtki.cell_types["QUADRATIC_TETRA"]: 141 from vtkmodules.vtkCommonDataModel import vtkQuadraticTetra 142 cell = vtkQuadraticTetra() 143 elif ct == vtki.cell_types["QUADRATIC_HEXAHEDRON"]: 144 from vtkmodules.vtkCommonDataModel import vtkQuadraticHexahedron 145 cell = vtkQuadraticHexahedron() 146 elif ct == vtki.cell_types["QUADRATIC_WEDGE"]: 147 from vtkmodules.vtkCommonDataModel import vtkQuadraticWedge 148 cell = vtkQuadraticWedge() 149 elif ct == vtki.cell_types["QUADRATIC_PYRAMID"]: 150 from vtkmodules.vtkCommonDataModel import vtkQuadraticPyramid 151 cell = vtkQuadraticPyramid() 152 elif ct == vtki.cell_types["QUADRATIC_LINEAR_QUAD"]: 153 from vtkmodules.vtkCommonDataModel import vtkQuadraticLinearQuad 154 cell = vtkQuadraticLinearQuad() 155 elif ct == vtki.cell_types["QUADRATIC_LINEAR_WEDGE"]: 156 from vtkmodules.vtkCommonDataModel import vtkQuadraticLinearWedge 157 cell = vtkQuadraticLinearWedge() 158 elif ct == vtki.cell_types["BIQUADRATIC_QUADRATIC_WEDGE"]: 159 from vtkmodules.vtkCommonDataModel import vtkBiQuadraticQuadraticWedge 160 cell = vtkBiQuadraticQuadraticWedge() 161 elif ct == vtki.cell_types["BIQUADRATIC_QUADRATIC_HEXAHEDRON"]: 162 from vtkmodules.vtkCommonDataModel import vtkBiQuadraticQuadraticHexahedron 163 cell = vtkBiQuadraticQuadraticHexahedron() 164 elif ct == vtki.cell_types["BIQUADRATIC_TRIANGLE"]: 165 from vtkmodules.vtkCommonDataModel import vtkBiQuadraticTriangle 166 cell = vtkBiQuadraticTriangle() 167 elif ct == vtki.cell_types["CUBIC_LINE"]: 168 from vtkmodules.vtkCommonDataModel import vtkCubicLine 169 cell = vtkCubicLine() 170 elif ct == vtki.cell_types["CONVEX_POINT_SET"]: 171 from vtkmodules.vtkCommonDataModel import vtkConvexPointSet 172 cell = vtkConvexPointSet() 173 elif ct == vtki.cell_types["POLYHEDRON"]: 174 from vtkmodules.vtkCommonDataModel import vtkPolyhedron 175 cell = vtkPolyhedron() 176 elif ct == vtki.cell_types["HIGHER_ORDER_TRIANGLE"]: 177 from vtkmodules.vtkCommonDataModel import vtkHigherOrderTriangle 178 cell = vtkHigherOrderTriangle() 179 elif ct == vtki.cell_types["HIGHER_ORDER_QUAD"]: 180 from vtkmodules.vtkCommonDataModel import vtkHigherOrderQuadrilateral 181 cell = vtkHigherOrderQuadrilateral() 182 elif ct == vtki.cell_types["HIGHER_ORDER_TETRAHEDRON"]: 183 from vtkmodules.vtkCommonDataModel import vtkHigherOrderTetra 184 cell = vtkHigherOrderTetra() 185 elif ct == vtki.cell_types["HIGHER_ORDER_WEDGE"]: 186 from vtkmodules.vtkCommonDataModel import vtkHigherOrderWedge 187 cell = vtkHigherOrderWedge() 188 elif ct == vtki.cell_types["HIGHER_ORDER_HEXAHEDRON"]: 189 from vtkmodules.vtkCommonDataModel import vtkHigherOrderHexahedron 190 cell = vtkHigherOrderHexahedron() 191 elif ct == vtki.cell_types["LAGRANGE_CURVE"]: 192 from vtkmodules.vtkCommonDataModel import vtkLagrangeCurve 193 cell = vtkLagrangeCurve() 194 elif ct == vtki.cell_types["LAGRANGE_TRIANGLE"]: 195 from vtkmodules.vtkCommonDataModel import vtkLagrangeTriangle 196 cell = vtkLagrangeTriangle() 197 elif ct == vtki.cell_types["LAGRANGE_QUADRILATERAL"]: 198 from vtkmodules.vtkCommonDataModel import vtkLagrangeQuadrilateral 199 cell = vtkLagrangeQuadrilateral() 200 elif ct == vtki.cell_types["LAGRANGE_TETRAHEDRON"]: 201 from vtkmodules.vtkCommonDataModel import vtkLagrangeTetra 202 cell = vtkLagrangeTetra() 203 elif ct == vtki.cell_types["LAGRANGE_HEXAHEDRON"]: 204 from vtkmodules.vtkCommonDataModel import vtkLagrangeHexahedron 205 cell = vtkLagrangeHexahedron() 206 elif ct == vtki.cell_types["LAGRANGE_WEDGE"]: 207 from vtkmodules.vtkCommonDataModel import vtkLagrangeWedge 208 cell = vtkLagrangeWedge() 209 elif ct == vtki.cell_types["BEZIER_CURVE"]: 210 from vtkmodules.vtkCommonDataModel import vtkBezierCurve 211 cell = vtkBezierCurve() 212 elif ct == vtki.cell_types["BEZIER_TRIANGLE"]: 213 from vtkmodules.vtkCommonDataModel import vtkBezierTriangle 214 cell = vtkBezierTriangle() 215 elif ct == vtki.cell_types["BEZIER_QUADRILATERAL"]: 216 from vtkmodules.vtkCommonDataModel import vtkBezierQuadrilateral 217 cell = vtkBezierQuadrilateral() 218 elif ct == vtki.cell_types["BEZIER_TETRAHEDRON"]: 219 from vtkmodules.vtkCommonDataModel import vtkBezierTetra 220 cell = vtkBezierTetra() 221 elif ct == vtki.cell_types["BEZIER_HEXAHEDRON"]: 222 from vtkmodules.vtkCommonDataModel import vtkBezierHexahedron 223 cell = vtkBezierHexahedron() 224 elif ct == vtki.cell_types["BEZIER_WEDGE"]: 225 from vtkmodules.vtkCommonDataModel import vtkBezierWedge 226 cell = vtkBezierWedge() 227 else: 228 vedo.logger.error( 229 f"UnstructuredGrid: cell type {ct} not supported. Skip.") 230 continue 231 232 cpids = cell.GetPointIds() 233 cell_conn = cells[i] 234 for j, pid in enumerate(cell_conn): 235 cpids.SetId(j, pid) 236 self.dataset.InsertNextCell(ct, cpids) 237 238 elif "UnstructuredGrid" in inputtype: 239 self.dataset = inputobj 240 241 elif isinstance(inputobj, str): 242 if "https://" in inputobj: 243 inputobj = download(inputobj, verbose=False) 244 self.filename = inputobj 245 if inputobj.endswith(".vtu"): 246 reader = vtki.new("XMLUnstructuredGridReader") 247 else: 248 reader = vtki.new("UnstructuredGridReader") 249 self.filename = inputobj 250 reader.SetFileName(inputobj) 251 reader.Update() 252 self.dataset = reader.GetOutput() 253 254 else: 255 # this converts other types of vtk objects to UnstructuredGrid 256 apf = vtki.new("AppendFilter") 257 try: 258 apf.AddInputData(inputobj) 259 except TypeError: 260 apf.AddInputData(inputobj.dataset) 261 apf.Update() 262 self.dataset = apf.GetOutput() 263 264 self.properties.SetColor(0.89, 0.455, 0.671) # pink7 265 266 self.pipeline = utils.OperationNode( 267 self, comment=f"#cells {self.dataset.GetNumberOfCells()}", c="#4cc9f0" 268 ) 269 270 # ------------------------------------------------------------------ 271 def __str__(self): 272 """Print a string summary of the `UnstructuredGrid` object.""" 273 module = self.__class__.__module__ 274 name = self.__class__.__name__ 275 out = vedo.printc( 276 f"{module}.{name} at ({hex(self.memory_address())})".ljust(75), 277 c="m", bold=True, invert=True, return_string=True, 278 ) 279 out += "\x1b[0m\u001b[35m" 280 281 out += "nr. of verts".ljust(14) + ": " + str(self.npoints) + "\n" 282 out += "nr. of cells".ljust(14) + ": " + str(self.ncells) + "\n" 283 ct_arr = np.unique(self.cell_types_array) 284 cnames = [k for k, v in vtki.cell_types.items() if v in ct_arr] 285 out += "cell types".ljust(14) + ": " + str(cnames) + "\n" 286 287 if self.npoints: 288 out+="size".ljust(14)+ ": average=" + utils.precision(self.average_size(),6) 289 out+=", diagonal="+ utils.precision(self.diagonal_size(), 6)+ "\n" 290 out+="center of mass".ljust(14) + ": " + utils.precision(self.center_of_mass(),6)+"\n" 291 292 bnds = self.bounds() 293 bx1, bx2 = utils.precision(bnds[0], 3), utils.precision(bnds[1], 3) 294 by1, by2 = utils.precision(bnds[2], 3), utils.precision(bnds[3], 3) 295 bz1, bz2 = utils.precision(bnds[4], 3), utils.precision(bnds[5], 3) 296 out += "bounds".ljust(14) + ":" 297 out += " x=(" + bx1 + ", " + bx2 + ")," 298 out += " y=(" + by1 + ", " + by2 + ")," 299 out += " z=(" + bz1 + ", " + bz2 + ")\n" 300 301 for key in self.pointdata.keys(): 302 arr = self.pointdata[key] 303 rng = utils.precision(arr.min(), 3) + ", " + utils.precision(arr.max(), 3) 304 mark_active = "pointdata" 305 a_scalars = self.dataset.GetPointData().GetScalars() 306 a_vectors = self.dataset.GetPointData().GetVectors() 307 a_tensors = self.dataset.GetPointData().GetTensors() 308 if a_scalars and a_scalars.GetName() == key: 309 mark_active += " *" 310 elif a_vectors and a_vectors.GetName() == key: 311 mark_active += " **" 312 elif a_tensors and a_tensors.GetName() == key: 313 mark_active += " ***" 314 out += mark_active.ljust(14) + f': "{key}" ({arr.dtype}), ndim={arr.ndim}' 315 out += f", range=({rng})\n" 316 317 for key in self.celldata.keys(): 318 arr = self.celldata[key] 319 rng = utils.precision(arr.min(), 3) + ", " + utils.precision(arr.max(), 3) 320 mark_active = "celldata" 321 a_scalars = self.dataset.GetCellData().GetScalars() 322 a_vectors = self.dataset.GetCellData().GetVectors() 323 a_tensors = self.dataset.GetCellData().GetTensors() 324 if a_scalars and a_scalars.GetName() == key: 325 mark_active += " *" 326 elif a_vectors and a_vectors.GetName() == key: 327 mark_active += " **" 328 elif a_tensors and a_tensors.GetName() == key: 329 mark_active += " ***" 330 out += mark_active.ljust(14) + f': "{key}" ({arr.dtype}), ndim={arr.ndim}' 331 out += f", range=({rng})\n" 332 333 for key in self.metadata.keys(): 334 arr = self.metadata[key] 335 out += "metadata".ljust(14) + ": " + f'"{key}" ({len(arr)} values)\n' 336 337 return out.rstrip() + "\x1b[0m" 338 339 def _repr_html_(self): 340 """ 341 HTML representation of the UnstructuredGrid object for Jupyter Notebooks. 342 343 Returns: 344 HTML text with the image and some properties. 345 """ 346 import io 347 import base64 348 from PIL import Image 349 350 library_name = "vedo.grids.UnstructuredGrid" 351 help_url = "https://vedo.embl.es/docs/vedo/grids.html#UnstructuredGrid" 352 353 arr = self.thumbnail() 354 im = Image.fromarray(arr) 355 buffered = io.BytesIO() 356 im.save(buffered, format="PNG", quality=100) 357 encoded = base64.b64encode(buffered.getvalue()).decode("utf-8") 358 url = "data:image/png;base64," + encoded 359 image = f"<img src='{url}'></img>" 360 361 bounds = "<br/>".join( 362 [ 363 utils.precision(min_x, 4) + " ... " + utils.precision(max_x, 4) 364 for min_x, max_x in zip(self.bounds()[::2], self.bounds()[1::2]) 365 ] 366 ) 367 368 help_text = "" 369 if self.name: 370 help_text += f"<b> {self.name}:   </b>" 371 help_text += '<b><a href="' + help_url + '" target="_blank">' + library_name + "</a></b>" 372 if self.filename: 373 dots = "" 374 if len(self.filename) > 30: 375 dots = "..." 376 help_text += f"<br/><code><i>({dots}{self.filename[-30:]})</i></code>" 377 378 pdata = "" 379 if self.dataset.GetPointData().GetScalars(): 380 if self.dataset.GetPointData().GetScalars().GetName(): 381 name = self.dataset.GetPointData().GetScalars().GetName() 382 pdata = "<tr><td><b> point data array </b></td><td>" + name + "</td></tr>" 383 384 cdata = "" 385 if self.dataset.GetCellData().GetScalars(): 386 if self.dataset.GetCellData().GetScalars().GetName(): 387 name = self.dataset.GetCellData().GetScalars().GetName() 388 cdata = "<tr><td><b> cell data array </b></td><td>" + name + "</td></tr>" 389 390 pts = self.vertices 391 cm = np.mean(pts, axis=0) 392 393 all = [ 394 "<table>", 395 "<tr>", 396 "<td>", image, "</td>", 397 "<td style='text-align: center; vertical-align: center;'><br/>", help_text, 398 "<table>", 399 "<tr><td><b> bounds </b> <br/> (x/y/z) </td><td>" + str(bounds) + "</td></tr>", 400 "<tr><td><b> center of mass </b></td><td>" + utils.precision(cm,3) + "</td></tr>", 401 # "<tr><td><b> average size </b></td><td>" + str(average_size) + "</td></tr>", 402 "<tr><td><b> nr. points / cells </b></td><td>" 403 + str(self.npoints) + " / " + str(self.ncells) + "</td></tr>", 404 pdata, 405 cdata, 406 "</table>", 407 "</table>", 408 ] 409 return "\n".join(all) 410 411 @property 412 def actor(self): 413 """Return the `vtkActor` of the object.""" 414 # print("building actor") 415 gf = vtki.new("GeometryFilter") 416 gf.SetInputData(self.dataset) 417 gf.Update() 418 out = gf.GetOutput() 419 self.mapper.SetInputData(out) 420 self.mapper.Modified() 421 return self._actor 422 423 @actor.setter 424 def actor(self, _): 425 pass 426 427 def _update(self, data, reset_locators=False): 428 self.dataset = data 429 if reset_locators: 430 self.cell_locator = None 431 self.point_locator = None 432 return self 433 434 def merge(self, *others) -> Self: 435 """ 436 Merge multiple datasets into one single `UnstrcturedGrid`. 437 """ 438 apf = vtki.new("AppendFilter") 439 for o in others: 440 if isinstance(o, UnstructuredGrid): 441 apf.AddInputData(o.dataset) 442 elif isinstance(o, vtki.vtkUnstructuredGrid): 443 apf.AddInputData(o) 444 else: 445 vedo.printc("Error: cannot merge type", type(o), c="r") 446 apf.Update() 447 self._update(apf.GetOutput()) 448 self.pipeline = utils.OperationNode( 449 "merge", parents=[self, *others], c="#9e2a2b" 450 ) 451 return self 452 453 def copy(self, deep=True) -> "UnstructuredGrid": 454 """Return a copy of the object. Alias of `clone()`.""" 455 return self.clone(deep=deep) 456 457 def clone(self, deep=True) -> "UnstructuredGrid": 458 """Clone the UnstructuredGrid object to yield an exact copy.""" 459 ug = vtki.vtkUnstructuredGrid() 460 if deep: 461 ug.DeepCopy(self.dataset) 462 else: 463 ug.ShallowCopy(self.dataset) 464 if isinstance(self, vedo.UnstructuredGrid): 465 cloned = vedo.UnstructuredGrid(ug) 466 else: 467 cloned = vedo.TetMesh(ug) 468 469 cloned.copy_properties_from(self) 470 471 cloned.pipeline = utils.OperationNode( 472 "clone", parents=[self], shape="diamond", c="#bbe1ed" 473 ) 474 return cloned 475 476 def bounds(self) -> np.ndarray: 477 """ 478 Get the object bounds. 479 Returns a list in format `[xmin,xmax, ymin,ymax, zmin,zmax]`. 480 """ 481 # OVERRIDE CommonAlgorithms.bounds() which is too slow 482 return np.array(self.dataset.GetBounds()) 483 484 def threshold(self, name=None, above=None, below=None, on="cells") -> Self: 485 """ 486 Threshold the tetrahedral mesh by a cell scalar value. 487 Reduce to only tets which satisfy the threshold limits. 488 489 - if `above = below` will only select tets with that specific value. 490 - if `above > below` selection range is flipped. 491 492 Set keyword "on" to either "cells" or "points". 493 """ 494 th = vtki.new("Threshold") 495 th.SetInputData(self.dataset) 496 497 if name is None: 498 if self.celldata.keys(): 499 name = self.celldata.keys()[0] 500 th.SetInputArrayToProcess(0, 0, 0, 1, name) 501 elif self.pointdata.keys(): 502 name = self.pointdata.keys()[0] 503 th.SetInputArrayToProcess(0, 0, 0, 0, name) 504 if name is None: 505 vedo.logger.warning("cannot find active array. Skip.") 506 return self 507 else: 508 if on.startswith("c"): 509 th.SetInputArrayToProcess(0, 0, 0, 1, name) 510 else: 511 th.SetInputArrayToProcess(0, 0, 0, 0, name) 512 513 if above is not None: 514 th.SetLowerThreshold(above) 515 516 if below is not None: 517 th.SetUpperThreshold(below) 518 519 th.Update() 520 return self._update(th.GetOutput()) 521 522 def isosurface(self, value=None, flying_edges=False) -> "vedo.mesh.Mesh": 523 """ 524 Return an `Mesh` isosurface extracted from the `Volume` object. 525 526 Set `value` as single float or list of values to draw the isosurface(s). 527 Use flying_edges for faster results (but sometimes can interfere with `smooth()`). 528 529 Examples: 530 - [isosurfaces1.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/isosurfaces1.py) 531 532  533 """ 534 scrange = self.dataset.GetScalarRange() 535 536 if flying_edges: 537 cf = vtki.new("FlyingEdges3D") 538 cf.InterpolateAttributesOn() 539 else: 540 cf = vtki.new("ContourFilter") 541 cf.UseScalarTreeOn() 542 543 cf.SetInputData(self.dataset) 544 cf.ComputeNormalsOn() 545 546 if utils.is_sequence(value): 547 cf.SetNumberOfContours(len(value)) 548 for i, t in enumerate(value): 549 cf.SetValue(i, t) 550 else: 551 if value is None: 552 value = (2 * scrange[0] + scrange[1]) / 3.0 553 # print("automatic isosurface value =", value) 554 cf.SetValue(0, value) 555 556 cf.Update() 557 poly = cf.GetOutput() 558 559 out = vedo.mesh.Mesh(poly, c=None).flat() 560 out.mapper.SetScalarRange(scrange[0], scrange[1]) 561 562 out.pipeline = utils.OperationNode( 563 "isosurface", 564 parents=[self], 565 comment=f"#pts {out.dataset.GetNumberOfPoints()}", 566 c="#4cc9f0:#e9c46a", 567 ) 568 return out 569 570 def shrink(self, fraction=0.8) -> Self: 571 """ 572 Shrink the individual cells. 573 574  575 """ 576 sf = vtki.new("ShrinkFilter") 577 sf.SetInputData(self.dataset) 578 sf.SetShrinkFactor(fraction) 579 sf.Update() 580 out = sf.GetOutput() 581 self._update(out) 582 self.pipeline = utils.OperationNode( 583 "shrink", comment=f"by {fraction}", parents=[self], c="#9e2a2b" 584 ) 585 return self 586 587 def tomesh(self, fill=False, shrink=1.0) -> "vedo.mesh.Mesh": 588 """ 589 Build a polygonal `Mesh` from the current object. 590 591 If `fill=True`, the interior faces of all the cells are created. 592 (setting a `shrink` value slightly smaller than the default 1.0 593 can avoid flickering due to internal adjacent faces). 594 595 If `fill=False`, only the boundary faces will be generated. 596 """ 597 gf = vtki.new("GeometryFilter") 598 if fill: 599 sf = vtki.new("ShrinkFilter") 600 sf.SetInputData(self.dataset) 601 sf.SetShrinkFactor(shrink) 602 sf.Update() 603 gf.SetInputData(sf.GetOutput()) 604 gf.Update() 605 poly = gf.GetOutput() 606 else: 607 gf.SetInputData(self.dataset) 608 gf.Update() 609 poly = gf.GetOutput() 610 611 msh = vedo.mesh.Mesh(poly) 612 msh.copy_properties_from(self) 613 msh.pipeline = utils.OperationNode( 614 "tomesh", parents=[self], comment=f"fill={fill}", c="#9e2a2b:#e9c46a" 615 ) 616 return msh 617 618 @property 619 def cell_types_array(self): 620 """Return the list of cell types in the dataset.""" 621 uarr = self.dataset.GetCellTypesArray() 622 return utils.vtk2numpy(uarr) 623 624 def extract_cells_by_type(self, ctype) -> "UnstructuredGrid": 625 """Extract a specific cell type and return a new `UnstructuredGrid`.""" 626 if isinstance(ctype, str): 627 try: 628 ctype = vtki.cell_types[ctype.upper()] 629 except KeyError: 630 vedo.logger.error(f"extract_cells_by_type: cell type {ctype} does not exist. Skip.") 631 return self 632 uarr = self.dataset.GetCellTypesArray() 633 ctarrtyp = np.where(utils.vtk2numpy(uarr) == ctype)[0] 634 uarrtyp = utils.numpy2vtk(ctarrtyp, deep=False, dtype="id") 635 selection_node = vtki.new("SelectionNode") 636 selection_node.SetFieldType(vtki.get_class("SelectionNode").CELL) 637 selection_node.SetContentType(vtki.get_class("SelectionNode").INDICES) 638 selection_node.SetSelectionList(uarrtyp) 639 selection = vtki.new("Selection") 640 selection.AddNode(selection_node) 641 es = vtki.new("ExtractSelection") 642 es.SetInputData(0, self.dataset) 643 es.SetInputData(1, selection) 644 es.Update() 645 646 ug = UnstructuredGrid(es.GetOutput()) 647 ug.pipeline = utils.OperationNode( 648 "extract_cell_type", comment=f"type {ctype}", c="#edabab", parents=[self] 649 ) 650 return ug 651 652 def extract_cells_by_id(self, idlist, use_point_ids=False) -> "UnstructuredGrid": 653 """Return a new `UnstructuredGrid` composed of the specified subset of indices.""" 654 selection_node = vtki.new("SelectionNode") 655 if use_point_ids: 656 selection_node.SetFieldType(vtki.get_class("SelectionNode").POINT) 657 contcells = vtki.get_class("SelectionNode").CONTAINING_CELLS() 658 selection_node.GetProperties().Set(contcells, 1) 659 else: 660 selection_node.SetFieldType(vtki.get_class("SelectionNode").CELL) 661 selection_node.SetContentType(vtki.get_class("SelectionNode").INDICES) 662 vidlist = utils.numpy2vtk(idlist, dtype="id") 663 selection_node.SetSelectionList(vidlist) 664 selection = vtki.new("Selection") 665 selection.AddNode(selection_node) 666 es = vtki.new("ExtractSelection") 667 es.SetInputData(0, self) 668 es.SetInputData(1, selection) 669 es.Update() 670 671 ug = UnstructuredGrid(es.GetOutput()) 672 pr = vtki.vtkProperty() 673 pr.DeepCopy(self.properties) 674 ug.actor.SetProperty(pr) 675 ug.properties = pr 676 677 ug.mapper.SetLookupTable(utils.ctf2lut(self)) 678 ug.pipeline = utils.OperationNode( 679 "extract_cells_by_id", 680 parents=[self], 681 comment=f"#cells {self.dataset.GetNumberOfCells()}", 682 c="#9e2a2b", 683 ) 684 return ug 685 686 def find_cell(self, p: list) -> int: 687 """Locate the cell that contains a point and return the cell ID.""" 688 cell = vtki.vtkTetra() 689 cell_id = vtki.mutable(0) 690 tol2 = vtki.mutable(0) 691 sub_id = vtki.mutable(0) 692 pcoords = [0, 0, 0] 693 weights = [0, 0, 0] 694 cid = self.dataset.FindCell(p, cell, cell_id, tol2, sub_id, pcoords, weights) 695 return cid 696 697 def clean(self) -> Self: 698 """ 699 Cleanup unused points and empty cells 700 """ 701 cl = vtki.new("StaticCleanUnstructuredGrid") 702 cl.SetInputData(self.dataset) 703 cl.RemoveUnusedPointsOn() 704 cl.ProduceMergeMapOff() 705 cl.AveragePointDataOff() 706 cl.Update() 707 708 self._update(cl.GetOutput()) 709 self.pipeline = utils.OperationNode( 710 "clean", 711 parents=[self], 712 comment=f"#cells {self.dataset.GetNumberOfCells()}", 713 c="#9e2a2b", 714 ) 715 return self 716 717 def extract_cells_on_plane(self, origin: tuple, normal: tuple) -> Self: 718 """ 719 Extract cells that are lying of the specified surface. 720 """ 721 bf = vtki.new("3DLinearGridCrinkleExtractor") 722 bf.SetInputData(self.dataset) 723 bf.CopyPointDataOn() 724 bf.CopyCellDataOn() 725 bf.RemoveUnusedPointsOff() 726 727 plane = vtki.new("Plane") 728 plane.SetOrigin(origin) 729 plane.SetNormal(normal) 730 bf.SetImplicitFunction(plane) 731 bf.Update() 732 733 self._update(bf.GetOutput(), reset_locators=False) 734 self.pipeline = utils.OperationNode( 735 "extract_cells_on_plane", 736 parents=[self], 737 comment=f"#cells {self.dataset.GetNumberOfCells()}", 738 c="#9e2a2b", 739 ) 740 return self 741 742 def extract_cells_on_sphere(self, center: tuple, radius: tuple) -> Self: 743 """ 744 Extract cells that are lying of the specified surface. 745 """ 746 bf = vtki.new("3DLinearGridCrinkleExtractor") 747 bf.SetInputData(self.dataset) 748 bf.CopyPointDataOn() 749 bf.CopyCellDataOn() 750 bf.RemoveUnusedPointsOff() 751 752 sph = vtki.new("Sphere") 753 sph.SetRadius(radius) 754 sph.SetCenter(center) 755 bf.SetImplicitFunction(sph) 756 bf.Update() 757 758 self._update(bf.GetOutput()) 759 self.pipeline = utils.OperationNode( 760 "extract_cells_on_sphere", 761 parents=[self], 762 comment=f"#cells {self.dataset.GetNumberOfCells()}", 763 c="#9e2a2b", 764 ) 765 return self 766 767 def extract_cells_on_cylinder(self, center: tuple, axis: tuple, radius: float) -> Self: 768 """ 769 Extract cells that are lying of the specified surface. 770 """ 771 bf = vtki.new("3DLinearGridCrinkleExtractor") 772 bf.SetInputData(self.dataset) 773 bf.CopyPointDataOn() 774 bf.CopyCellDataOn() 775 bf.RemoveUnusedPointsOff() 776 777 cyl = vtki.new("Cylinder") 778 cyl.SetRadius(radius) 779 cyl.SetCenter(center) 780 cyl.SetAxis(axis) 781 bf.SetImplicitFunction(cyl) 782 bf.Update() 783 784 self.pipeline = utils.OperationNode( 785 "extract_cells_on_cylinder", 786 parents=[self], 787 comment=f"#cells {self.dataset.GetNumberOfCells()}", 788 c="#9e2a2b", 789 ) 790 self._update(bf.GetOutput()) 791 return self 792 793 def cut_with_plane(self, origin=(0, 0, 0), normal="x") -> "UnstructuredGrid": 794 """ 795 Cut the object with the plane defined by a point and a normal. 796 797 Arguments: 798 origin : (list) 799 the cutting plane goes through this point 800 normal : (list, str) 801 normal vector to the cutting plane 802 """ 803 # if isinstance(self, vedo.Volume): 804 # raise RuntimeError("cut_with_plane() is not applicable to Volume objects.") 805 806 strn = str(normal) 807 if strn == "x": normal = (1, 0, 0) 808 elif strn == "y": normal = (0, 1, 0) 809 elif strn == "z": normal = (0, 0, 1) 810 elif strn == "-x": normal = (-1, 0, 0) 811 elif strn == "-y": normal = (0, -1, 0) 812 elif strn == "-z": normal = (0, 0, -1) 813 plane = vtki.new("Plane") 814 plane.SetOrigin(origin) 815 plane.SetNormal(normal) 816 clipper = vtki.new("ClipDataSet") 817 clipper.SetInputData(self.dataset) 818 clipper.SetClipFunction(plane) 819 clipper.GenerateClipScalarsOff() 820 clipper.GenerateClippedOutputOff() 821 clipper.SetValue(0) 822 clipper.Update() 823 cout = clipper.GetOutput() 824 825 if isinstance(cout, vtki.vtkUnstructuredGrid): 826 ug = vedo.UnstructuredGrid(cout) 827 if isinstance(self, vedo.UnstructuredGrid): 828 self._update(cout) 829 self.pipeline = utils.OperationNode("cut_with_plane", parents=[self], c="#9e2a2b") 830 return self 831 ug.pipeline = utils.OperationNode("cut_with_plane", parents=[self], c="#9e2a2b") 832 return ug 833 834 else: 835 self._update(cout) 836 self.pipeline = utils.OperationNode("cut_with_plane", parents=[self], c="#9e2a2b") 837 return self 838 839 def cut_with_box(self, box: Any) -> "UnstructuredGrid": 840 """ 841 Cut the grid with the specified bounding box. 842 843 Parameter box has format [xmin, xmax, ymin, ymax, zmin, zmax]. 844 If an object is passed, its bounding box are used. 845 846 This method always returns a TetMesh object. 847 848 Example: 849 ```python 850 from vedo import * 851 tmesh = TetMesh(dataurl+'limb_ugrid.vtk') 852 tmesh.color('rainbow') 853 cu = Cube(side=500).x(500) # any Mesh works 854 tmesh.cut_with_box(cu).show(axes=1) 855 ``` 856 857  858 """ 859 bc = vtki.new("BoxClipDataSet") 860 bc.SetInputData(self.dataset) 861 try: 862 boxb = box.bounds() 863 except AttributeError: 864 boxb = box 865 866 bc.SetBoxClip(*boxb) 867 bc.Update() 868 cout = bc.GetOutput() 869 870 # output of vtkBoxClipDataSet is always tetrahedrons 871 tm = vedo.TetMesh(cout) 872 tm.pipeline = utils.OperationNode("cut_with_box", parents=[self], c="#9e2a2b") 873 return tm 874 875 def cut_with_mesh(self, mesh: Mesh, invert=False, whole_cells=False, on_boundary=False) -> "UnstructuredGrid": 876 """ 877 Cut a `UnstructuredGrid` or `TetMesh` with a `Mesh`. 878 879 Use `invert` to return cut off part of the input object. 880 """ 881 ug = self.dataset 882 883 ippd = vtki.new("ImplicitPolyDataDistance") 884 ippd.SetInput(mesh.dataset) 885 886 if whole_cells or on_boundary: 887 clipper = vtki.new("ExtractGeometry") 888 clipper.SetInputData(ug) 889 clipper.SetImplicitFunction(ippd) 890 clipper.SetExtractInside(not invert) 891 clipper.SetExtractBoundaryCells(False) 892 if on_boundary: 893 clipper.SetExtractBoundaryCells(True) 894 clipper.SetExtractOnlyBoundaryCells(True) 895 else: 896 signed_dists = vtki.vtkFloatArray() 897 signed_dists.SetNumberOfComponents(1) 898 signed_dists.SetName("SignedDistance") 899 for pointId in range(ug.GetNumberOfPoints()): 900 p = ug.GetPoint(pointId) 901 signed_dist = ippd.EvaluateFunction(p) 902 signed_dists.InsertNextValue(signed_dist) 903 ug.GetPointData().AddArray(signed_dists) 904 ug.GetPointData().SetActiveScalars("SignedDistance") # NEEDED 905 clipper = vtki.new("ClipDataSet") 906 clipper.SetInputData(ug) 907 clipper.SetInsideOut(not invert) 908 clipper.SetValue(0.0) 909 910 clipper.Update() 911 912 out = vedo.UnstructuredGrid(clipper.GetOutput()) 913 out.pipeline = utils.OperationNode("cut_with_mesh", parents=[self], c="#9e2a2b") 914 return out 915 916 917########################################################################## 918class TetMesh(UnstructuredGrid): 919 """The class describing tetrahedral meshes.""" 920 921 def __init__(self, inputobj=None): 922 """ 923 Arguments: 924 inputobj : (vtkUnstructuredGrid, list, str, tetgenpy.TetgenIO) 925 list of points and tet indices, or filename 926 """ 927 super().__init__() 928 929 self.dataset = None 930 931 self.mapper = vtki.new("PolyDataMapper") 932 self._actor = vtki.vtkActor() 933 self._actor.retrieve_object = weak_ref_to(self) 934 self._actor.SetMapper(self.mapper) 935 self.properties = self._actor.GetProperty() 936 937 self.name = "TetMesh" 938 939 # print('TetMesh inputtype', type(inputobj)) 940 941 ################### 942 if inputobj is None: 943 self.dataset = vtki.vtkUnstructuredGrid() 944 945 elif isinstance(inputobj, vtki.vtkUnstructuredGrid): 946 self.dataset = inputobj 947 948 elif isinstance(inputobj, UnstructuredGrid): 949 self.dataset = inputobj.dataset 950 951 elif "TetgenIO" in str(type(inputobj)): # tetgenpy object 952 inputobj = [inputobj.points(), inputobj.tetrahedra()] 953 954 elif isinstance(inputobj, vtki.vtkRectilinearGrid): 955 r2t = vtki.new("RectilinearGridToTetrahedra") 956 r2t.SetInputData(inputobj) 957 r2t.RememberVoxelIdOn() 958 r2t.SetTetraPerCellTo6() 959 r2t.Update() 960 self.dataset = r2t.GetOutput() 961 962 elif isinstance(inputobj, vtki.vtkDataSet): 963 r2t = vtki.new("DataSetTriangleFilter") 964 r2t.SetInputData(inputobj) 965 r2t.TetrahedraOnlyOn() 966 r2t.Update() 967 self.dataset = r2t.GetOutput() 968 969 elif isinstance(inputobj, str): 970 if "https://" in inputobj: 971 inputobj = download(inputobj, verbose=False) 972 if inputobj.endswith(".vtu"): 973 reader = vtki.new("XMLUnstructuredGridReader") 974 else: 975 reader = vtki.new("UnstructuredGridReader") 976 977 if not os.path.isfile(inputobj): 978 # for some reason vtk Reader does not complain 979 vedo.logger.error(f"file {inputobj} not found") 980 raise FileNotFoundError 981 982 self.filename = inputobj 983 reader.SetFileName(inputobj) 984 reader.Update() 985 ug = reader.GetOutput() 986 987 tt = vtki.new("DataSetTriangleFilter") 988 tt.SetInputData(ug) 989 tt.SetTetrahedraOnly(True) 990 tt.Update() 991 self.dataset = tt.GetOutput() 992 993 ############################### 994 if utils.is_sequence(inputobj): 995 self.dataset = vtki.vtkUnstructuredGrid() 996 997 points, cells = inputobj 998 if len(points) == 0: 999 return 1000 if not utils.is_sequence(points[0]): 1001 return 1002 if len(cells) == 0: 1003 return 1004 1005 if not utils.is_sequence(cells[0]): 1006 tets = [] 1007 nf = cells[0] + 1 1008 for i, cl in enumerate(cells): 1009 if i in (nf, 0): 1010 k = i + 1 1011 nf = cl + k 1012 cell = [cells[j + k] for j in range(cl)] 1013 tets.append(cell) 1014 cells = tets 1015 1016 source_points = vtki.vtkPoints() 1017 varr = utils.numpy2vtk(points, dtype=np.float32) 1018 source_points.SetData(varr) 1019 self.dataset.SetPoints(source_points) 1020 1021 source_tets = vtki.vtkCellArray() 1022 for f in cells: 1023 ele = vtki.vtkTetra() 1024 pid = ele.GetPointIds() 1025 for i, fi in enumerate(f): 1026 pid.SetId(i, fi) 1027 source_tets.InsertNextCell(ele) 1028 self.dataset.SetCells(vtki.cell_types["TETRA"], source_tets) 1029 1030 if not self.dataset: 1031 vedo.logger.error(f"cannot understand input type {type(inputobj)}") 1032 return 1033 1034 self.properties.SetColor(0.352, 0.612, 0.996) # blue7 1035 1036 self.pipeline = utils.OperationNode( 1037 self, comment=f"#tets {self.dataset.GetNumberOfCells()}", c="#9e2a2b" 1038 ) 1039 1040 ################################################################## 1041 def __str__(self): 1042 """Print a string summary of the `TetMesh` object.""" 1043 module = self.__class__.__module__ 1044 name = self.__class__.__name__ 1045 out = vedo.printc( 1046 f"{module}.{name} at ({hex(self.memory_address())})".ljust(75), 1047 c="c", bold=True, invert=True, return_string=True, 1048 ) 1049 out += "\x1b[0m\u001b[36m" 1050 1051 out += "nr. of verts".ljust(14) + ": " + str(self.npoints) + "\n" 1052 out += "nr. of tetras".ljust(14) + ": " + str(self.ncells) + "\n" 1053 1054 if self.npoints: 1055 out+="size".ljust(14)+ ": average=" + utils.precision(self.average_size(),6) 1056 out+=", diagonal="+ utils.precision(self.diagonal_size(), 6)+ "\n" 1057 out+="center of mass".ljust(14) + ": " + utils.precision(self.center_of_mass(),6)+"\n" 1058 1059 bnds = self.bounds() 1060 bx1, bx2 = utils.precision(bnds[0], 3), utils.precision(bnds[1], 3) 1061 by1, by2 = utils.precision(bnds[2], 3), utils.precision(bnds[3], 3) 1062 bz1, bz2 = utils.precision(bnds[4], 3), utils.precision(bnds[5], 3) 1063 out += "bounds".ljust(14) + ":" 1064 out += " x=(" + bx1 + ", " + bx2 + ")," 1065 out += " y=(" + by1 + ", " + by2 + ")," 1066 out += " z=(" + bz1 + ", " + bz2 + ")\n" 1067 1068 for key in self.pointdata.keys(): 1069 arr = self.pointdata[key] 1070 rng = utils.precision(arr.min(), 3) + ", " + utils.precision(arr.max(), 3) 1071 mark_active = "pointdata" 1072 a_scalars = self.dataset.GetPointData().GetScalars() 1073 a_vectors = self.dataset.GetPointData().GetVectors() 1074 a_tensors = self.dataset.GetPointData().GetTensors() 1075 if a_scalars and a_scalars.GetName() == key: 1076 mark_active += " *" 1077 elif a_vectors and a_vectors.GetName() == key: 1078 mark_active += " **" 1079 elif a_tensors and a_tensors.GetName() == key: 1080 mark_active += " ***" 1081 out += mark_active.ljust(14) + f': "{key}" ({arr.dtype}), ndim={arr.ndim}' 1082 out += f", range=({rng})\n" 1083 1084 for key in self.celldata.keys(): 1085 arr = self.celldata[key] 1086 rng = utils.precision(arr.min(), 3) + ", " + utils.precision(arr.max(), 3) 1087 mark_active = "celldata" 1088 a_scalars = self.dataset.GetCellData().GetScalars() 1089 a_vectors = self.dataset.GetCellData().GetVectors() 1090 a_tensors = self.dataset.GetCellData().GetTensors() 1091 if a_scalars and a_scalars.GetName() == key: 1092 mark_active += " *" 1093 elif a_vectors and a_vectors.GetName() == key: 1094 mark_active += " **" 1095 elif a_tensors and a_tensors.GetName() == key: 1096 mark_active += " ***" 1097 out += mark_active.ljust(14) + f': "{key}" ({arr.dtype}), ndim={arr.ndim}' 1098 out += f", range=({rng})\n" 1099 1100 for key in self.metadata.keys(): 1101 arr = self.metadata[key] 1102 out += "metadata".ljust(14) + ": " + f'"{key}" ({len(arr)} values)\n' 1103 1104 return out.rstrip() + "\x1b[0m" 1105 1106 def _repr_html_(self): 1107 """ 1108 HTML representation of the TetMesh object for Jupyter Notebooks. 1109 1110 Returns: 1111 HTML text with the image and some properties. 1112 """ 1113 import io 1114 import base64 1115 from PIL import Image 1116 1117 library_name = "vedo.grids.TetMesh" 1118 help_url = "https://vedo.embl.es/docs/vedo/grids.html#TetMesh" 1119 1120 arr = self.thumbnail() 1121 im = Image.fromarray(arr) 1122 buffered = io.BytesIO() 1123 im.save(buffered, format="PNG", quality=100) 1124 encoded = base64.b64encode(buffered.getvalue()).decode("utf-8") 1125 url = "data:image/png;base64," + encoded 1126 image = f"<img src='{url}'></img>" 1127 1128 bounds = "<br/>".join( 1129 [ 1130 utils.precision(min_x, 4) + " ... " + utils.precision(max_x, 4) 1131 for min_x, max_x in zip(self.bounds()[::2], self.bounds()[1::2]) 1132 ] 1133 ) 1134 1135 help_text = "" 1136 if self.name: 1137 help_text += f"<b> {self.name}:   </b>" 1138 help_text += '<b><a href="' + help_url + '" target="_blank">' + library_name + "</a></b>" 1139 if self.filename: 1140 dots = "" 1141 if len(self.filename) > 30: 1142 dots = "..." 1143 help_text += f"<br/><code><i>({dots}{self.filename[-30:]})</i></code>" 1144 1145 pdata = "" 1146 if self.dataset.GetPointData().GetScalars(): 1147 if self.dataset.GetPointData().GetScalars().GetName(): 1148 name = self.dataset.GetPointData().GetScalars().GetName() 1149 pdata = "<tr><td><b> point data array </b></td><td>" + name + "</td></tr>" 1150 1151 cdata = "" 1152 if self.dataset.GetCellData().GetScalars(): 1153 if self.dataset.GetCellData().GetScalars().GetName(): 1154 name = self.dataset.GetCellData().GetScalars().GetName() 1155 cdata = "<tr><td><b> cell data array </b></td><td>" + name + "</td></tr>" 1156 1157 pts = self.vertices 1158 cm = np.mean(pts, axis=0) 1159 1160 allt = [ 1161 "<table>", 1162 "<tr>", 1163 "<td>", image, "</td>", 1164 "<td style='text-align: center; vertical-align: center;'><br/>", help_text, 1165 "<table>", 1166 "<tr><td><b> bounds </b> <br/> (x/y/z) </td><td>" + str(bounds) + "</td></tr>", 1167 "<tr><td><b> center of mass </b></td><td>" + utils.precision(cm,3) + "</td></tr>", 1168 "<tr><td><b> nr. points / tets </b></td><td>" 1169 + str(self.npoints) + " / " + str(self.ncells) + "</td></tr>", 1170 pdata, 1171 cdata, 1172 "</table>", 1173 "</table>", 1174 ] 1175 return "\n".join(allt) 1176 1177 def compute_quality(self, metric=7) -> np.ndarray: 1178 """ 1179 Calculate functions of quality for the elements of a tetrahedral mesh. 1180 This method adds to the mesh a cell array named "Quality". 1181 1182 Arguments: 1183 metric : (int) 1184 type of estimators: 1185 - EDGE RATIO, 0 1186 - ASPECT RATIO, 1 1187 - RADIUS RATIO, 2 1188 - ASPECT FROBENIUS, 3 1189 - MIN_ANGLE, 4 1190 - COLLAPSE RATIO, 5 1191 - ASPECT GAMMA, 6 1192 - VOLUME, 7 1193 - ... 1194 1195 See class [vtkMeshQuality](https://vtk.org/doc/nightly/html/classvtkMeshQuality.html) 1196 for an explanation of the meaning of each metric.. 1197 """ 1198 qf = vtki.new("MeshQuality") 1199 qf.SetInputData(self.dataset) 1200 qf.SetTetQualityMeasure(metric) 1201 qf.SaveCellQualityOn() 1202 qf.Update() 1203 self._update(qf.GetOutput()) 1204 return utils.vtk2numpy(qf.GetOutput().GetCellData().GetArray("Quality")) 1205 1206 def check_validity(self, tol=0) -> np.ndarray: 1207 """ 1208 Return an array of possible problematic tets following this convention: 1209 ```python 1210 Valid = 0 1211 WrongNumberOfPoints = 01 1212 IntersectingEdges = 02 1213 IntersectingFaces = 04 1214 NoncontiguousEdges = 08 1215 Nonconvex = 10 1216 OrientedIncorrectly = 20 1217 ``` 1218 1219 Arguments: 1220 tol : (float) 1221 This value is used as an epsilon for floating point 1222 equality checks throughout the cell checking process. 1223 """ 1224 vald = vtki.new("CellValidator") 1225 if tol: 1226 vald.SetTolerance(tol) 1227 vald.SetInputData(self.dataset) 1228 vald.Update() 1229 varr = vald.GetOutput().GetCellData().GetArray("ValidityState") 1230 return utils.vtk2numpy(varr) 1231 1232 def decimate(self, scalars_name: str, fraction=0.5, n=0) -> Self: 1233 """ 1234 Downsample the number of tets in a TetMesh to a specified fraction. 1235 Either `fraction` or `n` must be set. 1236 1237 Arguments: 1238 fraction : (float) 1239 the desired final fraction of the total. 1240 n : (int) 1241 the desired number of final tets 1242 1243 .. note:: setting `fraction=0.1` leaves 10% of the original nr of tets. 1244 """ 1245 decimate = vtki.new("UnstructuredGridQuadricDecimation") 1246 decimate.SetInputData(self.dataset) 1247 decimate.SetScalarsName(scalars_name) 1248 1249 if n: # n = desired number of points 1250 decimate.SetNumberOfTetsOutput(n) 1251 else: 1252 decimate.SetTargetReduction(1 - fraction) 1253 decimate.Update() 1254 self._update(decimate.GetOutput()) 1255 self.pipeline = utils.OperationNode( 1256 "decimate", comment=f"array: {scalars_name}", c="#edabab", parents=[self] 1257 ) 1258 return self 1259 1260 def subdvide(self) -> Self: 1261 """ 1262 Increase the number of tetrahedrons of a `TetMesh`. 1263 Subdivides each tetrahedron into twelve smaller tetras. 1264 """ 1265 sd = vtki.new("SubdivideTetra") 1266 sd.SetInputData(self.dataset) 1267 sd.Update() 1268 self._update(sd.GetOutput()) 1269 self.pipeline = utils.OperationNode("subdvide", c="#edabab", parents=[self]) 1270 return self 1271 1272 def generate_random_points(self, n, min_radius=0) -> "vedo.Points": 1273 """ 1274 Generate `n` uniformly distributed random points 1275 inside the tetrahedral mesh. 1276 1277 A new point data array is added to the output points 1278 called "OriginalCellID" which contains the index of 1279 the cell ID in which the point was generated. 1280 1281 Arguments: 1282 n : (int) 1283 number of points to generate. 1284 min_radius: (float) 1285 impose a minimum distance between points. 1286 If `min_radius` is set to 0, the points are 1287 generated uniformly at random inside the mesh. 1288 If `min_radius` is set to a positive value, 1289 the points are generated uniformly at random 1290 inside the mesh, but points closer than `min_radius` 1291 to any other point are discarded. 1292 1293 Returns a `vedo.Points` object. 1294 1295 Note: 1296 Consider using `points.probe(msh)` to interpolate 1297 any existing mesh data onto the points. 1298 1299 Example: 1300 ```python 1301 from vedo import * 1302 tmesh = TetMesh(dataurl + "limb.vtu").alpha(0.2) 1303 pts = tmesh.generate_random_points(20000, min_radius=10) 1304 print(pts.pointdata["OriginalCellID"]) 1305 show(pts, tmesh, axes=1).close() 1306 ``` 1307 """ 1308 cmesh = self.compute_cell_size() 1309 tets = cmesh.cells 1310 verts = cmesh.vertices 1311 cumul = np.cumsum(np.abs(cmesh.celldata["Volume"])) 1312 1313 out_pts = [] 1314 orig_cell = [] 1315 for _ in range(n): 1316 random_area = np.random.random() * cumul[-1] 1317 it = np.searchsorted(cumul, random_area) 1318 A, B, C, D = verts[tets[it]] 1319 r1, r2, r3 = sorted(np.random.random(3)) 1320 p = r1 * A + (r2 - r1) * B + (r3 - r2) * C + (1 - r3) * D 1321 out_pts.append(p) 1322 orig_cell.append(it) 1323 orig_cellnp = np.array(orig_cell, dtype=np.uint32) 1324 1325 vpts = vedo.pointcloud.Points(out_pts) 1326 vpts.pointdata["OriginalCellID"] = orig_cellnp 1327 1328 if min_radius > 0: 1329 vpts.subsample(min_radius, absolute=True) 1330 1331 vpts.point_size(5).color("k1") 1332 vpts.name = "RandomPoints" 1333 vpts.pipeline = utils.OperationNode( 1334 "generate_random_points", c="#edabab", parents=[self]) 1335 return vpts 1336 1337 def isosurface(self, value=True, flying_edges=None) -> "vedo.Mesh": 1338 """ 1339 Return a `vedo.Mesh` isosurface. 1340 The "isosurface" is the surface of the region of points 1341 whose values equal to `value`. 1342 1343 Set `value` to a single value or list of values to compute the isosurface(s). 1344 1345 Note that flying_edges option is not available for `TetMesh`. 1346 """ 1347 if flying_edges is not None: 1348 vedo.logger.warning("flying_edges option is not available for TetMesh.") 1349 1350 if not self.dataset.GetPointData().GetScalars(): 1351 vedo.logger.warning( 1352 "in isosurface() no scalar pointdata found. " 1353 "Mappping cells to points." 1354 ) 1355 self.map_cells_to_points() 1356 scrange = self.dataset.GetPointData().GetScalars().GetRange() 1357 cf = vtki.new("ContourFilter") # vtki.new("ContourGrid") 1358 cf.SetInputData(self.dataset) 1359 1360 if utils.is_sequence(value): 1361 cf.SetNumberOfContours(len(value)) 1362 for i, t in enumerate(value): 1363 cf.SetValue(i, t) 1364 cf.Update() 1365 else: 1366 if value is True: 1367 value = (2 * scrange[0] + scrange[1]) / 3.0 1368 cf.SetValue(0, value) 1369 cf.Update() 1370 1371 msh = Mesh(cf.GetOutput(), c=None) 1372 msh.copy_properties_from(self) 1373 msh.pipeline = utils.OperationNode("isosurface", c="#edabab", parents=[self]) 1374 return msh 1375 1376 def slice(self, origin=(0, 0, 0), normal=(1, 0, 0)) -> "vedo.Mesh": 1377 """ 1378 Return a 2D slice of the mesh by a plane passing through origin and assigned normal. 1379 """ 1380 strn = str(normal) 1381 if strn == "x": normal = (1, 0, 0) 1382 elif strn == "y": normal = (0, 1, 0) 1383 elif strn == "z": normal = (0, 0, 1) 1384 elif strn == "-x": normal = (-1, 0, 0) 1385 elif strn == "-y": normal = (0, -1, 0) 1386 elif strn == "-z": normal = (0, 0, -1) 1387 plane = vtki.new("Plane") 1388 plane.SetOrigin(origin) 1389 plane.SetNormal(normal) 1390 1391 cc = vtki.new("Cutter") 1392 cc.SetInputData(self.dataset) 1393 cc.SetCutFunction(plane) 1394 cc.Update() 1395 msh = Mesh(cc.GetOutput()).flat().lighting("ambient") 1396 msh.copy_properties_from(self) 1397 msh.pipeline = utils.OperationNode("slice", c="#edabab", parents=[self]) 1398 return msh 1399 1400 1401########################################################################## 1402class RectilinearGrid(PointAlgorithms, MeshVisual): 1403 """ 1404 Build a rectilinear grid. 1405 """ 1406 1407 def __init__(self, inputobj=None): 1408 """ 1409 A RectilinearGrid is a dataset where edges are parallel to the coordinate axes. 1410 It can be thought of as a tessellation of a box in 3D space, similar to a `Volume` 1411 except that the cells are not necessarily cubes, but they can have different lengths 1412 along each axis. 1413 This can be useful to describe a volume with variable resolution where one needs 1414 to represent a region with higher detail with respect to another region. 1415 1416 Arguments: 1417 inputobj : (vtkRectilinearGrid, list, str) 1418 list of points and tet indices, or filename 1419 1420 Example: 1421 ```python 1422 from vedo import RectilinearGrid, show 1423 1424 xcoords = 7 + np.sqrt(np.arange(0,2500,25)) 1425 ycoords = np.arange(0, 20) 1426 zcoords = np.arange(0, 20) 1427 1428 rgrid = RectilinearGrid([xcoords, ycoords, zcoords]) 1429 1430 print(rgrid) 1431 print(rgrid.x_coordinates().shape) 1432 print(rgrid.compute_structured_coords([20,10,11])) 1433 1434 msh = rgrid.tomesh().lw(1) 1435 1436 show(msh, axes=1, viewup="z") 1437 ``` 1438 """ 1439 1440 super().__init__() 1441 1442 self.dataset = None 1443 1444 self.mapper = vtki.new("PolyDataMapper") 1445 self._actor = vtki.vtkActor() 1446 self._actor.retrieve_object = weak_ref_to(self) 1447 self._actor.SetMapper(self.mapper) 1448 self.properties = self._actor.GetProperty() 1449 1450 self.transform = LinearTransform() 1451 self.point_locator = None 1452 self.cell_locator = None 1453 self.line_locator = None 1454 1455 self.name = "RectilinearGrid" 1456 self.filename = "" 1457 1458 self.info = {} 1459 self.time = time.time() 1460 1461 ############################### 1462 if inputobj is None: 1463 self.dataset = vtki.vtkRectilinearGrid() 1464 1465 elif isinstance(inputobj, vtki.vtkRectilinearGrid): 1466 self.dataset = inputobj 1467 1468 elif isinstance(inputobj, RectilinearGrid): 1469 self.dataset = inputobj.dataset 1470 1471 elif isinstance(inputobj, str): 1472 if "https://" in inputobj: 1473 inputobj = download(inputobj, verbose=False) 1474 if inputobj.endswith(".vtr"): 1475 reader = vtki.new("XMLRectilinearGridReader") 1476 else: 1477 reader = vtki.new("RectilinearGridReader") 1478 self.filename = inputobj 1479 reader.SetFileName(inputobj) 1480 reader.Update() 1481 self.dataset = reader.GetOutput() 1482 1483 elif utils.is_sequence(inputobj): 1484 self.dataset = vtki.vtkRectilinearGrid() 1485 xcoords, ycoords, zcoords = inputobj 1486 nx, ny, nz = len(xcoords), len(ycoords), len(zcoords) 1487 self.dataset.SetDimensions(nx, ny, nz) 1488 self.dataset.SetXCoordinates(utils.numpy2vtk(xcoords)) 1489 self.dataset.SetYCoordinates(utils.numpy2vtk(ycoords)) 1490 self.dataset.SetZCoordinates(utils.numpy2vtk(zcoords)) 1491 1492 ############################### 1493 1494 if not self.dataset: 1495 vedo.logger.error(f"RectilinearGrid: cannot understand input type {type(inputobj)}") 1496 return 1497 1498 self.properties.SetColor(0.352, 0.612, 0.996) # blue7 1499 1500 self.pipeline = utils.OperationNode( 1501 self, comment=f"#tets {self.dataset.GetNumberOfCells()}", c="#9e2a2b" 1502 ) 1503 1504 @property 1505 def actor(self): 1506 """Return the `vtkActor` of the object.""" 1507 gf = vtki.new("GeometryFilter") 1508 gf.SetInputData(self.dataset) 1509 gf.Update() 1510 self.mapper.SetInputData(gf.GetOutput()) 1511 self.mapper.Modified() 1512 return self._actor 1513 1514 @actor.setter 1515 def actor(self, _): 1516 pass 1517 1518 def _update(self, data, reset_locators=False): 1519 self.dataset = data 1520 if reset_locators: 1521 self.cell_locator = None 1522 self.point_locator = None 1523 return self 1524 1525 ################################################################## 1526 def __str__(self): 1527 """Print a summary for the `RectilinearGrid` object.""" 1528 module = self.__class__.__module__ 1529 name = self.__class__.__name__ 1530 out = vedo.printc( 1531 f"{module}.{name} at ({hex(self.memory_address())})".ljust(75), 1532 c="c", bold=True, invert=True, return_string=True, 1533 ) 1534 out += "\x1b[0m\x1b[36;1m" 1535 1536 out += "name".ljust(14) + ": " + str(self.name) + "\n" 1537 if self.filename: 1538 out += "filename".ljust(14) + ": " + str(self.filename) + "\n" 1539 1540 out += "dimensions".ljust(14) + ": " + str(self.dataset.GetDimensions()) + "\n" 1541 1542 out += "center".ljust(14) + ": " 1543 out += utils.precision(self.dataset.GetCenter(), 6) + "\n" 1544 1545 bnds = self.bounds() 1546 bx1, bx2 = utils.precision(bnds[0], 3), utils.precision(bnds[1], 3) 1547 by1, by2 = utils.precision(bnds[2], 3), utils.precision(bnds[3], 3) 1548 bz1, bz2 = utils.precision(bnds[4], 3), utils.precision(bnds[5], 3) 1549 out += "bounds".ljust(14) + ":" 1550 out += " x=(" + bx1 + ", " + bx2 + ")," 1551 out += " y=(" + by1 + ", " + by2 + ")," 1552 out += " z=(" + bz1 + ", " + bz2 + ")\n" 1553 1554 out += "memory size".ljust(14) + ": " 1555 out += utils.precision(self.dataset.GetActualMemorySize() / 1024, 2) + " MB\n" 1556 1557 for key in self.pointdata.keys(): 1558 arr = self.pointdata[key] 1559 rng = utils.precision(arr.min(), 3) + ", " + utils.precision(arr.max(), 3) 1560 mark_active = "pointdata" 1561 a_scalars = self.dataset.GetPointData().GetScalars() 1562 a_vectors = self.dataset.GetPointData().GetVectors() 1563 a_tensors = self.dataset.GetPointData().GetTensors() 1564 if a_scalars and a_scalars.GetName() == key: 1565 mark_active += " *" 1566 elif a_vectors and a_vectors.GetName() == key: 1567 mark_active += " **" 1568 elif a_tensors and a_tensors.GetName() == key: 1569 mark_active += " ***" 1570 out += mark_active.ljust(14) + f': "{key}" ({arr.dtype}), ndim={arr.ndim}' 1571 out += f", range=({rng})\n" 1572 1573 for key in self.celldata.keys(): 1574 arr = self.celldata[key] 1575 rng = utils.precision(arr.min(), 3) + ", " + utils.precision(arr.max(), 3) 1576 mark_active = "celldata" 1577 a_scalars = self.dataset.GetCellData().GetScalars() 1578 a_vectors = self.dataset.GetCellData().GetVectors() 1579 a_tensors = self.dataset.GetCellData().GetTensors() 1580 if a_scalars and a_scalars.GetName() == key: 1581 mark_active += " *" 1582 elif a_vectors and a_vectors.GetName() == key: 1583 mark_active += " **" 1584 elif a_tensors and a_tensors.GetName() == key: 1585 mark_active += " ***" 1586 out += mark_active.ljust(14) + f': "{key}" ({arr.dtype}), ndim={arr.ndim}' 1587 out += f", range=({rng})\n" 1588 1589 for key in self.metadata.keys(): 1590 arr = self.metadata[key] 1591 out += "metadata".ljust(14) + ": " + f'"{key}" ({len(arr)} values)\n' 1592 1593 return out.rstrip() + "\x1b[0m" 1594 1595 def _repr_html_(self): 1596 """ 1597 HTML representation of the RectilinearGrid object for Jupyter Notebooks. 1598 1599 Returns: 1600 HTML text with the image and some properties. 1601 """ 1602 import io 1603 import base64 1604 from PIL import Image 1605 1606 library_name = "vedo.grids.RectilinearGrid" 1607 help_url = "https://vedo.embl.es/docs/vedo/grids.html#RectilinearGrid" 1608 1609 m = self.tomesh().linewidth(1).lighting("off") 1610 arr= m.thumbnail(zoom=1, elevation=-30, azimuth=-30) 1611 1612 im = Image.fromarray(arr) 1613 buffered = io.BytesIO() 1614 im.save(buffered, format="PNG", quality=100) 1615 encoded = base64.b64encode(buffered.getvalue()).decode("utf-8") 1616 url = "data:image/png;base64," + encoded 1617 image = f"<img src='{url}'></img>" 1618 1619 bounds = "<br/>".join( 1620 [ 1621 utils.precision(min_x, 4) + " ... " + utils.precision(max_x, 4) 1622 for min_x, max_x in zip(self.bounds()[::2], self.bounds()[1::2]) 1623 ] 1624 ) 1625 1626 help_text = "" 1627 if self.name: 1628 help_text += f"<b> {self.name}:   </b>" 1629 help_text += '<b><a href="' + help_url + '" target="_blank">' + library_name + "</a></b>" 1630 if self.filename: 1631 dots = "" 1632 if len(self.filename) > 30: 1633 dots = "..." 1634 help_text += f"<br/><code><i>({dots}{self.filename[-30:]})</i></code>" 1635 1636 pdata = "" 1637 if self.dataset.GetPointData().GetScalars(): 1638 if self.dataset.GetPointData().GetScalars().GetName(): 1639 name = self.dataset.GetPointData().GetScalars().GetName() 1640 pdata = "<tr><td><b> point data array </b></td><td>" + name + "</td></tr>" 1641 1642 cdata = "" 1643 if self.dataset.GetCellData().GetScalars(): 1644 if self.dataset.GetCellData().GetScalars().GetName(): 1645 name = self.dataset.GetCellData().GetScalars().GetName() 1646 cdata = "<tr><td><b> cell data array </b></td><td>" + name + "</td></tr>" 1647 1648 pts = self.vertices 1649 cm = np.mean(pts, axis=0) 1650 1651 all = [ 1652 "<table>", 1653 "<tr>", 1654 "<td>", image, "</td>", 1655 "<td style='text-align: center; vertical-align: center;'><br/>", help_text, 1656 "<table>", 1657 "<tr><td><b> bounds </b> <br/> (x/y/z) </td><td>" + str(bounds) + "</td></tr>", 1658 "<tr><td><b> center of mass </b></td><td>" + utils.precision(cm,3) + "</td></tr>", 1659 "<tr><td><b> nr. points / cells </b></td><td>" 1660 + str(self.npoints) + " / " + str(self.ncells) + "</td></tr>", 1661 pdata, 1662 cdata, 1663 "</table>", 1664 "</table>", 1665 ] 1666 return "\n".join(all) 1667 1668 def dimensions(self) -> np.ndarray: 1669 """Return the number of points in the x, y and z directions.""" 1670 return np.array(self.dataset.GetDimensions()) 1671 1672 def x_coordinates(self) -> np.ndarray: 1673 """Return the x-coordinates of the grid.""" 1674 return utils.vtk2numpy(self.dataset.GetXCoordinates()) 1675 1676 def y_coordinates(self) -> np.ndarray: 1677 """Return the y-coordinates of the grid.""" 1678 return utils.vtk2numpy(self.dataset.GetYCoordinates()) 1679 1680 def z_coordinates(self) -> np.ndarray: 1681 """Return the z-coordinates of the grid.""" 1682 return utils.vtk2numpy(self.dataset.GetZCoordinates()) 1683 1684 def is_point_visible(self, pid: int) -> bool: 1685 """Return True if point `pid` is visible.""" 1686 return self.dataset.IsPointVisible(pid) 1687 1688 def is_cell_visible(self, cid: int) -> bool: 1689 """Return True if cell `cid` is visible.""" 1690 return self.dataset.IsCellVisible(cid) 1691 1692 def has_blank_points(self) -> bool: 1693 """Return True if the grid has blank points.""" 1694 return self.dataset.HasAnyBlankPoints() 1695 1696 def has_blank_cells(self) -> bool: 1697 """Return True if the grid has blank cells.""" 1698 return self.dataset.HasAnyBlankCells() 1699 1700 def compute_structured_coords(self, x: list) -> dict: 1701 """ 1702 Convenience function computes the structured coordinates for a point `x`. 1703 1704 This method returns a dictionary with keys `ijk`, `pcoords` and `inside`. 1705 The cell is specified by the array `ijk`. 1706 and the parametric coordinates in the cell are specified with `pcoords`. 1707 Value of `inside` is False if the point x is outside of the grid. 1708 """ 1709 ijk = [0, 0, 0] 1710 pcoords = [0., 0., 0.] 1711 inout = self.dataset.ComputeStructuredCoordinates(x, ijk, pcoords) 1712 return {"ijk": np.array(ijk), "pcoords": np.array(pcoords), "inside": bool(inout)} 1713 1714 def compute_pointid(self, ijk: int) -> int: 1715 """Given a location in structured coordinates (i-j-k), return the point id.""" 1716 return self.dataset.ComputePointId(ijk) 1717 1718 def compute_cellid(self, ijk: int) -> int: 1719 """Given a location in structured coordinates (i-j-k), return the cell id.""" 1720 return self.dataset.ComputeCellId(ijk) 1721 1722 def find_point(self, x: list) -> int: 1723 """Given a position `x`, return the id of the closest point.""" 1724 return self.dataset.FindPoint(x) 1725 1726 def find_cell(self, x: list) -> dict: 1727 """Given a position `x`, return the id of the closest cell.""" 1728 cell = vtki.vtkHexagonalPrism() 1729 cellid = vtki.mutable(0) 1730 tol2 = 0.001 # vtki.mutable(0) 1731 subid = vtki.mutable(0) 1732 pcoords = [0.0, 0.0, 0.0] 1733 weights = [0.0, 0.0, 0.0] 1734 res = self.dataset.FindCell(x, cell, cellid, tol2, subid, pcoords, weights) 1735 result = {} 1736 result["cellid"] = cellid 1737 result["subid"] = subid 1738 result["pcoords"] = pcoords 1739 result["weights"] = weights 1740 result["status"] = res 1741 return result 1742 1743 def clone(self, deep=True) -> "RectilinearGrid": 1744 """Return a clone copy of the RectilinearGrid. Alias of `copy()`.""" 1745 if deep: 1746 newrg = vtki.vtkRectilinearGrid() 1747 newrg.CopyStructure(self.dataset) 1748 newrg.CopyAttributes(self.dataset) 1749 newvol = RectilinearGrid(newrg) 1750 else: 1751 newvol = RectilinearGrid(self.dataset) 1752 1753 prop = vtki.vtkProperty() 1754 prop.DeepCopy(self.properties) 1755 newvol.actor.SetProperty(prop) 1756 newvol.properties = prop 1757 newvol.pipeline = utils.OperationNode("clone", parents=[self], c="#bbd0ff", shape="diamond") 1758 return newvol 1759 1760 def bounds(self) -> np.ndarray: 1761 """ 1762 Get the object bounds. 1763 Returns a list in format `[xmin,xmax, ymin,ymax, zmin,zmax]`. 1764 """ 1765 # OVERRIDE CommonAlgorithms.bounds() which is too slow 1766 return np.array(self.dataset.GetBounds()) 1767 1768 def isosurface(self, value=None) -> "vedo.Mesh": 1769 """ 1770 Return a `Mesh` isosurface extracted from the object. 1771 1772 Set `value` as single float or list of values to draw the isosurface(s). 1773 """ 1774 scrange = self.dataset.GetScalarRange() 1775 1776 cf = vtki.new("ContourFilter") 1777 cf.UseScalarTreeOn() 1778 cf.SetInputData(self.dataset) 1779 cf.ComputeNormalsOn() 1780 1781 if value is None: 1782 value = (2 * scrange[0] + scrange[1]) / 3.0 1783 # print("automatic isosurface value =", value) 1784 cf.SetValue(0, value) 1785 else: 1786 if utils.is_sequence(value): 1787 cf.SetNumberOfContours(len(value)) 1788 for i, t in enumerate(value): 1789 cf.SetValue(i, t) 1790 else: 1791 cf.SetValue(0, value) 1792 1793 cf.Update() 1794 poly = cf.GetOutput() 1795 1796 out = vedo.mesh.Mesh(poly, c=None).phong() 1797 out.mapper.SetScalarRange(scrange[0], scrange[1]) 1798 1799 out.pipeline = utils.OperationNode( 1800 "isosurface", 1801 parents=[self], 1802 comment=f"#pts {out.dataset.GetNumberOfPoints()}", 1803 c="#4cc9f0:#e9c46a", 1804 ) 1805 return out 1806 1807 def cut_with_plane(self, origin=(0, 0, 0), normal="x") -> "UnstructuredGrid": 1808 """ 1809 Cut the object with the plane defined by a point and a normal. 1810 1811 Arguments: 1812 origin : (list) 1813 the cutting plane goes through this point 1814 normal : (list, str) 1815 normal vector to the cutting plane 1816 """ 1817 strn = str(normal) 1818 if strn == "x": normal = (1, 0, 0) 1819 elif strn == "y": normal = (0, 1, 0) 1820 elif strn == "z": normal = (0, 0, 1) 1821 elif strn == "-x": normal = (-1, 0, 0) 1822 elif strn == "-y": normal = (0, -1, 0) 1823 elif strn == "-z": normal = (0, 0, -1) 1824 plane = vtki.new("Plane") 1825 plane.SetOrigin(origin) 1826 plane.SetNormal(normal) 1827 clipper = vtki.new("ClipDataSet") 1828 clipper.SetInputData(self.dataset) 1829 clipper.SetClipFunction(plane) 1830 clipper.GenerateClipScalarsOff() 1831 clipper.GenerateClippedOutputOff() 1832 clipper.SetValue(0) 1833 clipper.Update() 1834 cout = clipper.GetOutput() 1835 ug = vedo.UnstructuredGrid(cout) 1836 if isinstance(self, UnstructuredGrid): 1837 self._update(cout) 1838 self.pipeline = utils.OperationNode("cut_with_plane", parents=[self], c="#9e2a2b") 1839 return self 1840 ug.pipeline = utils.OperationNode("cut_with_plane", parents=[self], c="#9e2a2b") 1841 return ug 1842 1843 def cut_with_mesh(self, mesh, invert=False, whole_cells=False, on_boundary=False) -> "UnstructuredGrid": 1844 """ 1845 Cut a `RectilinearGrid` with a `Mesh`. 1846 1847 Use `invert` to return cut off part of the input object. 1848 """ 1849 ug = self.dataset 1850 1851 ippd = vtki.new("ImplicitPolyDataDistance") 1852 ippd.SetInput(mesh.dataset) 1853 1854 if whole_cells or on_boundary: 1855 clipper = vtki.new("ExtractGeometry") 1856 clipper.SetInputData(ug) 1857 clipper.SetImplicitFunction(ippd) 1858 clipper.SetExtractInside(not invert) 1859 clipper.SetExtractBoundaryCells(False) 1860 if on_boundary: 1861 clipper.SetExtractBoundaryCells(True) 1862 clipper.SetExtractOnlyBoundaryCells(True) 1863 else: 1864 signed_dists = vtki.vtkFloatArray() 1865 signed_dists.SetNumberOfComponents(1) 1866 signed_dists.SetName("SignedDistance") 1867 for pointId in range(ug.GetNumberOfPoints()): 1868 p = ug.GetPoint(pointId) 1869 signed_dist = ippd.EvaluateFunction(p) 1870 signed_dists.InsertNextValue(signed_dist) 1871 ug.GetPointData().AddArray(signed_dists) 1872 ug.GetPointData().SetActiveScalars("SignedDistance") # NEEDED 1873 clipper = vtki.new("ClipDataSet") 1874 clipper.SetInputData(ug) 1875 clipper.SetInsideOut(not invert) 1876 clipper.SetValue(0.0) 1877 1878 clipper.Update() 1879 1880 out = UnstructuredGrid(clipper.GetOutput()) 1881 out.pipeline = utils.OperationNode("cut_with_mesh", parents=[self], c="#9e2a2b") 1882 return out 1883 1884 1885########################################################################## 1886class StructuredGrid(PointAlgorithms, MeshVisual): 1887 """ 1888 Build a structured grid. 1889 """ 1890 1891 def __init__(self, inputobj=None): 1892 """ 1893 A StructuredGrid is a dataset where edges of the hexahedrons are 1894 not necessarily parallel to the coordinate axes. 1895 It can be thought of as a tessellation of a block of 3D space, 1896 similar to a `RectilinearGrid` 1897 except that the cells are not necessarily cubes, they can have different 1898 orientations but are connected in the same way as a `RectilinearGrid`. 1899 1900 Arguments: 1901 inputobj : (vtkStructuredGrid, list, str) 1902 list of points and tet indices, or filename 1903 1904 Example: 1905 ```python 1906 from vedo import * 1907 sgrid = StructuredGrid(dataurl+"structgrid.vts") 1908 print(sgrid) 1909 msh = sgrid.tomesh().lw(1) 1910 show(msh, axes=1, viewup="z") 1911 ``` 1912 1913 ```python 1914 from vedo import * 1915 1916 cx = np.sqrt(np.linspace(100, 400, 10)) 1917 cy = np.linspace(30, 40, 20) 1918 cz = np.linspace(40, 50, 30) 1919 x, y, z = np.meshgrid(cx, cy, cz) 1920 1921 sgrid1 = StructuredGrid([x, y, z]) 1922 sgrid1.cmap("viridis", sgrid1.vertices[:, 0]) 1923 print(sgrid1) 1924 1925 sgrid2 = sgrid1.clone().cut_with_plane(normal=(-1,1,1), origin=[14,34,44]) 1926 msh2 = sgrid2.tomesh(shrink=0.9).lw(1).cmap("viridis") 1927 1928 show( 1929 [["StructuredGrid", sgrid1], ["Shrinked Mesh", msh2]], 1930 N=2, axes=1, viewup="z", 1931 ) 1932 ``` 1933 """ 1934 1935 super().__init__() 1936 1937 self.dataset = None 1938 1939 self.mapper = vtki.new("PolyDataMapper") 1940 self._actor = vtki.vtkActor() 1941 self._actor.retrieve_object = weak_ref_to(self) 1942 self._actor.SetMapper(self.mapper) 1943 self.properties = self._actor.GetProperty() 1944 1945 self.transform = LinearTransform() 1946 self.point_locator = None 1947 self.cell_locator = None 1948 self.line_locator = None 1949 1950 self.name = "StructuredGrid" 1951 self.filename = "" 1952 1953 self.info = {} 1954 self.time = time.time() 1955 1956 ############################### 1957 if inputobj is None: 1958 self.dataset = vtki.vtkStructuredGrid() 1959 1960 elif isinstance(inputobj, vtki.vtkStructuredGrid): 1961 self.dataset = inputobj 1962 1963 elif isinstance(inputobj, StructuredGrid): 1964 self.dataset = inputobj.dataset 1965 1966 elif isinstance(inputobj, str): 1967 if "https://" in inputobj: 1968 inputobj = download(inputobj, verbose=False) 1969 if inputobj.endswith(".vts"): 1970 reader = vtki.new("XMLStructuredGridReader") 1971 else: 1972 reader = vtki.new("StructuredGridReader") 1973 self.filename = inputobj 1974 reader.SetFileName(inputobj) 1975 reader.Update() 1976 self.dataset = reader.GetOutput() 1977 1978 elif utils.is_sequence(inputobj): 1979 self.dataset = vtki.vtkStructuredGrid() 1980 x, y, z = inputobj 1981 xyz = np.vstack(( 1982 x.flatten(order="F"), 1983 y.flatten(order="F"), 1984 z.flatten(order="F")) 1985 ).T 1986 dims = x.shape 1987 self.dataset.SetDimensions(dims) 1988 # self.dataset.SetDimensions(dims[1], dims[0], dims[2]) 1989 vpoints = vtki.vtkPoints() 1990 vpoints.SetData(utils.numpy2vtk(xyz)) 1991 self.dataset.SetPoints(vpoints) 1992 1993 1994 ############################### 1995 if not self.dataset: 1996 vedo.logger.error(f"StructuredGrid: cannot understand input type {type(inputobj)}") 1997 return 1998 1999 self.properties.SetColor(0.352, 0.612, 0.996) # blue7 2000 2001 self.pipeline = utils.OperationNode( 2002 self, comment=f"#tets {self.dataset.GetNumberOfCells()}", c="#9e2a2b" 2003 ) 2004 2005 @property 2006 def actor(self): 2007 """Return the `vtkActor` of the object.""" 2008 gf = vtki.new("GeometryFilter") 2009 gf.SetInputData(self.dataset) 2010 gf.Update() 2011 self.mapper.SetInputData(gf.GetOutput()) 2012 self.mapper.Modified() 2013 return self._actor 2014 2015 @actor.setter 2016 def actor(self, _): 2017 pass 2018 2019 def _update(self, data, reset_locators=False): 2020 self.dataset = data 2021 if reset_locators: 2022 self.cell_locator = None 2023 self.point_locator = None 2024 return self 2025 2026 ################################################################## 2027 def __str__(self): 2028 """Print a summary for the `StructuredGrid` object.""" 2029 module = self.__class__.__module__ 2030 name = self.__class__.__name__ 2031 out = vedo.printc( 2032 f"{module}.{name} at ({hex(self.memory_address())})".ljust(75), 2033 c="c", bold=True, invert=True, return_string=True, 2034 ) 2035 out += "\x1b[0m\x1b[36;1m" 2036 2037 out += "name".ljust(14) + ": " + str(self.name) + "\n" 2038 if self.filename: 2039 out += "filename".ljust(14) + ": " + str(self.filename) + "\n" 2040 2041 out += "dimensions".ljust(14) + ": " + str(self.dataset.GetDimensions()) + "\n" 2042 2043 out += "center".ljust(14) + ": " 2044 out += utils.precision(self.dataset.GetCenter(), 6) + "\n" 2045 2046 bnds = self.bounds() 2047 bx1, bx2 = utils.precision(bnds[0], 3), utils.precision(bnds[1], 3) 2048 by1, by2 = utils.precision(bnds[2], 3), utils.precision(bnds[3], 3) 2049 bz1, bz2 = utils.precision(bnds[4], 3), utils.precision(bnds[5], 3) 2050 out += "bounds".ljust(14) + ":" 2051 out += " x=(" + bx1 + ", " + bx2 + ")," 2052 out += " y=(" + by1 + ", " + by2 + ")," 2053 out += " z=(" + bz1 + ", " + bz2 + ")\n" 2054 2055 out += "memory size".ljust(14) + ": " 2056 out += utils.precision(self.dataset.GetActualMemorySize() / 1024, 2) + " MB\n" 2057 2058 for key in self.pointdata.keys(): 2059 arr = self.pointdata[key] 2060 rng = utils.precision(arr.min(), 3) + ", " + utils.precision(arr.max(), 3) 2061 mark_active = "pointdata" 2062 a_scalars = self.dataset.GetPointData().GetScalars() 2063 a_vectors = self.dataset.GetPointData().GetVectors() 2064 a_tensors = self.dataset.GetPointData().GetTensors() 2065 if a_scalars and a_scalars.GetName() == key: 2066 mark_active += " *" 2067 elif a_vectors and a_vectors.GetName() == key: 2068 mark_active += " **" 2069 elif a_tensors and a_tensors.GetName() == key: 2070 mark_active += " ***" 2071 out += mark_active.ljust(14) + f': "{key}" ({arr.dtype}), ndim={arr.ndim}' 2072 out += f", range=({rng})\n" 2073 2074 for key in self.celldata.keys(): 2075 arr = self.celldata[key] 2076 rng = utils.precision(arr.min(), 3) + ", " + utils.precision(arr.max(), 3) 2077 mark_active = "celldata" 2078 a_scalars = self.dataset.GetCellData().GetScalars() 2079 a_vectors = self.dataset.GetCellData().GetVectors() 2080 a_tensors = self.dataset.GetCellData().GetTensors() 2081 if a_scalars and a_scalars.GetName() == key: 2082 mark_active += " *" 2083 elif a_vectors and a_vectors.GetName() == key: 2084 mark_active += " **" 2085 elif a_tensors and a_tensors.GetName() == key: 2086 mark_active += " ***" 2087 out += mark_active.ljust(14) + f': "{key}" ({arr.dtype}), ndim={arr.ndim}' 2088 out += f", range=({rng})\n" 2089 2090 for key in self.metadata.keys(): 2091 arr = self.metadata[key] 2092 out += "metadata".ljust(14) + ": " + f'"{key}" ({len(arr)} values)\n' 2093 2094 return out.rstrip() + "\x1b[0m" 2095 2096 def _repr_html_(self): 2097 """ 2098 HTML representation of the StructuredGrid object for Jupyter Notebooks. 2099 2100 Returns: 2101 HTML text with the image and some properties. 2102 """ 2103 import io 2104 import base64 2105 from PIL import Image 2106 2107 library_name = "vedo.grids.StructuredGrid" 2108 help_url = "https://vedo.embl.es/docs/vedo/grids.html#StructuredGrid" 2109 2110 m = self.tomesh().linewidth(1).lighting("off") 2111 arr= m.thumbnail(zoom=1, elevation=-30, azimuth=-30) 2112 2113 im = Image.fromarray(arr) 2114 buffered = io.BytesIO() 2115 im.save(buffered, format="PNG", quality=100) 2116 encoded = base64.b64encode(buffered.getvalue()).decode("utf-8") 2117 url = "data:image/png;base64," + encoded 2118 image = f"<img src='{url}'></img>" 2119 2120 bounds = "<br/>".join( 2121 [ 2122 utils.precision(min_x, 4) + " ... " + utils.precision(max_x, 4) 2123 for min_x, max_x in zip(self.bounds()[::2], self.bounds()[1::2]) 2124 ] 2125 ) 2126 2127 help_text = "" 2128 if self.name: 2129 help_text += f"<b> {self.name}:   </b>" 2130 help_text += '<b><a href="' + help_url + '" target="_blank">' + library_name + "</a></b>" 2131 if self.filename: 2132 dots = "" 2133 if len(self.filename) > 30: 2134 dots = "..." 2135 help_text += f"<br/><code><i>({dots}{self.filename[-30:]})</i></code>" 2136 2137 pdata = "" 2138 if self.dataset.GetPointData().GetScalars(): 2139 if self.dataset.GetPointData().GetScalars().GetName(): 2140 name = self.dataset.GetPointData().GetScalars().GetName() 2141 pdata = "<tr><td><b> point data array </b></td><td>" + name + "</td></tr>" 2142 2143 cdata = "" 2144 if self.dataset.GetCellData().GetScalars(): 2145 if self.dataset.GetCellData().GetScalars().GetName(): 2146 name = self.dataset.GetCellData().GetScalars().GetName() 2147 cdata = "<tr><td><b> cell data array </b></td><td>" + name + "</td></tr>" 2148 2149 pts = self.vertices 2150 cm = np.mean(pts, axis=0) 2151 2152 all = [ 2153 "<table>", 2154 "<tr>", 2155 "<td>", image, "</td>", 2156 "<td style='text-align: center; vertical-align: center;'><br/>", help_text, 2157 "<table>", 2158 "<tr><td><b> bounds </b> <br/> (x/y/z) </td><td>" + str(bounds) + "</td></tr>", 2159 "<tr><td><b> center of mass </b></td><td>" + utils.precision(cm,3) + "</td></tr>", 2160 "<tr><td><b> nr. points / cells </b></td><td>" 2161 + str(self.npoints) + " / " + str(self.ncells) + "</td></tr>", 2162 pdata, 2163 cdata, 2164 "</table>", 2165 "</table>", 2166 ] 2167 return "\n".join(all) 2168 2169 def dimensions(self) -> np.ndarray: 2170 """Return the number of points in the x, y and z directions.""" 2171 return np.array(self.dataset.GetDimensions()) 2172 2173 def clone(self, deep=True) -> "StructuredGrid": 2174 """Return a clone copy of the StructuredGrid. Alias of `copy()`.""" 2175 if deep: 2176 newrg = vtki.vtkStructuredGrid() 2177 newrg.CopyStructure(self.dataset) 2178 newrg.CopyAttributes(self.dataset) 2179 newvol = StructuredGrid(newrg) 2180 else: 2181 newvol = StructuredGrid(self.dataset) 2182 2183 prop = vtki.vtkProperty() 2184 prop.DeepCopy(self.properties) 2185 newvol.actor.SetProperty(prop) 2186 newvol.properties = prop 2187 newvol.pipeline = utils.OperationNode("clone", parents=[self], c="#bbd0ff", shape="diamond") 2188 return newvol 2189 2190 def find_point(self, x: list) -> int: 2191 """Given a position `x`, return the id of the closest point.""" 2192 return self.dataset.FindPoint(x) 2193 2194 def find_cell(self, x: list) -> dict: 2195 """Given a position `x`, return the id of the closest cell.""" 2196 cell = vtki.vtkHexagonalPrism() 2197 cellid = vtki.mutable(0) 2198 tol2 = 0.001 # vtki.mutable(0) 2199 subid = vtki.mutable(0) 2200 pcoords = [0.0, 0.0, 0.0] 2201 weights = [0.0, 0.0, 0.0] 2202 res = self.dataset.FindCell(x, cell, cellid, tol2, subid, pcoords, weights) 2203 result = {} 2204 result["cellid"] = cellid 2205 result["subid"] = subid 2206 result["pcoords"] = pcoords 2207 result["weights"] = weights 2208 result["status"] = res 2209 return result 2210 2211 def cut_with_plane(self, origin=(0, 0, 0), normal="x") -> "vedo.UnstructuredGrid": 2212 """ 2213 Cut the object with the plane defined by a point and a normal. 2214 2215 Arguments: 2216 origin : (list) 2217 the cutting plane goes through this point 2218 normal : (list, str) 2219 normal vector to the cutting plane 2220 2221 Returns an `UnstructuredGrid` object. 2222 """ 2223 strn = str(normal) 2224 if strn == "x": normal = (1, 0, 0) 2225 elif strn == "y": normal = (0, 1, 0) 2226 elif strn == "z": normal = (0, 0, 1) 2227 elif strn == "-x": normal = (-1, 0, 0) 2228 elif strn == "-y": normal = (0, -1, 0) 2229 elif strn == "-z": normal = (0, 0, -1) 2230 plane = vtki.new("Plane") 2231 plane.SetOrigin(origin) 2232 plane.SetNormal(normal) 2233 clipper = vtki.new("ClipDataSet") 2234 clipper.SetInputData(self.dataset) 2235 clipper.SetClipFunction(plane) 2236 clipper.GenerateClipScalarsOff() 2237 clipper.GenerateClippedOutputOff() 2238 clipper.SetValue(0) 2239 clipper.Update() 2240 cout = clipper.GetOutput() 2241 ug = vedo.UnstructuredGrid(cout) 2242 if isinstance(self, vedo.UnstructuredGrid): 2243 self._update(cout) 2244 self.pipeline = utils.OperationNode("cut_with_plane", parents=[self], c="#9e2a2b") 2245 return self 2246 ug.pipeline = utils.OperationNode("cut_with_plane", parents=[self], c="#9e2a2b") 2247 return ug 2248 2249 def cut_with_mesh(self, mesh: Mesh, invert=False, whole_cells=False, on_boundary=False) -> "UnstructuredGrid": 2250 """ 2251 Cut a `RectilinearGrid` with a `Mesh`. 2252 2253 Use `invert` to return cut off part of the input object. 2254 2255 Returns an `UnstructuredGrid` object. 2256 """ 2257 ug = self.dataset 2258 2259 ippd = vtki.new("ImplicitPolyDataDistance") 2260 ippd.SetInput(mesh.dataset) 2261 2262 if whole_cells or on_boundary: 2263 clipper = vtki.new("ExtractGeometry") 2264 clipper.SetInputData(ug) 2265 clipper.SetImplicitFunction(ippd) 2266 clipper.SetExtractInside(not invert) 2267 clipper.SetExtractBoundaryCells(False) 2268 if on_boundary: 2269 clipper.SetExtractBoundaryCells(True) 2270 clipper.SetExtractOnlyBoundaryCells(True) 2271 else: 2272 signed_dists = vtki.vtkFloatArray() 2273 signed_dists.SetNumberOfComponents(1) 2274 signed_dists.SetName("SignedDistance") 2275 for pointId in range(ug.GetNumberOfPoints()): 2276 p = ug.GetPoint(pointId) 2277 signed_dist = ippd.EvaluateFunction(p) 2278 signed_dists.InsertNextValue(signed_dist) 2279 ug.GetPointData().AddArray(signed_dists) 2280 ug.GetPointData().SetActiveScalars("SignedDistance") # NEEDED 2281 clipper = vtki.new("ClipDataSet") 2282 clipper.SetInputData(ug) 2283 clipper.SetInsideOut(not invert) 2284 clipper.SetValue(0.0) 2285 2286 clipper.Update() 2287 2288 out = UnstructuredGrid(clipper.GetOutput()) 2289 out.pipeline = utils.OperationNode("cut_with_mesh", parents=[self], c="#9e2a2b") 2290 return out 2291 2292 def isosurface(self, value=None) -> "vedo.Mesh": 2293 """ 2294 Return a `Mesh` isosurface extracted from the object. 2295 2296 Set `value` as single float or list of values to draw the isosurface(s). 2297 """ 2298 scrange = self.dataset.GetScalarRange() 2299 2300 cf = vtki.new("ContourFilter") 2301 cf.UseScalarTreeOn() 2302 cf.SetInputData(self.dataset) 2303 cf.ComputeNormalsOn() 2304 2305 if value is None: 2306 value = (2 * scrange[0] + scrange[1]) / 3.0 2307 # print("automatic isosurface value =", value) 2308 cf.SetValue(0, value) 2309 else: 2310 if utils.is_sequence(value): 2311 cf.SetNumberOfContours(len(value)) 2312 for i, t in enumerate(value): 2313 cf.SetValue(i, t) 2314 else: 2315 cf.SetValue(0, value) 2316 2317 cf.Update() 2318 poly = cf.GetOutput() 2319 2320 out = vedo.mesh.Mesh(poly, c=None).phong() 2321 out.mapper.SetScalarRange(scrange[0], scrange[1]) 2322 2323 out.pipeline = utils.OperationNode( 2324 "isosurface", 2325 parents=[self], 2326 comment=f"#pts {out.dataset.GetNumberOfPoints()}", 2327 c="#4cc9f0:#e9c46a", 2328 ) 2329 return out
38class UnstructuredGrid(PointAlgorithms, MeshVisual): 39 """Support for UnstructuredGrid objects.""" 40 41 def __init__(self, inputobj=None): 42 """ 43 Support for UnstructuredGrid objects. 44 45 Arguments: 46 inputobj : (list, vtkUnstructuredGrid, str) 47 A list in the form `[points, cells, celltypes]`, 48 or a vtkUnstructuredGrid object, or a filename 49 50 Celltypes are identified by the following 51 [convention](https://vtk.org/doc/nightly/html/vtkCellType_8h_source.html). 52 """ 53 super().__init__() 54 55 self.dataset = None 56 57 self.mapper = vtki.new("PolyDataMapper") 58 self._actor = vtki.vtkActor() 59 self._actor.retrieve_object = weak_ref_to(self) 60 self._actor.SetMapper(self.mapper) 61 self.properties = self._actor.GetProperty() 62 63 self.transform = LinearTransform() 64 self.point_locator = None 65 self.cell_locator = None 66 self.line_locator = None 67 68 self.name = "UnstructuredGrid" 69 self.filename = "" 70 self.file_size = "" 71 72 self.info = {} 73 self.time = time.time() 74 self.rendered_at = set() 75 76 ###################################### 77 inputtype = str(type(inputobj)) 78 79 if inputobj is None: 80 self.dataset = vtki.vtkUnstructuredGrid() 81 82 elif utils.is_sequence(inputobj): 83 84 pts, cells, celltypes = inputobj 85 assert len(cells) == len(celltypes) 86 87 self.dataset = vtki.vtkUnstructuredGrid() 88 89 if not utils.is_sequence(cells[0]): 90 tets = [] 91 nf = cells[0] + 1 92 for i, cl in enumerate(cells): 93 if i in (nf, 0): 94 k = i + 1 95 nf = cl + k 96 cell = [cells[j + k] for j in range(cl)] 97 tets.append(cell) 98 cells = tets 99 100 # This would fill the points and use those to define orientation 101 vpts = utils.numpy2vtk(pts, dtype=np.float32) 102 points = vtki.vtkPoints() 103 points.SetData(vpts) 104 self.dataset.SetPoints(points) 105 106 # Fill cells 107 # https://vtk.org/doc/nightly/html/vtkCellType_8h_source.html 108 for i, ct in enumerate(celltypes): 109 if ct == vtki.cell_types["VERTEX"]: 110 cell = vtki.vtkVertex() 111 elif ct == vtki.cell_types["POLY_VERTEX"]: 112 cell = vtki.vtkPolyVertex() 113 elif ct == vtki.cell_types["TETRA"]: 114 cell = vtki.vtkTetra() 115 elif ct == vtki.cell_types["WEDGE"]: 116 cell = vtki.vtkWedge() 117 elif ct == vtki.cell_types["LINE"]: 118 cell = vtki.vtkLine() 119 elif ct == vtki.cell_types["POLY_LINE"]: 120 cell = vtki.vtkPolyLine() 121 elif ct == vtki.cell_types["TRIANGLE"]: 122 cell = vtki.vtkTriangle() 123 elif ct == vtki.cell_types["TRIANGLE_STRIP"]: 124 cell = vtki.vtkTriangleStrip() 125 elif ct == vtki.cell_types["POLYGON"]: 126 cell = vtki.vtkPolygon() 127 elif ct == vtki.cell_types["PIXEL"]: 128 cell = vtki.vtkPixel() 129 elif ct == vtki.cell_types["QUAD"]: 130 cell = vtki.vtkQuad() 131 elif ct == vtki.cell_types["VOXEL"]: 132 cell = vtki.vtkVoxel() 133 elif ct == vtki.cell_types["PYRAMID"]: 134 cell = vtki.vtkPyramid() 135 elif ct == vtki.cell_types["HEXAHEDRON"]: 136 cell = vtki.vtkHexahedron() 137 elif ct == vtki.cell_types["HEXAGONAL_PRISM"]: 138 cell = vtki.vtkHexagonalPrism() 139 elif ct == vtki.cell_types["PENTAGONAL_PRISM"]: 140 cell = vtki.vtkPentagonalPrism() 141 elif ct == vtki.cell_types["QUADRATIC_TETRA"]: 142 from vtkmodules.vtkCommonDataModel import vtkQuadraticTetra 143 cell = vtkQuadraticTetra() 144 elif ct == vtki.cell_types["QUADRATIC_HEXAHEDRON"]: 145 from vtkmodules.vtkCommonDataModel import vtkQuadraticHexahedron 146 cell = vtkQuadraticHexahedron() 147 elif ct == vtki.cell_types["QUADRATIC_WEDGE"]: 148 from vtkmodules.vtkCommonDataModel import vtkQuadraticWedge 149 cell = vtkQuadraticWedge() 150 elif ct == vtki.cell_types["QUADRATIC_PYRAMID"]: 151 from vtkmodules.vtkCommonDataModel import vtkQuadraticPyramid 152 cell = vtkQuadraticPyramid() 153 elif ct == vtki.cell_types["QUADRATIC_LINEAR_QUAD"]: 154 from vtkmodules.vtkCommonDataModel import vtkQuadraticLinearQuad 155 cell = vtkQuadraticLinearQuad() 156 elif ct == vtki.cell_types["QUADRATIC_LINEAR_WEDGE"]: 157 from vtkmodules.vtkCommonDataModel import vtkQuadraticLinearWedge 158 cell = vtkQuadraticLinearWedge() 159 elif ct == vtki.cell_types["BIQUADRATIC_QUADRATIC_WEDGE"]: 160 from vtkmodules.vtkCommonDataModel import vtkBiQuadraticQuadraticWedge 161 cell = vtkBiQuadraticQuadraticWedge() 162 elif ct == vtki.cell_types["BIQUADRATIC_QUADRATIC_HEXAHEDRON"]: 163 from vtkmodules.vtkCommonDataModel import vtkBiQuadraticQuadraticHexahedron 164 cell = vtkBiQuadraticQuadraticHexahedron() 165 elif ct == vtki.cell_types["BIQUADRATIC_TRIANGLE"]: 166 from vtkmodules.vtkCommonDataModel import vtkBiQuadraticTriangle 167 cell = vtkBiQuadraticTriangle() 168 elif ct == vtki.cell_types["CUBIC_LINE"]: 169 from vtkmodules.vtkCommonDataModel import vtkCubicLine 170 cell = vtkCubicLine() 171 elif ct == vtki.cell_types["CONVEX_POINT_SET"]: 172 from vtkmodules.vtkCommonDataModel import vtkConvexPointSet 173 cell = vtkConvexPointSet() 174 elif ct == vtki.cell_types["POLYHEDRON"]: 175 from vtkmodules.vtkCommonDataModel import vtkPolyhedron 176 cell = vtkPolyhedron() 177 elif ct == vtki.cell_types["HIGHER_ORDER_TRIANGLE"]: 178 from vtkmodules.vtkCommonDataModel import vtkHigherOrderTriangle 179 cell = vtkHigherOrderTriangle() 180 elif ct == vtki.cell_types["HIGHER_ORDER_QUAD"]: 181 from vtkmodules.vtkCommonDataModel import vtkHigherOrderQuadrilateral 182 cell = vtkHigherOrderQuadrilateral() 183 elif ct == vtki.cell_types["HIGHER_ORDER_TETRAHEDRON"]: 184 from vtkmodules.vtkCommonDataModel import vtkHigherOrderTetra 185 cell = vtkHigherOrderTetra() 186 elif ct == vtki.cell_types["HIGHER_ORDER_WEDGE"]: 187 from vtkmodules.vtkCommonDataModel import vtkHigherOrderWedge 188 cell = vtkHigherOrderWedge() 189 elif ct == vtki.cell_types["HIGHER_ORDER_HEXAHEDRON"]: 190 from vtkmodules.vtkCommonDataModel import vtkHigherOrderHexahedron 191 cell = vtkHigherOrderHexahedron() 192 elif ct == vtki.cell_types["LAGRANGE_CURVE"]: 193 from vtkmodules.vtkCommonDataModel import vtkLagrangeCurve 194 cell = vtkLagrangeCurve() 195 elif ct == vtki.cell_types["LAGRANGE_TRIANGLE"]: 196 from vtkmodules.vtkCommonDataModel import vtkLagrangeTriangle 197 cell = vtkLagrangeTriangle() 198 elif ct == vtki.cell_types["LAGRANGE_QUADRILATERAL"]: 199 from vtkmodules.vtkCommonDataModel import vtkLagrangeQuadrilateral 200 cell = vtkLagrangeQuadrilateral() 201 elif ct == vtki.cell_types["LAGRANGE_TETRAHEDRON"]: 202 from vtkmodules.vtkCommonDataModel import vtkLagrangeTetra 203 cell = vtkLagrangeTetra() 204 elif ct == vtki.cell_types["LAGRANGE_HEXAHEDRON"]: 205 from vtkmodules.vtkCommonDataModel import vtkLagrangeHexahedron 206 cell = vtkLagrangeHexahedron() 207 elif ct == vtki.cell_types["LAGRANGE_WEDGE"]: 208 from vtkmodules.vtkCommonDataModel import vtkLagrangeWedge 209 cell = vtkLagrangeWedge() 210 elif ct == vtki.cell_types["BEZIER_CURVE"]: 211 from vtkmodules.vtkCommonDataModel import vtkBezierCurve 212 cell = vtkBezierCurve() 213 elif ct == vtki.cell_types["BEZIER_TRIANGLE"]: 214 from vtkmodules.vtkCommonDataModel import vtkBezierTriangle 215 cell = vtkBezierTriangle() 216 elif ct == vtki.cell_types["BEZIER_QUADRILATERAL"]: 217 from vtkmodules.vtkCommonDataModel import vtkBezierQuadrilateral 218 cell = vtkBezierQuadrilateral() 219 elif ct == vtki.cell_types["BEZIER_TETRAHEDRON"]: 220 from vtkmodules.vtkCommonDataModel import vtkBezierTetra 221 cell = vtkBezierTetra() 222 elif ct == vtki.cell_types["BEZIER_HEXAHEDRON"]: 223 from vtkmodules.vtkCommonDataModel import vtkBezierHexahedron 224 cell = vtkBezierHexahedron() 225 elif ct == vtki.cell_types["BEZIER_WEDGE"]: 226 from vtkmodules.vtkCommonDataModel import vtkBezierWedge 227 cell = vtkBezierWedge() 228 else: 229 vedo.logger.error( 230 f"UnstructuredGrid: cell type {ct} not supported. Skip.") 231 continue 232 233 cpids = cell.GetPointIds() 234 cell_conn = cells[i] 235 for j, pid in enumerate(cell_conn): 236 cpids.SetId(j, pid) 237 self.dataset.InsertNextCell(ct, cpids) 238 239 elif "UnstructuredGrid" in inputtype: 240 self.dataset = inputobj 241 242 elif isinstance(inputobj, str): 243 if "https://" in inputobj: 244 inputobj = download(inputobj, verbose=False) 245 self.filename = inputobj 246 if inputobj.endswith(".vtu"): 247 reader = vtki.new("XMLUnstructuredGridReader") 248 else: 249 reader = vtki.new("UnstructuredGridReader") 250 self.filename = inputobj 251 reader.SetFileName(inputobj) 252 reader.Update() 253 self.dataset = reader.GetOutput() 254 255 else: 256 # this converts other types of vtk objects to UnstructuredGrid 257 apf = vtki.new("AppendFilter") 258 try: 259 apf.AddInputData(inputobj) 260 except TypeError: 261 apf.AddInputData(inputobj.dataset) 262 apf.Update() 263 self.dataset = apf.GetOutput() 264 265 self.properties.SetColor(0.89, 0.455, 0.671) # pink7 266 267 self.pipeline = utils.OperationNode( 268 self, comment=f"#cells {self.dataset.GetNumberOfCells()}", c="#4cc9f0" 269 ) 270 271 # ------------------------------------------------------------------ 272 def __str__(self): 273 """Print a string summary of the `UnstructuredGrid` object.""" 274 module = self.__class__.__module__ 275 name = self.__class__.__name__ 276 out = vedo.printc( 277 f"{module}.{name} at ({hex(self.memory_address())})".ljust(75), 278 c="m", bold=True, invert=True, return_string=True, 279 ) 280 out += "\x1b[0m\u001b[35m" 281 282 out += "nr. of verts".ljust(14) + ": " + str(self.npoints) + "\n" 283 out += "nr. of cells".ljust(14) + ": " + str(self.ncells) + "\n" 284 ct_arr = np.unique(self.cell_types_array) 285 cnames = [k for k, v in vtki.cell_types.items() if v in ct_arr] 286 out += "cell types".ljust(14) + ": " + str(cnames) + "\n" 287 288 if self.npoints: 289 out+="size".ljust(14)+ ": average=" + utils.precision(self.average_size(),6) 290 out+=", diagonal="+ utils.precision(self.diagonal_size(), 6)+ "\n" 291 out+="center of mass".ljust(14) + ": " + utils.precision(self.center_of_mass(),6)+"\n" 292 293 bnds = self.bounds() 294 bx1, bx2 = utils.precision(bnds[0], 3), utils.precision(bnds[1], 3) 295 by1, by2 = utils.precision(bnds[2], 3), utils.precision(bnds[3], 3) 296 bz1, bz2 = utils.precision(bnds[4], 3), utils.precision(bnds[5], 3) 297 out += "bounds".ljust(14) + ":" 298 out += " x=(" + bx1 + ", " + bx2 + ")," 299 out += " y=(" + by1 + ", " + by2 + ")," 300 out += " z=(" + bz1 + ", " + bz2 + ")\n" 301 302 for key in self.pointdata.keys(): 303 arr = self.pointdata[key] 304 rng = utils.precision(arr.min(), 3) + ", " + utils.precision(arr.max(), 3) 305 mark_active = "pointdata" 306 a_scalars = self.dataset.GetPointData().GetScalars() 307 a_vectors = self.dataset.GetPointData().GetVectors() 308 a_tensors = self.dataset.GetPointData().GetTensors() 309 if a_scalars and a_scalars.GetName() == key: 310 mark_active += " *" 311 elif a_vectors and a_vectors.GetName() == key: 312 mark_active += " **" 313 elif a_tensors and a_tensors.GetName() == key: 314 mark_active += " ***" 315 out += mark_active.ljust(14) + f': "{key}" ({arr.dtype}), ndim={arr.ndim}' 316 out += f", range=({rng})\n" 317 318 for key in self.celldata.keys(): 319 arr = self.celldata[key] 320 rng = utils.precision(arr.min(), 3) + ", " + utils.precision(arr.max(), 3) 321 mark_active = "celldata" 322 a_scalars = self.dataset.GetCellData().GetScalars() 323 a_vectors = self.dataset.GetCellData().GetVectors() 324 a_tensors = self.dataset.GetCellData().GetTensors() 325 if a_scalars and a_scalars.GetName() == key: 326 mark_active += " *" 327 elif a_vectors and a_vectors.GetName() == key: 328 mark_active += " **" 329 elif a_tensors and a_tensors.GetName() == key: 330 mark_active += " ***" 331 out += mark_active.ljust(14) + f': "{key}" ({arr.dtype}), ndim={arr.ndim}' 332 out += f", range=({rng})\n" 333 334 for key in self.metadata.keys(): 335 arr = self.metadata[key] 336 out += "metadata".ljust(14) + ": " + f'"{key}" ({len(arr)} values)\n' 337 338 return out.rstrip() + "\x1b[0m" 339 340 def _repr_html_(self): 341 """ 342 HTML representation of the UnstructuredGrid object for Jupyter Notebooks. 343 344 Returns: 345 HTML text with the image and some properties. 346 """ 347 import io 348 import base64 349 from PIL import Image 350 351 library_name = "vedo.grids.UnstructuredGrid" 352 help_url = "https://vedo.embl.es/docs/vedo/grids.html#UnstructuredGrid" 353 354 arr = self.thumbnail() 355 im = Image.fromarray(arr) 356 buffered = io.BytesIO() 357 im.save(buffered, format="PNG", quality=100) 358 encoded = base64.b64encode(buffered.getvalue()).decode("utf-8") 359 url = "data:image/png;base64," + encoded 360 image = f"<img src='{url}'></img>" 361 362 bounds = "<br/>".join( 363 [ 364 utils.precision(min_x, 4) + " ... " + utils.precision(max_x, 4) 365 for min_x, max_x in zip(self.bounds()[::2], self.bounds()[1::2]) 366 ] 367 ) 368 369 help_text = "" 370 if self.name: 371 help_text += f"<b> {self.name}:   </b>" 372 help_text += '<b><a href="' + help_url + '" target="_blank">' + library_name + "</a></b>" 373 if self.filename: 374 dots = "" 375 if len(self.filename) > 30: 376 dots = "..." 377 help_text += f"<br/><code><i>({dots}{self.filename[-30:]})</i></code>" 378 379 pdata = "" 380 if self.dataset.GetPointData().GetScalars(): 381 if self.dataset.GetPointData().GetScalars().GetName(): 382 name = self.dataset.GetPointData().GetScalars().GetName() 383 pdata = "<tr><td><b> point data array </b></td><td>" + name + "</td></tr>" 384 385 cdata = "" 386 if self.dataset.GetCellData().GetScalars(): 387 if self.dataset.GetCellData().GetScalars().GetName(): 388 name = self.dataset.GetCellData().GetScalars().GetName() 389 cdata = "<tr><td><b> cell data array </b></td><td>" + name + "</td></tr>" 390 391 pts = self.vertices 392 cm = np.mean(pts, axis=0) 393 394 all = [ 395 "<table>", 396 "<tr>", 397 "<td>", image, "</td>", 398 "<td style='text-align: center; vertical-align: center;'><br/>", help_text, 399 "<table>", 400 "<tr><td><b> bounds </b> <br/> (x/y/z) </td><td>" + str(bounds) + "</td></tr>", 401 "<tr><td><b> center of mass </b></td><td>" + utils.precision(cm,3) + "</td></tr>", 402 # "<tr><td><b> average size </b></td><td>" + str(average_size) + "</td></tr>", 403 "<tr><td><b> nr. points / cells </b></td><td>" 404 + str(self.npoints) + " / " + str(self.ncells) + "</td></tr>", 405 pdata, 406 cdata, 407 "</table>", 408 "</table>", 409 ] 410 return "\n".join(all) 411 412 @property 413 def actor(self): 414 """Return the `vtkActor` of the object.""" 415 # print("building actor") 416 gf = vtki.new("GeometryFilter") 417 gf.SetInputData(self.dataset) 418 gf.Update() 419 out = gf.GetOutput() 420 self.mapper.SetInputData(out) 421 self.mapper.Modified() 422 return self._actor 423 424 @actor.setter 425 def actor(self, _): 426 pass 427 428 def _update(self, data, reset_locators=False): 429 self.dataset = data 430 if reset_locators: 431 self.cell_locator = None 432 self.point_locator = None 433 return self 434 435 def merge(self, *others) -> Self: 436 """ 437 Merge multiple datasets into one single `UnstrcturedGrid`. 438 """ 439 apf = vtki.new("AppendFilter") 440 for o in others: 441 if isinstance(o, UnstructuredGrid): 442 apf.AddInputData(o.dataset) 443 elif isinstance(o, vtki.vtkUnstructuredGrid): 444 apf.AddInputData(o) 445 else: 446 vedo.printc("Error: cannot merge type", type(o), c="r") 447 apf.Update() 448 self._update(apf.GetOutput()) 449 self.pipeline = utils.OperationNode( 450 "merge", parents=[self, *others], c="#9e2a2b" 451 ) 452 return self 453 454 def copy(self, deep=True) -> "UnstructuredGrid": 455 """Return a copy of the object. Alias of `clone()`.""" 456 return self.clone(deep=deep) 457 458 def clone(self, deep=True) -> "UnstructuredGrid": 459 """Clone the UnstructuredGrid object to yield an exact copy.""" 460 ug = vtki.vtkUnstructuredGrid() 461 if deep: 462 ug.DeepCopy(self.dataset) 463 else: 464 ug.ShallowCopy(self.dataset) 465 if isinstance(self, vedo.UnstructuredGrid): 466 cloned = vedo.UnstructuredGrid(ug) 467 else: 468 cloned = vedo.TetMesh(ug) 469 470 cloned.copy_properties_from(self) 471 472 cloned.pipeline = utils.OperationNode( 473 "clone", parents=[self], shape="diamond", c="#bbe1ed" 474 ) 475 return cloned 476 477 def bounds(self) -> np.ndarray: 478 """ 479 Get the object bounds. 480 Returns a list in format `[xmin,xmax, ymin,ymax, zmin,zmax]`. 481 """ 482 # OVERRIDE CommonAlgorithms.bounds() which is too slow 483 return np.array(self.dataset.GetBounds()) 484 485 def threshold(self, name=None, above=None, below=None, on="cells") -> Self: 486 """ 487 Threshold the tetrahedral mesh by a cell scalar value. 488 Reduce to only tets which satisfy the threshold limits. 489 490 - if `above = below` will only select tets with that specific value. 491 - if `above > below` selection range is flipped. 492 493 Set keyword "on" to either "cells" or "points". 494 """ 495 th = vtki.new("Threshold") 496 th.SetInputData(self.dataset) 497 498 if name is None: 499 if self.celldata.keys(): 500 name = self.celldata.keys()[0] 501 th.SetInputArrayToProcess(0, 0, 0, 1, name) 502 elif self.pointdata.keys(): 503 name = self.pointdata.keys()[0] 504 th.SetInputArrayToProcess(0, 0, 0, 0, name) 505 if name is None: 506 vedo.logger.warning("cannot find active array. Skip.") 507 return self 508 else: 509 if on.startswith("c"): 510 th.SetInputArrayToProcess(0, 0, 0, 1, name) 511 else: 512 th.SetInputArrayToProcess(0, 0, 0, 0, name) 513 514 if above is not None: 515 th.SetLowerThreshold(above) 516 517 if below is not None: 518 th.SetUpperThreshold(below) 519 520 th.Update() 521 return self._update(th.GetOutput()) 522 523 def isosurface(self, value=None, flying_edges=False) -> "vedo.mesh.Mesh": 524 """ 525 Return an `Mesh` isosurface extracted from the `Volume` object. 526 527 Set `value` as single float or list of values to draw the isosurface(s). 528 Use flying_edges for faster results (but sometimes can interfere with `smooth()`). 529 530 Examples: 531 - [isosurfaces1.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/isosurfaces1.py) 532 533  534 """ 535 scrange = self.dataset.GetScalarRange() 536 537 if flying_edges: 538 cf = vtki.new("FlyingEdges3D") 539 cf.InterpolateAttributesOn() 540 else: 541 cf = vtki.new("ContourFilter") 542 cf.UseScalarTreeOn() 543 544 cf.SetInputData(self.dataset) 545 cf.ComputeNormalsOn() 546 547 if utils.is_sequence(value): 548 cf.SetNumberOfContours(len(value)) 549 for i, t in enumerate(value): 550 cf.SetValue(i, t) 551 else: 552 if value is None: 553 value = (2 * scrange[0] + scrange[1]) / 3.0 554 # print("automatic isosurface value =", value) 555 cf.SetValue(0, value) 556 557 cf.Update() 558 poly = cf.GetOutput() 559 560 out = vedo.mesh.Mesh(poly, c=None).flat() 561 out.mapper.SetScalarRange(scrange[0], scrange[1]) 562 563 out.pipeline = utils.OperationNode( 564 "isosurface", 565 parents=[self], 566 comment=f"#pts {out.dataset.GetNumberOfPoints()}", 567 c="#4cc9f0:#e9c46a", 568 ) 569 return out 570 571 def shrink(self, fraction=0.8) -> Self: 572 """ 573 Shrink the individual cells. 574 575  576 """ 577 sf = vtki.new("ShrinkFilter") 578 sf.SetInputData(self.dataset) 579 sf.SetShrinkFactor(fraction) 580 sf.Update() 581 out = sf.GetOutput() 582 self._update(out) 583 self.pipeline = utils.OperationNode( 584 "shrink", comment=f"by {fraction}", parents=[self], c="#9e2a2b" 585 ) 586 return self 587 588 def tomesh(self, fill=False, shrink=1.0) -> "vedo.mesh.Mesh": 589 """ 590 Build a polygonal `Mesh` from the current object. 591 592 If `fill=True`, the interior faces of all the cells are created. 593 (setting a `shrink` value slightly smaller than the default 1.0 594 can avoid flickering due to internal adjacent faces). 595 596 If `fill=False`, only the boundary faces will be generated. 597 """ 598 gf = vtki.new("GeometryFilter") 599 if fill: 600 sf = vtki.new("ShrinkFilter") 601 sf.SetInputData(self.dataset) 602 sf.SetShrinkFactor(shrink) 603 sf.Update() 604 gf.SetInputData(sf.GetOutput()) 605 gf.Update() 606 poly = gf.GetOutput() 607 else: 608 gf.SetInputData(self.dataset) 609 gf.Update() 610 poly = gf.GetOutput() 611 612 msh = vedo.mesh.Mesh(poly) 613 msh.copy_properties_from(self) 614 msh.pipeline = utils.OperationNode( 615 "tomesh", parents=[self], comment=f"fill={fill}", c="#9e2a2b:#e9c46a" 616 ) 617 return msh 618 619 @property 620 def cell_types_array(self): 621 """Return the list of cell types in the dataset.""" 622 uarr = self.dataset.GetCellTypesArray() 623 return utils.vtk2numpy(uarr) 624 625 def extract_cells_by_type(self, ctype) -> "UnstructuredGrid": 626 """Extract a specific cell type and return a new `UnstructuredGrid`.""" 627 if isinstance(ctype, str): 628 try: 629 ctype = vtki.cell_types[ctype.upper()] 630 except KeyError: 631 vedo.logger.error(f"extract_cells_by_type: cell type {ctype} does not exist. Skip.") 632 return self 633 uarr = self.dataset.GetCellTypesArray() 634 ctarrtyp = np.where(utils.vtk2numpy(uarr) == ctype)[0] 635 uarrtyp = utils.numpy2vtk(ctarrtyp, deep=False, dtype="id") 636 selection_node = vtki.new("SelectionNode") 637 selection_node.SetFieldType(vtki.get_class("SelectionNode").CELL) 638 selection_node.SetContentType(vtki.get_class("SelectionNode").INDICES) 639 selection_node.SetSelectionList(uarrtyp) 640 selection = vtki.new("Selection") 641 selection.AddNode(selection_node) 642 es = vtki.new("ExtractSelection") 643 es.SetInputData(0, self.dataset) 644 es.SetInputData(1, selection) 645 es.Update() 646 647 ug = UnstructuredGrid(es.GetOutput()) 648 ug.pipeline = utils.OperationNode( 649 "extract_cell_type", comment=f"type {ctype}", c="#edabab", parents=[self] 650 ) 651 return ug 652 653 def extract_cells_by_id(self, idlist, use_point_ids=False) -> "UnstructuredGrid": 654 """Return a new `UnstructuredGrid` composed of the specified subset of indices.""" 655 selection_node = vtki.new("SelectionNode") 656 if use_point_ids: 657 selection_node.SetFieldType(vtki.get_class("SelectionNode").POINT) 658 contcells = vtki.get_class("SelectionNode").CONTAINING_CELLS() 659 selection_node.GetProperties().Set(contcells, 1) 660 else: 661 selection_node.SetFieldType(vtki.get_class("SelectionNode").CELL) 662 selection_node.SetContentType(vtki.get_class("SelectionNode").INDICES) 663 vidlist = utils.numpy2vtk(idlist, dtype="id") 664 selection_node.SetSelectionList(vidlist) 665 selection = vtki.new("Selection") 666 selection.AddNode(selection_node) 667 es = vtki.new("ExtractSelection") 668 es.SetInputData(0, self) 669 es.SetInputData(1, selection) 670 es.Update() 671 672 ug = UnstructuredGrid(es.GetOutput()) 673 pr = vtki.vtkProperty() 674 pr.DeepCopy(self.properties) 675 ug.actor.SetProperty(pr) 676 ug.properties = pr 677 678 ug.mapper.SetLookupTable(utils.ctf2lut(self)) 679 ug.pipeline = utils.OperationNode( 680 "extract_cells_by_id", 681 parents=[self], 682 comment=f"#cells {self.dataset.GetNumberOfCells()}", 683 c="#9e2a2b", 684 ) 685 return ug 686 687 def find_cell(self, p: list) -> int: 688 """Locate the cell that contains a point and return the cell ID.""" 689 cell = vtki.vtkTetra() 690 cell_id = vtki.mutable(0) 691 tol2 = vtki.mutable(0) 692 sub_id = vtki.mutable(0) 693 pcoords = [0, 0, 0] 694 weights = [0, 0, 0] 695 cid = self.dataset.FindCell(p, cell, cell_id, tol2, sub_id, pcoords, weights) 696 return cid 697 698 def clean(self) -> Self: 699 """ 700 Cleanup unused points and empty cells 701 """ 702 cl = vtki.new("StaticCleanUnstructuredGrid") 703 cl.SetInputData(self.dataset) 704 cl.RemoveUnusedPointsOn() 705 cl.ProduceMergeMapOff() 706 cl.AveragePointDataOff() 707 cl.Update() 708 709 self._update(cl.GetOutput()) 710 self.pipeline = utils.OperationNode( 711 "clean", 712 parents=[self], 713 comment=f"#cells {self.dataset.GetNumberOfCells()}", 714 c="#9e2a2b", 715 ) 716 return self 717 718 def extract_cells_on_plane(self, origin: tuple, normal: tuple) -> Self: 719 """ 720 Extract cells that are lying of the specified surface. 721 """ 722 bf = vtki.new("3DLinearGridCrinkleExtractor") 723 bf.SetInputData(self.dataset) 724 bf.CopyPointDataOn() 725 bf.CopyCellDataOn() 726 bf.RemoveUnusedPointsOff() 727 728 plane = vtki.new("Plane") 729 plane.SetOrigin(origin) 730 plane.SetNormal(normal) 731 bf.SetImplicitFunction(plane) 732 bf.Update() 733 734 self._update(bf.GetOutput(), reset_locators=False) 735 self.pipeline = utils.OperationNode( 736 "extract_cells_on_plane", 737 parents=[self], 738 comment=f"#cells {self.dataset.GetNumberOfCells()}", 739 c="#9e2a2b", 740 ) 741 return self 742 743 def extract_cells_on_sphere(self, center: tuple, radius: tuple) -> Self: 744 """ 745 Extract cells that are lying of the specified surface. 746 """ 747 bf = vtki.new("3DLinearGridCrinkleExtractor") 748 bf.SetInputData(self.dataset) 749 bf.CopyPointDataOn() 750 bf.CopyCellDataOn() 751 bf.RemoveUnusedPointsOff() 752 753 sph = vtki.new("Sphere") 754 sph.SetRadius(radius) 755 sph.SetCenter(center) 756 bf.SetImplicitFunction(sph) 757 bf.Update() 758 759 self._update(bf.GetOutput()) 760 self.pipeline = utils.OperationNode( 761 "extract_cells_on_sphere", 762 parents=[self], 763 comment=f"#cells {self.dataset.GetNumberOfCells()}", 764 c="#9e2a2b", 765 ) 766 return self 767 768 def extract_cells_on_cylinder(self, center: tuple, axis: tuple, radius: float) -> Self: 769 """ 770 Extract cells that are lying of the specified surface. 771 """ 772 bf = vtki.new("3DLinearGridCrinkleExtractor") 773 bf.SetInputData(self.dataset) 774 bf.CopyPointDataOn() 775 bf.CopyCellDataOn() 776 bf.RemoveUnusedPointsOff() 777 778 cyl = vtki.new("Cylinder") 779 cyl.SetRadius(radius) 780 cyl.SetCenter(center) 781 cyl.SetAxis(axis) 782 bf.SetImplicitFunction(cyl) 783 bf.Update() 784 785 self.pipeline = utils.OperationNode( 786 "extract_cells_on_cylinder", 787 parents=[self], 788 comment=f"#cells {self.dataset.GetNumberOfCells()}", 789 c="#9e2a2b", 790 ) 791 self._update(bf.GetOutput()) 792 return self 793 794 def cut_with_plane(self, origin=(0, 0, 0), normal="x") -> "UnstructuredGrid": 795 """ 796 Cut the object with the plane defined by a point and a normal. 797 798 Arguments: 799 origin : (list) 800 the cutting plane goes through this point 801 normal : (list, str) 802 normal vector to the cutting plane 803 """ 804 # if isinstance(self, vedo.Volume): 805 # raise RuntimeError("cut_with_plane() is not applicable to Volume objects.") 806 807 strn = str(normal) 808 if strn == "x": normal = (1, 0, 0) 809 elif strn == "y": normal = (0, 1, 0) 810 elif strn == "z": normal = (0, 0, 1) 811 elif strn == "-x": normal = (-1, 0, 0) 812 elif strn == "-y": normal = (0, -1, 0) 813 elif strn == "-z": normal = (0, 0, -1) 814 plane = vtki.new("Plane") 815 plane.SetOrigin(origin) 816 plane.SetNormal(normal) 817 clipper = vtki.new("ClipDataSet") 818 clipper.SetInputData(self.dataset) 819 clipper.SetClipFunction(plane) 820 clipper.GenerateClipScalarsOff() 821 clipper.GenerateClippedOutputOff() 822 clipper.SetValue(0) 823 clipper.Update() 824 cout = clipper.GetOutput() 825 826 if isinstance(cout, vtki.vtkUnstructuredGrid): 827 ug = vedo.UnstructuredGrid(cout) 828 if isinstance(self, vedo.UnstructuredGrid): 829 self._update(cout) 830 self.pipeline = utils.OperationNode("cut_with_plane", parents=[self], c="#9e2a2b") 831 return self 832 ug.pipeline = utils.OperationNode("cut_with_plane", parents=[self], c="#9e2a2b") 833 return ug 834 835 else: 836 self._update(cout) 837 self.pipeline = utils.OperationNode("cut_with_plane", parents=[self], c="#9e2a2b") 838 return self 839 840 def cut_with_box(self, box: Any) -> "UnstructuredGrid": 841 """ 842 Cut the grid with the specified bounding box. 843 844 Parameter box has format [xmin, xmax, ymin, ymax, zmin, zmax]. 845 If an object is passed, its bounding box are used. 846 847 This method always returns a TetMesh object. 848 849 Example: 850 ```python 851 from vedo import * 852 tmesh = TetMesh(dataurl+'limb_ugrid.vtk') 853 tmesh.color('rainbow') 854 cu = Cube(side=500).x(500) # any Mesh works 855 tmesh.cut_with_box(cu).show(axes=1) 856 ``` 857 858  859 """ 860 bc = vtki.new("BoxClipDataSet") 861 bc.SetInputData(self.dataset) 862 try: 863 boxb = box.bounds() 864 except AttributeError: 865 boxb = box 866 867 bc.SetBoxClip(*boxb) 868 bc.Update() 869 cout = bc.GetOutput() 870 871 # output of vtkBoxClipDataSet is always tetrahedrons 872 tm = vedo.TetMesh(cout) 873 tm.pipeline = utils.OperationNode("cut_with_box", parents=[self], c="#9e2a2b") 874 return tm 875 876 def cut_with_mesh(self, mesh: Mesh, invert=False, whole_cells=False, on_boundary=False) -> "UnstructuredGrid": 877 """ 878 Cut a `UnstructuredGrid` or `TetMesh` with a `Mesh`. 879 880 Use `invert` to return cut off part of the input object. 881 """ 882 ug = self.dataset 883 884 ippd = vtki.new("ImplicitPolyDataDistance") 885 ippd.SetInput(mesh.dataset) 886 887 if whole_cells or on_boundary: 888 clipper = vtki.new("ExtractGeometry") 889 clipper.SetInputData(ug) 890 clipper.SetImplicitFunction(ippd) 891 clipper.SetExtractInside(not invert) 892 clipper.SetExtractBoundaryCells(False) 893 if on_boundary: 894 clipper.SetExtractBoundaryCells(True) 895 clipper.SetExtractOnlyBoundaryCells(True) 896 else: 897 signed_dists = vtki.vtkFloatArray() 898 signed_dists.SetNumberOfComponents(1) 899 signed_dists.SetName("SignedDistance") 900 for pointId in range(ug.GetNumberOfPoints()): 901 p = ug.GetPoint(pointId) 902 signed_dist = ippd.EvaluateFunction(p) 903 signed_dists.InsertNextValue(signed_dist) 904 ug.GetPointData().AddArray(signed_dists) 905 ug.GetPointData().SetActiveScalars("SignedDistance") # NEEDED 906 clipper = vtki.new("ClipDataSet") 907 clipper.SetInputData(ug) 908 clipper.SetInsideOut(not invert) 909 clipper.SetValue(0.0) 910 911 clipper.Update() 912 913 out = vedo.UnstructuredGrid(clipper.GetOutput()) 914 out.pipeline = utils.OperationNode("cut_with_mesh", parents=[self], c="#9e2a2b") 915 return out
Support for UnstructuredGrid objects.
UnstructuredGrid(inputobj=None)
41 def __init__(self, inputobj=None): 42 """ 43 Support for UnstructuredGrid objects. 44 45 Arguments: 46 inputobj : (list, vtkUnstructuredGrid, str) 47 A list in the form `[points, cells, celltypes]`, 48 or a vtkUnstructuredGrid object, or a filename 49 50 Celltypes are identified by the following 51 [convention](https://vtk.org/doc/nightly/html/vtkCellType_8h_source.html). 52 """ 53 super().__init__() 54 55 self.dataset = None 56 57 self.mapper = vtki.new("PolyDataMapper") 58 self._actor = vtki.vtkActor() 59 self._actor.retrieve_object = weak_ref_to(self) 60 self._actor.SetMapper(self.mapper) 61 self.properties = self._actor.GetProperty() 62 63 self.transform = LinearTransform() 64 self.point_locator = None 65 self.cell_locator = None 66 self.line_locator = None 67 68 self.name = "UnstructuredGrid" 69 self.filename = "" 70 self.file_size = "" 71 72 self.info = {} 73 self.time = time.time() 74 self.rendered_at = set() 75 76 ###################################### 77 inputtype = str(type(inputobj)) 78 79 if inputobj is None: 80 self.dataset = vtki.vtkUnstructuredGrid() 81 82 elif utils.is_sequence(inputobj): 83 84 pts, cells, celltypes = inputobj 85 assert len(cells) == len(celltypes) 86 87 self.dataset = vtki.vtkUnstructuredGrid() 88 89 if not utils.is_sequence(cells[0]): 90 tets = [] 91 nf = cells[0] + 1 92 for i, cl in enumerate(cells): 93 if i in (nf, 0): 94 k = i + 1 95 nf = cl + k 96 cell = [cells[j + k] for j in range(cl)] 97 tets.append(cell) 98 cells = tets 99 100 # This would fill the points and use those to define orientation 101 vpts = utils.numpy2vtk(pts, dtype=np.float32) 102 points = vtki.vtkPoints() 103 points.SetData(vpts) 104 self.dataset.SetPoints(points) 105 106 # Fill cells 107 # https://vtk.org/doc/nightly/html/vtkCellType_8h_source.html 108 for i, ct in enumerate(celltypes): 109 if ct == vtki.cell_types["VERTEX"]: 110 cell = vtki.vtkVertex() 111 elif ct == vtki.cell_types["POLY_VERTEX"]: 112 cell = vtki.vtkPolyVertex() 113 elif ct == vtki.cell_types["TETRA"]: 114 cell = vtki.vtkTetra() 115 elif ct == vtki.cell_types["WEDGE"]: 116 cell = vtki.vtkWedge() 117 elif ct == vtki.cell_types["LINE"]: 118 cell = vtki.vtkLine() 119 elif ct == vtki.cell_types["POLY_LINE"]: 120 cell = vtki.vtkPolyLine() 121 elif ct == vtki.cell_types["TRIANGLE"]: 122 cell = vtki.vtkTriangle() 123 elif ct == vtki.cell_types["TRIANGLE_STRIP"]: 124 cell = vtki.vtkTriangleStrip() 125 elif ct == vtki.cell_types["POLYGON"]: 126 cell = vtki.vtkPolygon() 127 elif ct == vtki.cell_types["PIXEL"]: 128 cell = vtki.vtkPixel() 129 elif ct == vtki.cell_types["QUAD"]: 130 cell = vtki.vtkQuad() 131 elif ct == vtki.cell_types["VOXEL"]: 132 cell = vtki.vtkVoxel() 133 elif ct == vtki.cell_types["PYRAMID"]: 134 cell = vtki.vtkPyramid() 135 elif ct == vtki.cell_types["HEXAHEDRON"]: 136 cell = vtki.vtkHexahedron() 137 elif ct == vtki.cell_types["HEXAGONAL_PRISM"]: 138 cell = vtki.vtkHexagonalPrism() 139 elif ct == vtki.cell_types["PENTAGONAL_PRISM"]: 140 cell = vtki.vtkPentagonalPrism() 141 elif ct == vtki.cell_types["QUADRATIC_TETRA"]: 142 from vtkmodules.vtkCommonDataModel import vtkQuadraticTetra 143 cell = vtkQuadraticTetra() 144 elif ct == vtki.cell_types["QUADRATIC_HEXAHEDRON"]: 145 from vtkmodules.vtkCommonDataModel import vtkQuadraticHexahedron 146 cell = vtkQuadraticHexahedron() 147 elif ct == vtki.cell_types["QUADRATIC_WEDGE"]: 148 from vtkmodules.vtkCommonDataModel import vtkQuadraticWedge 149 cell = vtkQuadraticWedge() 150 elif ct == vtki.cell_types["QUADRATIC_PYRAMID"]: 151 from vtkmodules.vtkCommonDataModel import vtkQuadraticPyramid 152 cell = vtkQuadraticPyramid() 153 elif ct == vtki.cell_types["QUADRATIC_LINEAR_QUAD"]: 154 from vtkmodules.vtkCommonDataModel import vtkQuadraticLinearQuad 155 cell = vtkQuadraticLinearQuad() 156 elif ct == vtki.cell_types["QUADRATIC_LINEAR_WEDGE"]: 157 from vtkmodules.vtkCommonDataModel import vtkQuadraticLinearWedge 158 cell = vtkQuadraticLinearWedge() 159 elif ct == vtki.cell_types["BIQUADRATIC_QUADRATIC_WEDGE"]: 160 from vtkmodules.vtkCommonDataModel import vtkBiQuadraticQuadraticWedge 161 cell = vtkBiQuadraticQuadraticWedge() 162 elif ct == vtki.cell_types["BIQUADRATIC_QUADRATIC_HEXAHEDRON"]: 163 from vtkmodules.vtkCommonDataModel import vtkBiQuadraticQuadraticHexahedron 164 cell = vtkBiQuadraticQuadraticHexahedron() 165 elif ct == vtki.cell_types["BIQUADRATIC_TRIANGLE"]: 166 from vtkmodules.vtkCommonDataModel import vtkBiQuadraticTriangle 167 cell = vtkBiQuadraticTriangle() 168 elif ct == vtki.cell_types["CUBIC_LINE"]: 169 from vtkmodules.vtkCommonDataModel import vtkCubicLine 170 cell = vtkCubicLine() 171 elif ct == vtki.cell_types["CONVEX_POINT_SET"]: 172 from vtkmodules.vtkCommonDataModel import vtkConvexPointSet 173 cell = vtkConvexPointSet() 174 elif ct == vtki.cell_types["POLYHEDRON"]: 175 from vtkmodules.vtkCommonDataModel import vtkPolyhedron 176 cell = vtkPolyhedron() 177 elif ct == vtki.cell_types["HIGHER_ORDER_TRIANGLE"]: 178 from vtkmodules.vtkCommonDataModel import vtkHigherOrderTriangle 179 cell = vtkHigherOrderTriangle() 180 elif ct == vtki.cell_types["HIGHER_ORDER_QUAD"]: 181 from vtkmodules.vtkCommonDataModel import vtkHigherOrderQuadrilateral 182 cell = vtkHigherOrderQuadrilateral() 183 elif ct == vtki.cell_types["HIGHER_ORDER_TETRAHEDRON"]: 184 from vtkmodules.vtkCommonDataModel import vtkHigherOrderTetra 185 cell = vtkHigherOrderTetra() 186 elif ct == vtki.cell_types["HIGHER_ORDER_WEDGE"]: 187 from vtkmodules.vtkCommonDataModel import vtkHigherOrderWedge 188 cell = vtkHigherOrderWedge() 189 elif ct == vtki.cell_types["HIGHER_ORDER_HEXAHEDRON"]: 190 from vtkmodules.vtkCommonDataModel import vtkHigherOrderHexahedron 191 cell = vtkHigherOrderHexahedron() 192 elif ct == vtki.cell_types["LAGRANGE_CURVE"]: 193 from vtkmodules.vtkCommonDataModel import vtkLagrangeCurve 194 cell = vtkLagrangeCurve() 195 elif ct == vtki.cell_types["LAGRANGE_TRIANGLE"]: 196 from vtkmodules.vtkCommonDataModel import vtkLagrangeTriangle 197 cell = vtkLagrangeTriangle() 198 elif ct == vtki.cell_types["LAGRANGE_QUADRILATERAL"]: 199 from vtkmodules.vtkCommonDataModel import vtkLagrangeQuadrilateral 200 cell = vtkLagrangeQuadrilateral() 201 elif ct == vtki.cell_types["LAGRANGE_TETRAHEDRON"]: 202 from vtkmodules.vtkCommonDataModel import vtkLagrangeTetra 203 cell = vtkLagrangeTetra() 204 elif ct == vtki.cell_types["LAGRANGE_HEXAHEDRON"]: 205 from vtkmodules.vtkCommonDataModel import vtkLagrangeHexahedron 206 cell = vtkLagrangeHexahedron() 207 elif ct == vtki.cell_types["LAGRANGE_WEDGE"]: 208 from vtkmodules.vtkCommonDataModel import vtkLagrangeWedge 209 cell = vtkLagrangeWedge() 210 elif ct == vtki.cell_types["BEZIER_CURVE"]: 211 from vtkmodules.vtkCommonDataModel import vtkBezierCurve 212 cell = vtkBezierCurve() 213 elif ct == vtki.cell_types["BEZIER_TRIANGLE"]: 214 from vtkmodules.vtkCommonDataModel import vtkBezierTriangle 215 cell = vtkBezierTriangle() 216 elif ct == vtki.cell_types["BEZIER_QUADRILATERAL"]: 217 from vtkmodules.vtkCommonDataModel import vtkBezierQuadrilateral 218 cell = vtkBezierQuadrilateral() 219 elif ct == vtki.cell_types["BEZIER_TETRAHEDRON"]: 220 from vtkmodules.vtkCommonDataModel import vtkBezierTetra 221 cell = vtkBezierTetra() 222 elif ct == vtki.cell_types["BEZIER_HEXAHEDRON"]: 223 from vtkmodules.vtkCommonDataModel import vtkBezierHexahedron 224 cell = vtkBezierHexahedron() 225 elif ct == vtki.cell_types["BEZIER_WEDGE"]: 226 from vtkmodules.vtkCommonDataModel import vtkBezierWedge 227 cell = vtkBezierWedge() 228 else: 229 vedo.logger.error( 230 f"UnstructuredGrid: cell type {ct} not supported. Skip.") 231 continue 232 233 cpids = cell.GetPointIds() 234 cell_conn = cells[i] 235 for j, pid in enumerate(cell_conn): 236 cpids.SetId(j, pid) 237 self.dataset.InsertNextCell(ct, cpids) 238 239 elif "UnstructuredGrid" in inputtype: 240 self.dataset = inputobj 241 242 elif isinstance(inputobj, str): 243 if "https://" in inputobj: 244 inputobj = download(inputobj, verbose=False) 245 self.filename = inputobj 246 if inputobj.endswith(".vtu"): 247 reader = vtki.new("XMLUnstructuredGridReader") 248 else: 249 reader = vtki.new("UnstructuredGridReader") 250 self.filename = inputobj 251 reader.SetFileName(inputobj) 252 reader.Update() 253 self.dataset = reader.GetOutput() 254 255 else: 256 # this converts other types of vtk objects to UnstructuredGrid 257 apf = vtki.new("AppendFilter") 258 try: 259 apf.AddInputData(inputobj) 260 except TypeError: 261 apf.AddInputData(inputobj.dataset) 262 apf.Update() 263 self.dataset = apf.GetOutput() 264 265 self.properties.SetColor(0.89, 0.455, 0.671) # pink7 266 267 self.pipeline = utils.OperationNode( 268 self, comment=f"#cells {self.dataset.GetNumberOfCells()}", c="#4cc9f0" 269 )
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.
actor
412 @property 413 def actor(self): 414 """Return the `vtkActor` of the object.""" 415 # print("building actor") 416 gf = vtki.new("GeometryFilter") 417 gf.SetInputData(self.dataset) 418 gf.Update() 419 out = gf.GetOutput() 420 self.mapper.SetInputData(out) 421 self.mapper.Modified() 422 return self._actor
Return the vtkActor of the object.
def
merge(self, *others) -> Self:
435 def merge(self, *others) -> Self: 436 """ 437 Merge multiple datasets into one single `UnstrcturedGrid`. 438 """ 439 apf = vtki.new("AppendFilter") 440 for o in others: 441 if isinstance(o, UnstructuredGrid): 442 apf.AddInputData(o.dataset) 443 elif isinstance(o, vtki.vtkUnstructuredGrid): 444 apf.AddInputData(o) 445 else: 446 vedo.printc("Error: cannot merge type", type(o), c="r") 447 apf.Update() 448 self._update(apf.GetOutput()) 449 self.pipeline = utils.OperationNode( 450 "merge", parents=[self, *others], c="#9e2a2b" 451 ) 452 return self
Merge multiple datasets into one single UnstrcturedGrid.
454 def copy(self, deep=True) -> "UnstructuredGrid": 455 """Return a copy of the object. Alias of `clone()`.""" 456 return self.clone(deep=deep)
Return a copy of the object. Alias of clone().
458 def clone(self, deep=True) -> "UnstructuredGrid": 459 """Clone the UnstructuredGrid object to yield an exact copy.""" 460 ug = vtki.vtkUnstructuredGrid() 461 if deep: 462 ug.DeepCopy(self.dataset) 463 else: 464 ug.ShallowCopy(self.dataset) 465 if isinstance(self, vedo.UnstructuredGrid): 466 cloned = vedo.UnstructuredGrid(ug) 467 else: 468 cloned = vedo.TetMesh(ug) 469 470 cloned.copy_properties_from(self) 471 472 cloned.pipeline = utils.OperationNode( 473 "clone", parents=[self], shape="diamond", c="#bbe1ed" 474 ) 475 return cloned
Clone the UnstructuredGrid object to yield an exact copy.
def
bounds(self) -> numpy.ndarray:
477 def bounds(self) -> np.ndarray: 478 """ 479 Get the object bounds. 480 Returns a list in format `[xmin,xmax, ymin,ymax, zmin,zmax]`. 481 """ 482 # OVERRIDE CommonAlgorithms.bounds() which is too slow 483 return np.array(self.dataset.GetBounds())
Get the object bounds.
Returns a list in format [xmin,xmax, ymin,ymax, zmin,zmax].
def
threshold(self, name=None, above=None, below=None, on='cells') -> Self:
485 def threshold(self, name=None, above=None, below=None, on="cells") -> Self: 486 """ 487 Threshold the tetrahedral mesh by a cell scalar value. 488 Reduce to only tets which satisfy the threshold limits. 489 490 - if `above = below` will only select tets with that specific value. 491 - if `above > below` selection range is flipped. 492 493 Set keyword "on" to either "cells" or "points". 494 """ 495 th = vtki.new("Threshold") 496 th.SetInputData(self.dataset) 497 498 if name is None: 499 if self.celldata.keys(): 500 name = self.celldata.keys()[0] 501 th.SetInputArrayToProcess(0, 0, 0, 1, name) 502 elif self.pointdata.keys(): 503 name = self.pointdata.keys()[0] 504 th.SetInputArrayToProcess(0, 0, 0, 0, name) 505 if name is None: 506 vedo.logger.warning("cannot find active array. Skip.") 507 return self 508 else: 509 if on.startswith("c"): 510 th.SetInputArrayToProcess(0, 0, 0, 1, name) 511 else: 512 th.SetInputArrayToProcess(0, 0, 0, 0, name) 513 514 if above is not None: 515 th.SetLowerThreshold(above) 516 517 if below is not None: 518 th.SetUpperThreshold(below) 519 520 th.Update() 521 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 = belowwill only select tets with that specific value. - if
above > belowselection range is flipped.
Set keyword "on" to either "cells" or "points".
523 def isosurface(self, value=None, flying_edges=False) -> "vedo.mesh.Mesh": 524 """ 525 Return an `Mesh` isosurface extracted from the `Volume` object. 526 527 Set `value` as single float or list of values to draw the isosurface(s). 528 Use flying_edges for faster results (but sometimes can interfere with `smooth()`). 529 530 Examples: 531 - [isosurfaces1.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/isosurfaces1.py) 532 533  534 """ 535 scrange = self.dataset.GetScalarRange() 536 537 if flying_edges: 538 cf = vtki.new("FlyingEdges3D") 539 cf.InterpolateAttributesOn() 540 else: 541 cf = vtki.new("ContourFilter") 542 cf.UseScalarTreeOn() 543 544 cf.SetInputData(self.dataset) 545 cf.ComputeNormalsOn() 546 547 if utils.is_sequence(value): 548 cf.SetNumberOfContours(len(value)) 549 for i, t in enumerate(value): 550 cf.SetValue(i, t) 551 else: 552 if value is None: 553 value = (2 * scrange[0] + scrange[1]) / 3.0 554 # print("automatic isosurface value =", value) 555 cf.SetValue(0, value) 556 557 cf.Update() 558 poly = cf.GetOutput() 559 560 out = vedo.mesh.Mesh(poly, c=None).flat() 561 out.mapper.SetScalarRange(scrange[0], scrange[1]) 562 563 out.pipeline = utils.OperationNode( 564 "isosurface", 565 parents=[self], 566 comment=f"#pts {out.dataset.GetNumberOfPoints()}", 567 c="#4cc9f0:#e9c46a", 568 ) 569 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:
def
shrink(self, fraction=0.8) -> Self:
571 def shrink(self, fraction=0.8) -> Self: 572 """ 573 Shrink the individual cells. 574 575  576 """ 577 sf = vtki.new("ShrinkFilter") 578 sf.SetInputData(self.dataset) 579 sf.SetShrinkFactor(fraction) 580 sf.Update() 581 out = sf.GetOutput() 582 self._update(out) 583 self.pipeline = utils.OperationNode( 584 "shrink", comment=f"by {fraction}", parents=[self], c="#9e2a2b" 585 ) 586 return self
Shrink the individual cells.
588 def tomesh(self, fill=False, shrink=1.0) -> "vedo.mesh.Mesh": 589 """ 590 Build a polygonal `Mesh` from the current object. 591 592 If `fill=True`, the interior faces of all the cells are created. 593 (setting a `shrink` value slightly smaller than the default 1.0 594 can avoid flickering due to internal adjacent faces). 595 596 If `fill=False`, only the boundary faces will be generated. 597 """ 598 gf = vtki.new("GeometryFilter") 599 if fill: 600 sf = vtki.new("ShrinkFilter") 601 sf.SetInputData(self.dataset) 602 sf.SetShrinkFactor(shrink) 603 sf.Update() 604 gf.SetInputData(sf.GetOutput()) 605 gf.Update() 606 poly = gf.GetOutput() 607 else: 608 gf.SetInputData(self.dataset) 609 gf.Update() 610 poly = gf.GetOutput() 611 612 msh = vedo.mesh.Mesh(poly) 613 msh.copy_properties_from(self) 614 msh.pipeline = utils.OperationNode( 615 "tomesh", parents=[self], comment=f"fill={fill}", c="#9e2a2b:#e9c46a" 616 ) 617 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.
cell_types_array
619 @property 620 def cell_types_array(self): 621 """Return the list of cell types in the dataset.""" 622 uarr = self.dataset.GetCellTypesArray() 623 return utils.vtk2numpy(uarr)
Return the list of cell types in the dataset.
625 def extract_cells_by_type(self, ctype) -> "UnstructuredGrid": 626 """Extract a specific cell type and return a new `UnstructuredGrid`.""" 627 if isinstance(ctype, str): 628 try: 629 ctype = vtki.cell_types[ctype.upper()] 630 except KeyError: 631 vedo.logger.error(f"extract_cells_by_type: cell type {ctype} does not exist. Skip.") 632 return self 633 uarr = self.dataset.GetCellTypesArray() 634 ctarrtyp = np.where(utils.vtk2numpy(uarr) == ctype)[0] 635 uarrtyp = utils.numpy2vtk(ctarrtyp, deep=False, dtype="id") 636 selection_node = vtki.new("SelectionNode") 637 selection_node.SetFieldType(vtki.get_class("SelectionNode").CELL) 638 selection_node.SetContentType(vtki.get_class("SelectionNode").INDICES) 639 selection_node.SetSelectionList(uarrtyp) 640 selection = vtki.new("Selection") 641 selection.AddNode(selection_node) 642 es = vtki.new("ExtractSelection") 643 es.SetInputData(0, self.dataset) 644 es.SetInputData(1, selection) 645 es.Update() 646 647 ug = UnstructuredGrid(es.GetOutput()) 648 ug.pipeline = utils.OperationNode( 649 "extract_cell_type", comment=f"type {ctype}", c="#edabab", parents=[self] 650 ) 651 return ug
Extract a specific cell type and return a new UnstructuredGrid.
653 def extract_cells_by_id(self, idlist, use_point_ids=False) -> "UnstructuredGrid": 654 """Return a new `UnstructuredGrid` composed of the specified subset of indices.""" 655 selection_node = vtki.new("SelectionNode") 656 if use_point_ids: 657 selection_node.SetFieldType(vtki.get_class("SelectionNode").POINT) 658 contcells = vtki.get_class("SelectionNode").CONTAINING_CELLS() 659 selection_node.GetProperties().Set(contcells, 1) 660 else: 661 selection_node.SetFieldType(vtki.get_class("SelectionNode").CELL) 662 selection_node.SetContentType(vtki.get_class("SelectionNode").INDICES) 663 vidlist = utils.numpy2vtk(idlist, dtype="id") 664 selection_node.SetSelectionList(vidlist) 665 selection = vtki.new("Selection") 666 selection.AddNode(selection_node) 667 es = vtki.new("ExtractSelection") 668 es.SetInputData(0, self) 669 es.SetInputData(1, selection) 670 es.Update() 671 672 ug = UnstructuredGrid(es.GetOutput()) 673 pr = vtki.vtkProperty() 674 pr.DeepCopy(self.properties) 675 ug.actor.SetProperty(pr) 676 ug.properties = pr 677 678 ug.mapper.SetLookupTable(utils.ctf2lut(self)) 679 ug.pipeline = utils.OperationNode( 680 "extract_cells_by_id", 681 parents=[self], 682 comment=f"#cells {self.dataset.GetNumberOfCells()}", 683 c="#9e2a2b", 684 ) 685 return ug
Return a new UnstructuredGrid composed of the specified subset of indices.
def
find_cell(self, p: list) -> int:
687 def find_cell(self, p: list) -> int: 688 """Locate the cell that contains a point and return the cell ID.""" 689 cell = vtki.vtkTetra() 690 cell_id = vtki.mutable(0) 691 tol2 = vtki.mutable(0) 692 sub_id = vtki.mutable(0) 693 pcoords = [0, 0, 0] 694 weights = [0, 0, 0] 695 cid = self.dataset.FindCell(p, cell, cell_id, tol2, sub_id, pcoords, weights) 696 return cid
Locate the cell that contains a point and return the cell ID.
def
clean(self) -> Self:
698 def clean(self) -> Self: 699 """ 700 Cleanup unused points and empty cells 701 """ 702 cl = vtki.new("StaticCleanUnstructuredGrid") 703 cl.SetInputData(self.dataset) 704 cl.RemoveUnusedPointsOn() 705 cl.ProduceMergeMapOff() 706 cl.AveragePointDataOff() 707 cl.Update() 708 709 self._update(cl.GetOutput()) 710 self.pipeline = utils.OperationNode( 711 "clean", 712 parents=[self], 713 comment=f"#cells {self.dataset.GetNumberOfCells()}", 714 c="#9e2a2b", 715 ) 716 return self
Cleanup unused points and empty cells
def
extract_cells_on_plane(self, origin: tuple, normal: tuple) -> Self:
718 def extract_cells_on_plane(self, origin: tuple, normal: tuple) -> Self: 719 """ 720 Extract cells that are lying of the specified surface. 721 """ 722 bf = vtki.new("3DLinearGridCrinkleExtractor") 723 bf.SetInputData(self.dataset) 724 bf.CopyPointDataOn() 725 bf.CopyCellDataOn() 726 bf.RemoveUnusedPointsOff() 727 728 plane = vtki.new("Plane") 729 plane.SetOrigin(origin) 730 plane.SetNormal(normal) 731 bf.SetImplicitFunction(plane) 732 bf.Update() 733 734 self._update(bf.GetOutput(), reset_locators=False) 735 self.pipeline = utils.OperationNode( 736 "extract_cells_on_plane", 737 parents=[self], 738 comment=f"#cells {self.dataset.GetNumberOfCells()}", 739 c="#9e2a2b", 740 ) 741 return self
Extract cells that are lying of the specified surface.
def
extract_cells_on_sphere(self, center: tuple, radius: tuple) -> Self:
743 def extract_cells_on_sphere(self, center: tuple, radius: tuple) -> Self: 744 """ 745 Extract cells that are lying of the specified surface. 746 """ 747 bf = vtki.new("3DLinearGridCrinkleExtractor") 748 bf.SetInputData(self.dataset) 749 bf.CopyPointDataOn() 750 bf.CopyCellDataOn() 751 bf.RemoveUnusedPointsOff() 752 753 sph = vtki.new("Sphere") 754 sph.SetRadius(radius) 755 sph.SetCenter(center) 756 bf.SetImplicitFunction(sph) 757 bf.Update() 758 759 self._update(bf.GetOutput()) 760 self.pipeline = utils.OperationNode( 761 "extract_cells_on_sphere", 762 parents=[self], 763 comment=f"#cells {self.dataset.GetNumberOfCells()}", 764 c="#9e2a2b", 765 ) 766 return self
Extract cells that are lying of the specified surface.
def
extract_cells_on_cylinder(self, center: tuple, axis: tuple, radius: float) -> Self:
768 def extract_cells_on_cylinder(self, center: tuple, axis: tuple, radius: float) -> Self: 769 """ 770 Extract cells that are lying of the specified surface. 771 """ 772 bf = vtki.new("3DLinearGridCrinkleExtractor") 773 bf.SetInputData(self.dataset) 774 bf.CopyPointDataOn() 775 bf.CopyCellDataOn() 776 bf.RemoveUnusedPointsOff() 777 778 cyl = vtki.new("Cylinder") 779 cyl.SetRadius(radius) 780 cyl.SetCenter(center) 781 cyl.SetAxis(axis) 782 bf.SetImplicitFunction(cyl) 783 bf.Update() 784 785 self.pipeline = utils.OperationNode( 786 "extract_cells_on_cylinder", 787 parents=[self], 788 comment=f"#cells {self.dataset.GetNumberOfCells()}", 789 c="#9e2a2b", 790 ) 791 self._update(bf.GetOutput()) 792 return self
Extract cells that are lying of the specified surface.
794 def cut_with_plane(self, origin=(0, 0, 0), normal="x") -> "UnstructuredGrid": 795 """ 796 Cut the object with the plane defined by a point and a normal. 797 798 Arguments: 799 origin : (list) 800 the cutting plane goes through this point 801 normal : (list, str) 802 normal vector to the cutting plane 803 """ 804 # if isinstance(self, vedo.Volume): 805 # raise RuntimeError("cut_with_plane() is not applicable to Volume objects.") 806 807 strn = str(normal) 808 if strn == "x": normal = (1, 0, 0) 809 elif strn == "y": normal = (0, 1, 0) 810 elif strn == "z": normal = (0, 0, 1) 811 elif strn == "-x": normal = (-1, 0, 0) 812 elif strn == "-y": normal = (0, -1, 0) 813 elif strn == "-z": normal = (0, 0, -1) 814 plane = vtki.new("Plane") 815 plane.SetOrigin(origin) 816 plane.SetNormal(normal) 817 clipper = vtki.new("ClipDataSet") 818 clipper.SetInputData(self.dataset) 819 clipper.SetClipFunction(plane) 820 clipper.GenerateClipScalarsOff() 821 clipper.GenerateClippedOutputOff() 822 clipper.SetValue(0) 823 clipper.Update() 824 cout = clipper.GetOutput() 825 826 if isinstance(cout, vtki.vtkUnstructuredGrid): 827 ug = vedo.UnstructuredGrid(cout) 828 if isinstance(self, vedo.UnstructuredGrid): 829 self._update(cout) 830 self.pipeline = utils.OperationNode("cut_with_plane", parents=[self], c="#9e2a2b") 831 return self 832 ug.pipeline = utils.OperationNode("cut_with_plane", parents=[self], c="#9e2a2b") 833 return ug 834 835 else: 836 self._update(cout) 837 self.pipeline = utils.OperationNode("cut_with_plane", parents=[self], c="#9e2a2b") 838 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
840 def cut_with_box(self, box: Any) -> "UnstructuredGrid": 841 """ 842 Cut the grid with the specified bounding box. 843 844 Parameter box has format [xmin, xmax, ymin, ymax, zmin, zmax]. 845 If an object is passed, its bounding box are used. 846 847 This method always returns a TetMesh object. 848 849 Example: 850 ```python 851 from vedo import * 852 tmesh = TetMesh(dataurl+'limb_ugrid.vtk') 853 tmesh.color('rainbow') 854 cu = Cube(side=500).x(500) # any Mesh works 855 tmesh.cut_with_box(cu).show(axes=1) 856 ``` 857 858  859 """ 860 bc = vtki.new("BoxClipDataSet") 861 bc.SetInputData(self.dataset) 862 try: 863 boxb = box.bounds() 864 except AttributeError: 865 boxb = box 866 867 bc.SetBoxClip(*boxb) 868 bc.Update() 869 cout = bc.GetOutput() 870 871 # output of vtkBoxClipDataSet is always tetrahedrons 872 tm = vedo.TetMesh(cout) 873 tm.pipeline = utils.OperationNode("cut_with_box", parents=[self], c="#9e2a2b") 874 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 *
tmesh = TetMesh(dataurl+'limb_ugrid.vtk')
tmesh.color('rainbow')
cu = Cube(side=500).x(500) # any Mesh works
tmesh.cut_with_box(cu).show(axes=1)
def
cut_with_mesh( self, mesh: vedo.mesh.Mesh, invert=False, whole_cells=False, on_boundary=False) -> UnstructuredGrid:
876 def cut_with_mesh(self, mesh: Mesh, invert=False, whole_cells=False, on_boundary=False) -> "UnstructuredGrid": 877 """ 878 Cut a `UnstructuredGrid` or `TetMesh` with a `Mesh`. 879 880 Use `invert` to return cut off part of the input object. 881 """ 882 ug = self.dataset 883 884 ippd = vtki.new("ImplicitPolyDataDistance") 885 ippd.SetInput(mesh.dataset) 886 887 if whole_cells or on_boundary: 888 clipper = vtki.new("ExtractGeometry") 889 clipper.SetInputData(ug) 890 clipper.SetImplicitFunction(ippd) 891 clipper.SetExtractInside(not invert) 892 clipper.SetExtractBoundaryCells(False) 893 if on_boundary: 894 clipper.SetExtractBoundaryCells(True) 895 clipper.SetExtractOnlyBoundaryCells(True) 896 else: 897 signed_dists = vtki.vtkFloatArray() 898 signed_dists.SetNumberOfComponents(1) 899 signed_dists.SetName("SignedDistance") 900 for pointId in range(ug.GetNumberOfPoints()): 901 p = ug.GetPoint(pointId) 902 signed_dist = ippd.EvaluateFunction(p) 903 signed_dists.InsertNextValue(signed_dist) 904 ug.GetPointData().AddArray(signed_dists) 905 ug.GetPointData().SetActiveScalars("SignedDistance") # NEEDED 906 clipper = vtki.new("ClipDataSet") 907 clipper.SetInputData(ug) 908 clipper.SetInsideOut(not invert) 909 clipper.SetValue(0.0) 910 911 clipper.Update() 912 913 out = vedo.UnstructuredGrid(clipper.GetOutput()) 914 out.pipeline = utils.OperationNode("cut_with_mesh", parents=[self], c="#9e2a2b") 915 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
- update_dataset
- 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
- keep_cell_types
- 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
- generate_random_data
- integrate_data
- write
- signed_distance
- unsigned_distance
- smooth_data
- compute_streamlines
- 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
919class TetMesh(UnstructuredGrid): 920 """The class describing tetrahedral meshes.""" 921 922 def __init__(self, inputobj=None): 923 """ 924 Arguments: 925 inputobj : (vtkUnstructuredGrid, list, str, tetgenpy.TetgenIO) 926 list of points and tet indices, or filename 927 """ 928 super().__init__() 929 930 self.dataset = None 931 932 self.mapper = vtki.new("PolyDataMapper") 933 self._actor = vtki.vtkActor() 934 self._actor.retrieve_object = weak_ref_to(self) 935 self._actor.SetMapper(self.mapper) 936 self.properties = self._actor.GetProperty() 937 938 self.name = "TetMesh" 939 940 # print('TetMesh inputtype', type(inputobj)) 941 942 ################### 943 if inputobj is None: 944 self.dataset = vtki.vtkUnstructuredGrid() 945 946 elif isinstance(inputobj, vtki.vtkUnstructuredGrid): 947 self.dataset = inputobj 948 949 elif isinstance(inputobj, UnstructuredGrid): 950 self.dataset = inputobj.dataset 951 952 elif "TetgenIO" in str(type(inputobj)): # tetgenpy object 953 inputobj = [inputobj.points(), inputobj.tetrahedra()] 954 955 elif isinstance(inputobj, vtki.vtkRectilinearGrid): 956 r2t = vtki.new("RectilinearGridToTetrahedra") 957 r2t.SetInputData(inputobj) 958 r2t.RememberVoxelIdOn() 959 r2t.SetTetraPerCellTo6() 960 r2t.Update() 961 self.dataset = r2t.GetOutput() 962 963 elif isinstance(inputobj, vtki.vtkDataSet): 964 r2t = vtki.new("DataSetTriangleFilter") 965 r2t.SetInputData(inputobj) 966 r2t.TetrahedraOnlyOn() 967 r2t.Update() 968 self.dataset = r2t.GetOutput() 969 970 elif isinstance(inputobj, str): 971 if "https://" in inputobj: 972 inputobj = download(inputobj, verbose=False) 973 if inputobj.endswith(".vtu"): 974 reader = vtki.new("XMLUnstructuredGridReader") 975 else: 976 reader = vtki.new("UnstructuredGridReader") 977 978 if not os.path.isfile(inputobj): 979 # for some reason vtk Reader does not complain 980 vedo.logger.error(f"file {inputobj} not found") 981 raise FileNotFoundError 982 983 self.filename = inputobj 984 reader.SetFileName(inputobj) 985 reader.Update() 986 ug = reader.GetOutput() 987 988 tt = vtki.new("DataSetTriangleFilter") 989 tt.SetInputData(ug) 990 tt.SetTetrahedraOnly(True) 991 tt.Update() 992 self.dataset = tt.GetOutput() 993 994 ############################### 995 if utils.is_sequence(inputobj): 996 self.dataset = vtki.vtkUnstructuredGrid() 997 998 points, cells = inputobj 999 if len(points) == 0: 1000 return 1001 if not utils.is_sequence(points[0]): 1002 return 1003 if len(cells) == 0: 1004 return 1005 1006 if not utils.is_sequence(cells[0]): 1007 tets = [] 1008 nf = cells[0] + 1 1009 for i, cl in enumerate(cells): 1010 if i in (nf, 0): 1011 k = i + 1 1012 nf = cl + k 1013 cell = [cells[j + k] for j in range(cl)] 1014 tets.append(cell) 1015 cells = tets 1016 1017 source_points = vtki.vtkPoints() 1018 varr = utils.numpy2vtk(points, dtype=np.float32) 1019 source_points.SetData(varr) 1020 self.dataset.SetPoints(source_points) 1021 1022 source_tets = vtki.vtkCellArray() 1023 for f in cells: 1024 ele = vtki.vtkTetra() 1025 pid = ele.GetPointIds() 1026 for i, fi in enumerate(f): 1027 pid.SetId(i, fi) 1028 source_tets.InsertNextCell(ele) 1029 self.dataset.SetCells(vtki.cell_types["TETRA"], source_tets) 1030 1031 if not self.dataset: 1032 vedo.logger.error(f"cannot understand input type {type(inputobj)}") 1033 return 1034 1035 self.properties.SetColor(0.352, 0.612, 0.996) # blue7 1036 1037 self.pipeline = utils.OperationNode( 1038 self, comment=f"#tets {self.dataset.GetNumberOfCells()}", c="#9e2a2b" 1039 ) 1040 1041 ################################################################## 1042 def __str__(self): 1043 """Print a string summary of the `TetMesh` object.""" 1044 module = self.__class__.__module__ 1045 name = self.__class__.__name__ 1046 out = vedo.printc( 1047 f"{module}.{name} at ({hex(self.memory_address())})".ljust(75), 1048 c="c", bold=True, invert=True, return_string=True, 1049 ) 1050 out += "\x1b[0m\u001b[36m" 1051 1052 out += "nr. of verts".ljust(14) + ": " + str(self.npoints) + "\n" 1053 out += "nr. of tetras".ljust(14) + ": " + str(self.ncells) + "\n" 1054 1055 if self.npoints: 1056 out+="size".ljust(14)+ ": average=" + utils.precision(self.average_size(),6) 1057 out+=", diagonal="+ utils.precision(self.diagonal_size(), 6)+ "\n" 1058 out+="center of mass".ljust(14) + ": " + utils.precision(self.center_of_mass(),6)+"\n" 1059 1060 bnds = self.bounds() 1061 bx1, bx2 = utils.precision(bnds[0], 3), utils.precision(bnds[1], 3) 1062 by1, by2 = utils.precision(bnds[2], 3), utils.precision(bnds[3], 3) 1063 bz1, bz2 = utils.precision(bnds[4], 3), utils.precision(bnds[5], 3) 1064 out += "bounds".ljust(14) + ":" 1065 out += " x=(" + bx1 + ", " + bx2 + ")," 1066 out += " y=(" + by1 + ", " + by2 + ")," 1067 out += " z=(" + bz1 + ", " + bz2 + ")\n" 1068 1069 for key in self.pointdata.keys(): 1070 arr = self.pointdata[key] 1071 rng = utils.precision(arr.min(), 3) + ", " + utils.precision(arr.max(), 3) 1072 mark_active = "pointdata" 1073 a_scalars = self.dataset.GetPointData().GetScalars() 1074 a_vectors = self.dataset.GetPointData().GetVectors() 1075 a_tensors = self.dataset.GetPointData().GetTensors() 1076 if a_scalars and a_scalars.GetName() == key: 1077 mark_active += " *" 1078 elif a_vectors and a_vectors.GetName() == key: 1079 mark_active += " **" 1080 elif a_tensors and a_tensors.GetName() == key: 1081 mark_active += " ***" 1082 out += mark_active.ljust(14) + f': "{key}" ({arr.dtype}), ndim={arr.ndim}' 1083 out += f", range=({rng})\n" 1084 1085 for key in self.celldata.keys(): 1086 arr = self.celldata[key] 1087 rng = utils.precision(arr.min(), 3) + ", " + utils.precision(arr.max(), 3) 1088 mark_active = "celldata" 1089 a_scalars = self.dataset.GetCellData().GetScalars() 1090 a_vectors = self.dataset.GetCellData().GetVectors() 1091 a_tensors = self.dataset.GetCellData().GetTensors() 1092 if a_scalars and a_scalars.GetName() == key: 1093 mark_active += " *" 1094 elif a_vectors and a_vectors.GetName() == key: 1095 mark_active += " **" 1096 elif a_tensors and a_tensors.GetName() == key: 1097 mark_active += " ***" 1098 out += mark_active.ljust(14) + f': "{key}" ({arr.dtype}), ndim={arr.ndim}' 1099 out += f", range=({rng})\n" 1100 1101 for key in self.metadata.keys(): 1102 arr = self.metadata[key] 1103 out += "metadata".ljust(14) + ": " + f'"{key}" ({len(arr)} values)\n' 1104 1105 return out.rstrip() + "\x1b[0m" 1106 1107 def _repr_html_(self): 1108 """ 1109 HTML representation of the TetMesh object for Jupyter Notebooks. 1110 1111 Returns: 1112 HTML text with the image and some properties. 1113 """ 1114 import io 1115 import base64 1116 from PIL import Image 1117 1118 library_name = "vedo.grids.TetMesh" 1119 help_url = "https://vedo.embl.es/docs/vedo/grids.html#TetMesh" 1120 1121 arr = self.thumbnail() 1122 im = Image.fromarray(arr) 1123 buffered = io.BytesIO() 1124 im.save(buffered, format="PNG", quality=100) 1125 encoded = base64.b64encode(buffered.getvalue()).decode("utf-8") 1126 url = "data:image/png;base64," + encoded 1127 image = f"<img src='{url}'></img>" 1128 1129 bounds = "<br/>".join( 1130 [ 1131 utils.precision(min_x, 4) + " ... " + utils.precision(max_x, 4) 1132 for min_x, max_x in zip(self.bounds()[::2], self.bounds()[1::2]) 1133 ] 1134 ) 1135 1136 help_text = "" 1137 if self.name: 1138 help_text += f"<b> {self.name}:   </b>" 1139 help_text += '<b><a href="' + help_url + '" target="_blank">' + library_name + "</a></b>" 1140 if self.filename: 1141 dots = "" 1142 if len(self.filename) > 30: 1143 dots = "..." 1144 help_text += f"<br/><code><i>({dots}{self.filename[-30:]})</i></code>" 1145 1146 pdata = "" 1147 if self.dataset.GetPointData().GetScalars(): 1148 if self.dataset.GetPointData().GetScalars().GetName(): 1149 name = self.dataset.GetPointData().GetScalars().GetName() 1150 pdata = "<tr><td><b> point data array </b></td><td>" + name + "</td></tr>" 1151 1152 cdata = "" 1153 if self.dataset.GetCellData().GetScalars(): 1154 if self.dataset.GetCellData().GetScalars().GetName(): 1155 name = self.dataset.GetCellData().GetScalars().GetName() 1156 cdata = "<tr><td><b> cell data array </b></td><td>" + name + "</td></tr>" 1157 1158 pts = self.vertices 1159 cm = np.mean(pts, axis=0) 1160 1161 allt = [ 1162 "<table>", 1163 "<tr>", 1164 "<td>", image, "</td>", 1165 "<td style='text-align: center; vertical-align: center;'><br/>", help_text, 1166 "<table>", 1167 "<tr><td><b> bounds </b> <br/> (x/y/z) </td><td>" + str(bounds) + "</td></tr>", 1168 "<tr><td><b> center of mass </b></td><td>" + utils.precision(cm,3) + "</td></tr>", 1169 "<tr><td><b> nr. points / tets </b></td><td>" 1170 + str(self.npoints) + " / " + str(self.ncells) + "</td></tr>", 1171 pdata, 1172 cdata, 1173 "</table>", 1174 "</table>", 1175 ] 1176 return "\n".join(allt) 1177 1178 def compute_quality(self, metric=7) -> np.ndarray: 1179 """ 1180 Calculate functions of quality for the elements of a tetrahedral mesh. 1181 This method adds to the mesh a cell array named "Quality". 1182 1183 Arguments: 1184 metric : (int) 1185 type of estimators: 1186 - EDGE RATIO, 0 1187 - ASPECT RATIO, 1 1188 - RADIUS RATIO, 2 1189 - ASPECT FROBENIUS, 3 1190 - MIN_ANGLE, 4 1191 - COLLAPSE RATIO, 5 1192 - ASPECT GAMMA, 6 1193 - VOLUME, 7 1194 - ... 1195 1196 See class [vtkMeshQuality](https://vtk.org/doc/nightly/html/classvtkMeshQuality.html) 1197 for an explanation of the meaning of each metric.. 1198 """ 1199 qf = vtki.new("MeshQuality") 1200 qf.SetInputData(self.dataset) 1201 qf.SetTetQualityMeasure(metric) 1202 qf.SaveCellQualityOn() 1203 qf.Update() 1204 self._update(qf.GetOutput()) 1205 return utils.vtk2numpy(qf.GetOutput().GetCellData().GetArray("Quality")) 1206 1207 def check_validity(self, tol=0) -> np.ndarray: 1208 """ 1209 Return an array of possible problematic tets following this convention: 1210 ```python 1211 Valid = 0 1212 WrongNumberOfPoints = 01 1213 IntersectingEdges = 02 1214 IntersectingFaces = 04 1215 NoncontiguousEdges = 08 1216 Nonconvex = 10 1217 OrientedIncorrectly = 20 1218 ``` 1219 1220 Arguments: 1221 tol : (float) 1222 This value is used as an epsilon for floating point 1223 equality checks throughout the cell checking process. 1224 """ 1225 vald = vtki.new("CellValidator") 1226 if tol: 1227 vald.SetTolerance(tol) 1228 vald.SetInputData(self.dataset) 1229 vald.Update() 1230 varr = vald.GetOutput().GetCellData().GetArray("ValidityState") 1231 return utils.vtk2numpy(varr) 1232 1233 def decimate(self, scalars_name: str, fraction=0.5, n=0) -> Self: 1234 """ 1235 Downsample the number of tets in a TetMesh to a specified fraction. 1236 Either `fraction` or `n` must be set. 1237 1238 Arguments: 1239 fraction : (float) 1240 the desired final fraction of the total. 1241 n : (int) 1242 the desired number of final tets 1243 1244 .. note:: setting `fraction=0.1` leaves 10% of the original nr of tets. 1245 """ 1246 decimate = vtki.new("UnstructuredGridQuadricDecimation") 1247 decimate.SetInputData(self.dataset) 1248 decimate.SetScalarsName(scalars_name) 1249 1250 if n: # n = desired number of points 1251 decimate.SetNumberOfTetsOutput(n) 1252 else: 1253 decimate.SetTargetReduction(1 - fraction) 1254 decimate.Update() 1255 self._update(decimate.GetOutput()) 1256 self.pipeline = utils.OperationNode( 1257 "decimate", comment=f"array: {scalars_name}", c="#edabab", parents=[self] 1258 ) 1259 return self 1260 1261 def subdvide(self) -> Self: 1262 """ 1263 Increase the number of tetrahedrons of a `TetMesh`. 1264 Subdivides each tetrahedron into twelve smaller tetras. 1265 """ 1266 sd = vtki.new("SubdivideTetra") 1267 sd.SetInputData(self.dataset) 1268 sd.Update() 1269 self._update(sd.GetOutput()) 1270 self.pipeline = utils.OperationNode("subdvide", c="#edabab", parents=[self]) 1271 return self 1272 1273 def generate_random_points(self, n, min_radius=0) -> "vedo.Points": 1274 """ 1275 Generate `n` uniformly distributed random points 1276 inside the tetrahedral mesh. 1277 1278 A new point data array is added to the output points 1279 called "OriginalCellID" which contains the index of 1280 the cell ID in which the point was generated. 1281 1282 Arguments: 1283 n : (int) 1284 number of points to generate. 1285 min_radius: (float) 1286 impose a minimum distance between points. 1287 If `min_radius` is set to 0, the points are 1288 generated uniformly at random inside the mesh. 1289 If `min_radius` is set to a positive value, 1290 the points are generated uniformly at random 1291 inside the mesh, but points closer than `min_radius` 1292 to any other point are discarded. 1293 1294 Returns a `vedo.Points` object. 1295 1296 Note: 1297 Consider using `points.probe(msh)` to interpolate 1298 any existing mesh data onto the points. 1299 1300 Example: 1301 ```python 1302 from vedo import * 1303 tmesh = TetMesh(dataurl + "limb.vtu").alpha(0.2) 1304 pts = tmesh.generate_random_points(20000, min_radius=10) 1305 print(pts.pointdata["OriginalCellID"]) 1306 show(pts, tmesh, axes=1).close() 1307 ``` 1308 """ 1309 cmesh = self.compute_cell_size() 1310 tets = cmesh.cells 1311 verts = cmesh.vertices 1312 cumul = np.cumsum(np.abs(cmesh.celldata["Volume"])) 1313 1314 out_pts = [] 1315 orig_cell = [] 1316 for _ in range(n): 1317 random_area = np.random.random() * cumul[-1] 1318 it = np.searchsorted(cumul, random_area) 1319 A, B, C, D = verts[tets[it]] 1320 r1, r2, r3 = sorted(np.random.random(3)) 1321 p = r1 * A + (r2 - r1) * B + (r3 - r2) * C + (1 - r3) * D 1322 out_pts.append(p) 1323 orig_cell.append(it) 1324 orig_cellnp = np.array(orig_cell, dtype=np.uint32) 1325 1326 vpts = vedo.pointcloud.Points(out_pts) 1327 vpts.pointdata["OriginalCellID"] = orig_cellnp 1328 1329 if min_radius > 0: 1330 vpts.subsample(min_radius, absolute=True) 1331 1332 vpts.point_size(5).color("k1") 1333 vpts.name = "RandomPoints" 1334 vpts.pipeline = utils.OperationNode( 1335 "generate_random_points", c="#edabab", parents=[self]) 1336 return vpts 1337 1338 def isosurface(self, value=True, flying_edges=None) -> "vedo.Mesh": 1339 """ 1340 Return a `vedo.Mesh` isosurface. 1341 The "isosurface" is the surface of the region of points 1342 whose values equal to `value`. 1343 1344 Set `value` to a single value or list of values to compute the isosurface(s). 1345 1346 Note that flying_edges option is not available for `TetMesh`. 1347 """ 1348 if flying_edges is not None: 1349 vedo.logger.warning("flying_edges option is not available for TetMesh.") 1350 1351 if not self.dataset.GetPointData().GetScalars(): 1352 vedo.logger.warning( 1353 "in isosurface() no scalar pointdata found. " 1354 "Mappping cells to points." 1355 ) 1356 self.map_cells_to_points() 1357 scrange = self.dataset.GetPointData().GetScalars().GetRange() 1358 cf = vtki.new("ContourFilter") # vtki.new("ContourGrid") 1359 cf.SetInputData(self.dataset) 1360 1361 if utils.is_sequence(value): 1362 cf.SetNumberOfContours(len(value)) 1363 for i, t in enumerate(value): 1364 cf.SetValue(i, t) 1365 cf.Update() 1366 else: 1367 if value is True: 1368 value = (2 * scrange[0] + scrange[1]) / 3.0 1369 cf.SetValue(0, value) 1370 cf.Update() 1371 1372 msh = Mesh(cf.GetOutput(), c=None) 1373 msh.copy_properties_from(self) 1374 msh.pipeline = utils.OperationNode("isosurface", c="#edabab", parents=[self]) 1375 return msh 1376 1377 def slice(self, origin=(0, 0, 0), normal=(1, 0, 0)) -> "vedo.Mesh": 1378 """ 1379 Return a 2D slice of the mesh by a plane passing through origin and assigned normal. 1380 """ 1381 strn = str(normal) 1382 if strn == "x": normal = (1, 0, 0) 1383 elif strn == "y": normal = (0, 1, 0) 1384 elif strn == "z": normal = (0, 0, 1) 1385 elif strn == "-x": normal = (-1, 0, 0) 1386 elif strn == "-y": normal = (0, -1, 0) 1387 elif strn == "-z": normal = (0, 0, -1) 1388 plane = vtki.new("Plane") 1389 plane.SetOrigin(origin) 1390 plane.SetNormal(normal) 1391 1392 cc = vtki.new("Cutter") 1393 cc.SetInputData(self.dataset) 1394 cc.SetCutFunction(plane) 1395 cc.Update() 1396 msh = Mesh(cc.GetOutput()).flat().lighting("ambient") 1397 msh.copy_properties_from(self) 1398 msh.pipeline = utils.OperationNode("slice", c="#edabab", parents=[self]) 1399 return msh
The class describing tetrahedral meshes.
TetMesh(inputobj=None)
922 def __init__(self, inputobj=None): 923 """ 924 Arguments: 925 inputobj : (vtkUnstructuredGrid, list, str, tetgenpy.TetgenIO) 926 list of points and tet indices, or filename 927 """ 928 super().__init__() 929 930 self.dataset = None 931 932 self.mapper = vtki.new("PolyDataMapper") 933 self._actor = vtki.vtkActor() 934 self._actor.retrieve_object = weak_ref_to(self) 935 self._actor.SetMapper(self.mapper) 936 self.properties = self._actor.GetProperty() 937 938 self.name = "TetMesh" 939 940 # print('TetMesh inputtype', type(inputobj)) 941 942 ################### 943 if inputobj is None: 944 self.dataset = vtki.vtkUnstructuredGrid() 945 946 elif isinstance(inputobj, vtki.vtkUnstructuredGrid): 947 self.dataset = inputobj 948 949 elif isinstance(inputobj, UnstructuredGrid): 950 self.dataset = inputobj.dataset 951 952 elif "TetgenIO" in str(type(inputobj)): # tetgenpy object 953 inputobj = [inputobj.points(), inputobj.tetrahedra()] 954 955 elif isinstance(inputobj, vtki.vtkRectilinearGrid): 956 r2t = vtki.new("RectilinearGridToTetrahedra") 957 r2t.SetInputData(inputobj) 958 r2t.RememberVoxelIdOn() 959 r2t.SetTetraPerCellTo6() 960 r2t.Update() 961 self.dataset = r2t.GetOutput() 962 963 elif isinstance(inputobj, vtki.vtkDataSet): 964 r2t = vtki.new("DataSetTriangleFilter") 965 r2t.SetInputData(inputobj) 966 r2t.TetrahedraOnlyOn() 967 r2t.Update() 968 self.dataset = r2t.GetOutput() 969 970 elif isinstance(inputobj, str): 971 if "https://" in inputobj: 972 inputobj = download(inputobj, verbose=False) 973 if inputobj.endswith(".vtu"): 974 reader = vtki.new("XMLUnstructuredGridReader") 975 else: 976 reader = vtki.new("UnstructuredGridReader") 977 978 if not os.path.isfile(inputobj): 979 # for some reason vtk Reader does not complain 980 vedo.logger.error(f"file {inputobj} not found") 981 raise FileNotFoundError 982 983 self.filename = inputobj 984 reader.SetFileName(inputobj) 985 reader.Update() 986 ug = reader.GetOutput() 987 988 tt = vtki.new("DataSetTriangleFilter") 989 tt.SetInputData(ug) 990 tt.SetTetrahedraOnly(True) 991 tt.Update() 992 self.dataset = tt.GetOutput() 993 994 ############################### 995 if utils.is_sequence(inputobj): 996 self.dataset = vtki.vtkUnstructuredGrid() 997 998 points, cells = inputobj 999 if len(points) == 0: 1000 return 1001 if not utils.is_sequence(points[0]): 1002 return 1003 if len(cells) == 0: 1004 return 1005 1006 if not utils.is_sequence(cells[0]): 1007 tets = [] 1008 nf = cells[0] + 1 1009 for i, cl in enumerate(cells): 1010 if i in (nf, 0): 1011 k = i + 1 1012 nf = cl + k 1013 cell = [cells[j + k] for j in range(cl)] 1014 tets.append(cell) 1015 cells = tets 1016 1017 source_points = vtki.vtkPoints() 1018 varr = utils.numpy2vtk(points, dtype=np.float32) 1019 source_points.SetData(varr) 1020 self.dataset.SetPoints(source_points) 1021 1022 source_tets = vtki.vtkCellArray() 1023 for f in cells: 1024 ele = vtki.vtkTetra() 1025 pid = ele.GetPointIds() 1026 for i, fi in enumerate(f): 1027 pid.SetId(i, fi) 1028 source_tets.InsertNextCell(ele) 1029 self.dataset.SetCells(vtki.cell_types["TETRA"], source_tets) 1030 1031 if not self.dataset: 1032 vedo.logger.error(f"cannot understand input type {type(inputobj)}") 1033 return 1034 1035 self.properties.SetColor(0.352, 0.612, 0.996) # blue7 1036 1037 self.pipeline = utils.OperationNode( 1038 self, comment=f"#tets {self.dataset.GetNumberOfCells()}", c="#9e2a2b" 1039 )
Arguments:
- inputobj : (vtkUnstructuredGrid, list, str, tetgenpy.TetgenIO) list of points and tet indices, or filename
def
compute_quality(self, metric=7) -> numpy.ndarray:
1178 def compute_quality(self, metric=7) -> np.ndarray: 1179 """ 1180 Calculate functions of quality for the elements of a tetrahedral mesh. 1181 This method adds to the mesh a cell array named "Quality". 1182 1183 Arguments: 1184 metric : (int) 1185 type of estimators: 1186 - EDGE RATIO, 0 1187 - ASPECT RATIO, 1 1188 - RADIUS RATIO, 2 1189 - ASPECT FROBENIUS, 3 1190 - MIN_ANGLE, 4 1191 - COLLAPSE RATIO, 5 1192 - ASPECT GAMMA, 6 1193 - VOLUME, 7 1194 - ... 1195 1196 See class [vtkMeshQuality](https://vtk.org/doc/nightly/html/classvtkMeshQuality.html) 1197 for an explanation of the meaning of each metric.. 1198 """ 1199 qf = vtki.new("MeshQuality") 1200 qf.SetInputData(self.dataset) 1201 qf.SetTetQualityMeasure(metric) 1202 qf.SaveCellQualityOn() 1203 qf.Update() 1204 self._update(qf.GetOutput()) 1205 return utils.vtk2numpy(qf.GetOutput().GetCellData().GetArray("Quality"))
Calculate functions of quality for the elements of a tetrahedral 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..
def
check_validity(self, tol=0) -> numpy.ndarray:
1207 def check_validity(self, tol=0) -> np.ndarray: 1208 """ 1209 Return an array of possible problematic tets following this convention: 1210 ```python 1211 Valid = 0 1212 WrongNumberOfPoints = 01 1213 IntersectingEdges = 02 1214 IntersectingFaces = 04 1215 NoncontiguousEdges = 08 1216 Nonconvex = 10 1217 OrientedIncorrectly = 20 1218 ``` 1219 1220 Arguments: 1221 tol : (float) 1222 This value is used as an epsilon for floating point 1223 equality checks throughout the cell checking process. 1224 """ 1225 vald = vtki.new("CellValidator") 1226 if tol: 1227 vald.SetTolerance(tol) 1228 vald.SetInputData(self.dataset) 1229 vald.Update() 1230 varr = vald.GetOutput().GetCellData().GetArray("ValidityState") 1231 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.
def
decimate(self, scalars_name: str, fraction=0.5, n=0) -> Self:
1233 def decimate(self, scalars_name: str, fraction=0.5, n=0) -> Self: 1234 """ 1235 Downsample the number of tets in a TetMesh to a specified fraction. 1236 Either `fraction` or `n` must be set. 1237 1238 Arguments: 1239 fraction : (float) 1240 the desired final fraction of the total. 1241 n : (int) 1242 the desired number of final tets 1243 1244 .. note:: setting `fraction=0.1` leaves 10% of the original nr of tets. 1245 """ 1246 decimate = vtki.new("UnstructuredGridQuadricDecimation") 1247 decimate.SetInputData(self.dataset) 1248 decimate.SetScalarsName(scalars_name) 1249 1250 if n: # n = desired number of points 1251 decimate.SetNumberOfTetsOutput(n) 1252 else: 1253 decimate.SetTargetReduction(1 - fraction) 1254 decimate.Update() 1255 self._update(decimate.GetOutput()) 1256 self.pipeline = utils.OperationNode( 1257 "decimate", comment=f"array: {scalars_name}", c="#edabab", parents=[self] 1258 ) 1259 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.
def
subdvide(self) -> Self:
1261 def subdvide(self) -> Self: 1262 """ 1263 Increase the number of tetrahedrons of a `TetMesh`. 1264 Subdivides each tetrahedron into twelve smaller tetras. 1265 """ 1266 sd = vtki.new("SubdivideTetra") 1267 sd.SetInputData(self.dataset) 1268 sd.Update() 1269 self._update(sd.GetOutput()) 1270 self.pipeline = utils.OperationNode("subdvide", c="#edabab", parents=[self]) 1271 return self
Increase the number of tetrahedrons of a TetMesh.
Subdivides each tetrahedron into twelve smaller tetras.
1273 def generate_random_points(self, n, min_radius=0) -> "vedo.Points": 1274 """ 1275 Generate `n` uniformly distributed random points 1276 inside the tetrahedral mesh. 1277 1278 A new point data array is added to the output points 1279 called "OriginalCellID" which contains the index of 1280 the cell ID in which the point was generated. 1281 1282 Arguments: 1283 n : (int) 1284 number of points to generate. 1285 min_radius: (float) 1286 impose a minimum distance between points. 1287 If `min_radius` is set to 0, the points are 1288 generated uniformly at random inside the mesh. 1289 If `min_radius` is set to a positive value, 1290 the points are generated uniformly at random 1291 inside the mesh, but points closer than `min_radius` 1292 to any other point are discarded. 1293 1294 Returns a `vedo.Points` object. 1295 1296 Note: 1297 Consider using `points.probe(msh)` to interpolate 1298 any existing mesh data onto the points. 1299 1300 Example: 1301 ```python 1302 from vedo import * 1303 tmesh = TetMesh(dataurl + "limb.vtu").alpha(0.2) 1304 pts = tmesh.generate_random_points(20000, min_radius=10) 1305 print(pts.pointdata["OriginalCellID"]) 1306 show(pts, tmesh, axes=1).close() 1307 ``` 1308 """ 1309 cmesh = self.compute_cell_size() 1310 tets = cmesh.cells 1311 verts = cmesh.vertices 1312 cumul = np.cumsum(np.abs(cmesh.celldata["Volume"])) 1313 1314 out_pts = [] 1315 orig_cell = [] 1316 for _ in range(n): 1317 random_area = np.random.random() * cumul[-1] 1318 it = np.searchsorted(cumul, random_area) 1319 A, B, C, D = verts[tets[it]] 1320 r1, r2, r3 = sorted(np.random.random(3)) 1321 p = r1 * A + (r2 - r1) * B + (r3 - r2) * C + (1 - r3) * D 1322 out_pts.append(p) 1323 orig_cell.append(it) 1324 orig_cellnp = np.array(orig_cell, dtype=np.uint32) 1325 1326 vpts = vedo.pointcloud.Points(out_pts) 1327 vpts.pointdata["OriginalCellID"] = orig_cellnp 1328 1329 if min_radius > 0: 1330 vpts.subsample(min_radius, absolute=True) 1331 1332 vpts.point_size(5).color("k1") 1333 vpts.name = "RandomPoints" 1334 vpts.pipeline = utils.OperationNode( 1335 "generate_random_points", c="#edabab", parents=[self]) 1336 return vpts
Generate n uniformly distributed random points
inside the tetrahedral mesh.
A new point data array is added to the output points called "OriginalCellID" which contains the index of the cell ID in which the point was generated.
Arguments:
- n : (int) number of points to generate.
- min_radius: (float)
impose a minimum distance between points.
If
min_radiusis set to 0, the points are generated uniformly at random inside the mesh. Ifmin_radiusis set to a positive value, the points are generated uniformly at random inside the mesh, but points closer thanmin_radiusto any other point are discarded.
Returns a vedo.Points object.
Note:
Consider using
points.probe(msh)to interpolate any existing mesh data onto the points.
Example:
from vedo import *
tmesh = TetMesh(dataurl + "limb.vtu").alpha(0.2)
pts = tmesh.generate_random_points(20000, min_radius=10)
print(pts.pointdata["OriginalCellID"])
show(pts, tmesh, axes=1).close()
1338 def isosurface(self, value=True, flying_edges=None) -> "vedo.Mesh": 1339 """ 1340 Return a `vedo.Mesh` isosurface. 1341 The "isosurface" is the surface of the region of points 1342 whose values equal to `value`. 1343 1344 Set `value` to a single value or list of values to compute the isosurface(s). 1345 1346 Note that flying_edges option is not available for `TetMesh`. 1347 """ 1348 if flying_edges is not None: 1349 vedo.logger.warning("flying_edges option is not available for TetMesh.") 1350 1351 if not self.dataset.GetPointData().GetScalars(): 1352 vedo.logger.warning( 1353 "in isosurface() no scalar pointdata found. " 1354 "Mappping cells to points." 1355 ) 1356 self.map_cells_to_points() 1357 scrange = self.dataset.GetPointData().GetScalars().GetRange() 1358 cf = vtki.new("ContourFilter") # vtki.new("ContourGrid") 1359 cf.SetInputData(self.dataset) 1360 1361 if utils.is_sequence(value): 1362 cf.SetNumberOfContours(len(value)) 1363 for i, t in enumerate(value): 1364 cf.SetValue(i, t) 1365 cf.Update() 1366 else: 1367 if value is True: 1368 value = (2 * scrange[0] + scrange[1]) / 3.0 1369 cf.SetValue(0, value) 1370 cf.Update() 1371 1372 msh = Mesh(cf.GetOutput(), c=None) 1373 msh.copy_properties_from(self) 1374 msh.pipeline = utils.OperationNode("isosurface", c="#edabab", parents=[self]) 1375 return msh
Return a vedo.Mesh isosurface.
The "isosurface" is the surface of the region of points
whose values equal to value.
Set value to a single value or list of values to compute the isosurface(s).
Note that flying_edges option is not available for TetMesh.
1377 def slice(self, origin=(0, 0, 0), normal=(1, 0, 0)) -> "vedo.Mesh": 1378 """ 1379 Return a 2D slice of the mesh by a plane passing through origin and assigned normal. 1380 """ 1381 strn = str(normal) 1382 if strn == "x": normal = (1, 0, 0) 1383 elif strn == "y": normal = (0, 1, 0) 1384 elif strn == "z": normal = (0, 0, 1) 1385 elif strn == "-x": normal = (-1, 0, 0) 1386 elif strn == "-y": normal = (0, -1, 0) 1387 elif strn == "-z": normal = (0, 0, -1) 1388 plane = vtki.new("Plane") 1389 plane.SetOrigin(origin) 1390 plane.SetNormal(normal) 1391 1392 cc = vtki.new("Cutter") 1393 cc.SetInputData(self.dataset) 1394 cc.SetCutFunction(plane) 1395 cc.Update() 1396 msh = Mesh(cc.GetOutput()).flat().lighting("ambient") 1397 msh.copy_properties_from(self) 1398 msh.pipeline = utils.OperationNode("slice", c="#edabab", parents=[self]) 1399 return msh
Return a 2D slice of the mesh by a plane passing through origin and assigned normal.
Inherited Members
- UnstructuredGrid
- actor
- merge
- copy
- clone
- bounds
- threshold
- shrink
- tomesh
- cell_types_array
- extract_cells_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
- update_dataset
- 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
- keep_cell_types
- 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
- generate_random_data
- integrate_data
- write
- signed_distance
- unsigned_distance
- smooth_data
- compute_streamlines
- 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
1403class RectilinearGrid(PointAlgorithms, MeshVisual): 1404 """ 1405 Build a rectilinear grid. 1406 """ 1407 1408 def __init__(self, inputobj=None): 1409 """ 1410 A RectilinearGrid is a dataset where edges are parallel to the coordinate axes. 1411 It can be thought of as a tessellation of a box in 3D space, similar to a `Volume` 1412 except that the cells are not necessarily cubes, but they can have different lengths 1413 along each axis. 1414 This can be useful to describe a volume with variable resolution where one needs 1415 to represent a region with higher detail with respect to another region. 1416 1417 Arguments: 1418 inputobj : (vtkRectilinearGrid, list, str) 1419 list of points and tet indices, or filename 1420 1421 Example: 1422 ```python 1423 from vedo import RectilinearGrid, show 1424 1425 xcoords = 7 + np.sqrt(np.arange(0,2500,25)) 1426 ycoords = np.arange(0, 20) 1427 zcoords = np.arange(0, 20) 1428 1429 rgrid = RectilinearGrid([xcoords, ycoords, zcoords]) 1430 1431 print(rgrid) 1432 print(rgrid.x_coordinates().shape) 1433 print(rgrid.compute_structured_coords([20,10,11])) 1434 1435 msh = rgrid.tomesh().lw(1) 1436 1437 show(msh, axes=1, viewup="z") 1438 ``` 1439 """ 1440 1441 super().__init__() 1442 1443 self.dataset = None 1444 1445 self.mapper = vtki.new("PolyDataMapper") 1446 self._actor = vtki.vtkActor() 1447 self._actor.retrieve_object = weak_ref_to(self) 1448 self._actor.SetMapper(self.mapper) 1449 self.properties = self._actor.GetProperty() 1450 1451 self.transform = LinearTransform() 1452 self.point_locator = None 1453 self.cell_locator = None 1454 self.line_locator = None 1455 1456 self.name = "RectilinearGrid" 1457 self.filename = "" 1458 1459 self.info = {} 1460 self.time = time.time() 1461 1462 ############################### 1463 if inputobj is None: 1464 self.dataset = vtki.vtkRectilinearGrid() 1465 1466 elif isinstance(inputobj, vtki.vtkRectilinearGrid): 1467 self.dataset = inputobj 1468 1469 elif isinstance(inputobj, RectilinearGrid): 1470 self.dataset = inputobj.dataset 1471 1472 elif isinstance(inputobj, str): 1473 if "https://" in inputobj: 1474 inputobj = download(inputobj, verbose=False) 1475 if inputobj.endswith(".vtr"): 1476 reader = vtki.new("XMLRectilinearGridReader") 1477 else: 1478 reader = vtki.new("RectilinearGridReader") 1479 self.filename = inputobj 1480 reader.SetFileName(inputobj) 1481 reader.Update() 1482 self.dataset = reader.GetOutput() 1483 1484 elif utils.is_sequence(inputobj): 1485 self.dataset = vtki.vtkRectilinearGrid() 1486 xcoords, ycoords, zcoords = inputobj 1487 nx, ny, nz = len(xcoords), len(ycoords), len(zcoords) 1488 self.dataset.SetDimensions(nx, ny, nz) 1489 self.dataset.SetXCoordinates(utils.numpy2vtk(xcoords)) 1490 self.dataset.SetYCoordinates(utils.numpy2vtk(ycoords)) 1491 self.dataset.SetZCoordinates(utils.numpy2vtk(zcoords)) 1492 1493 ############################### 1494 1495 if not self.dataset: 1496 vedo.logger.error(f"RectilinearGrid: cannot understand input type {type(inputobj)}") 1497 return 1498 1499 self.properties.SetColor(0.352, 0.612, 0.996) # blue7 1500 1501 self.pipeline = utils.OperationNode( 1502 self, comment=f"#tets {self.dataset.GetNumberOfCells()}", c="#9e2a2b" 1503 ) 1504 1505 @property 1506 def actor(self): 1507 """Return the `vtkActor` of the object.""" 1508 gf = vtki.new("GeometryFilter") 1509 gf.SetInputData(self.dataset) 1510 gf.Update() 1511 self.mapper.SetInputData(gf.GetOutput()) 1512 self.mapper.Modified() 1513 return self._actor 1514 1515 @actor.setter 1516 def actor(self, _): 1517 pass 1518 1519 def _update(self, data, reset_locators=False): 1520 self.dataset = data 1521 if reset_locators: 1522 self.cell_locator = None 1523 self.point_locator = None 1524 return self 1525 1526 ################################################################## 1527 def __str__(self): 1528 """Print a summary for the `RectilinearGrid` object.""" 1529 module = self.__class__.__module__ 1530 name = self.__class__.__name__ 1531 out = vedo.printc( 1532 f"{module}.{name} at ({hex(self.memory_address())})".ljust(75), 1533 c="c", bold=True, invert=True, return_string=True, 1534 ) 1535 out += "\x1b[0m\x1b[36;1m" 1536 1537 out += "name".ljust(14) + ": " + str(self.name) + "\n" 1538 if self.filename: 1539 out += "filename".ljust(14) + ": " + str(self.filename) + "\n" 1540 1541 out += "dimensions".ljust(14) + ": " + str(self.dataset.GetDimensions()) + "\n" 1542 1543 out += "center".ljust(14) + ": " 1544 out += utils.precision(self.dataset.GetCenter(), 6) + "\n" 1545 1546 bnds = self.bounds() 1547 bx1, bx2 = utils.precision(bnds[0], 3), utils.precision(bnds[1], 3) 1548 by1, by2 = utils.precision(bnds[2], 3), utils.precision(bnds[3], 3) 1549 bz1, bz2 = utils.precision(bnds[4], 3), utils.precision(bnds[5], 3) 1550 out += "bounds".ljust(14) + ":" 1551 out += " x=(" + bx1 + ", " + bx2 + ")," 1552 out += " y=(" + by1 + ", " + by2 + ")," 1553 out += " z=(" + bz1 + ", " + bz2 + ")\n" 1554 1555 out += "memory size".ljust(14) + ": " 1556 out += utils.precision(self.dataset.GetActualMemorySize() / 1024, 2) + " MB\n" 1557 1558 for key in self.pointdata.keys(): 1559 arr = self.pointdata[key] 1560 rng = utils.precision(arr.min(), 3) + ", " + utils.precision(arr.max(), 3) 1561 mark_active = "pointdata" 1562 a_scalars = self.dataset.GetPointData().GetScalars() 1563 a_vectors = self.dataset.GetPointData().GetVectors() 1564 a_tensors = self.dataset.GetPointData().GetTensors() 1565 if a_scalars and a_scalars.GetName() == key: 1566 mark_active += " *" 1567 elif a_vectors and a_vectors.GetName() == key: 1568 mark_active += " **" 1569 elif a_tensors and a_tensors.GetName() == key: 1570 mark_active += " ***" 1571 out += mark_active.ljust(14) + f': "{key}" ({arr.dtype}), ndim={arr.ndim}' 1572 out += f", range=({rng})\n" 1573 1574 for key in self.celldata.keys(): 1575 arr = self.celldata[key] 1576 rng = utils.precision(arr.min(), 3) + ", " + utils.precision(arr.max(), 3) 1577 mark_active = "celldata" 1578 a_scalars = self.dataset.GetCellData().GetScalars() 1579 a_vectors = self.dataset.GetCellData().GetVectors() 1580 a_tensors = self.dataset.GetCellData().GetTensors() 1581 if a_scalars and a_scalars.GetName() == key: 1582 mark_active += " *" 1583 elif a_vectors and a_vectors.GetName() == key: 1584 mark_active += " **" 1585 elif a_tensors and a_tensors.GetName() == key: 1586 mark_active += " ***" 1587 out += mark_active.ljust(14) + f': "{key}" ({arr.dtype}), ndim={arr.ndim}' 1588 out += f", range=({rng})\n" 1589 1590 for key in self.metadata.keys(): 1591 arr = self.metadata[key] 1592 out += "metadata".ljust(14) + ": " + f'"{key}" ({len(arr)} values)\n' 1593 1594 return out.rstrip() + "\x1b[0m" 1595 1596 def _repr_html_(self): 1597 """ 1598 HTML representation of the RectilinearGrid object for Jupyter Notebooks. 1599 1600 Returns: 1601 HTML text with the image and some properties. 1602 """ 1603 import io 1604 import base64 1605 from PIL import Image 1606 1607 library_name = "vedo.grids.RectilinearGrid" 1608 help_url = "https://vedo.embl.es/docs/vedo/grids.html#RectilinearGrid" 1609 1610 m = self.tomesh().linewidth(1).lighting("off") 1611 arr= m.thumbnail(zoom=1, elevation=-30, azimuth=-30) 1612 1613 im = Image.fromarray(arr) 1614 buffered = io.BytesIO() 1615 im.save(buffered, format="PNG", quality=100) 1616 encoded = base64.b64encode(buffered.getvalue()).decode("utf-8") 1617 url = "data:image/png;base64," + encoded 1618 image = f"<img src='{url}'></img>" 1619 1620 bounds = "<br/>".join( 1621 [ 1622 utils.precision(min_x, 4) + " ... " + utils.precision(max_x, 4) 1623 for min_x, max_x in zip(self.bounds()[::2], self.bounds()[1::2]) 1624 ] 1625 ) 1626 1627 help_text = "" 1628 if self.name: 1629 help_text += f"<b> {self.name}:   </b>" 1630 help_text += '<b><a href="' + help_url + '" target="_blank">' + library_name + "</a></b>" 1631 if self.filename: 1632 dots = "" 1633 if len(self.filename) > 30: 1634 dots = "..." 1635 help_text += f"<br/><code><i>({dots}{self.filename[-30:]})</i></code>" 1636 1637 pdata = "" 1638 if self.dataset.GetPointData().GetScalars(): 1639 if self.dataset.GetPointData().GetScalars().GetName(): 1640 name = self.dataset.GetPointData().GetScalars().GetName() 1641 pdata = "<tr><td><b> point data array </b></td><td>" + name + "</td></tr>" 1642 1643 cdata = "" 1644 if self.dataset.GetCellData().GetScalars(): 1645 if self.dataset.GetCellData().GetScalars().GetName(): 1646 name = self.dataset.GetCellData().GetScalars().GetName() 1647 cdata = "<tr><td><b> cell data array </b></td><td>" + name + "</td></tr>" 1648 1649 pts = self.vertices 1650 cm = np.mean(pts, axis=0) 1651 1652 all = [ 1653 "<table>", 1654 "<tr>", 1655 "<td>", image, "</td>", 1656 "<td style='text-align: center; vertical-align: center;'><br/>", help_text, 1657 "<table>", 1658 "<tr><td><b> bounds </b> <br/> (x/y/z) </td><td>" + str(bounds) + "</td></tr>", 1659 "<tr><td><b> center of mass </b></td><td>" + utils.precision(cm,3) + "</td></tr>", 1660 "<tr><td><b> nr. points / cells </b></td><td>" 1661 + str(self.npoints) + " / " + str(self.ncells) + "</td></tr>", 1662 pdata, 1663 cdata, 1664 "</table>", 1665 "</table>", 1666 ] 1667 return "\n".join(all) 1668 1669 def dimensions(self) -> np.ndarray: 1670 """Return the number of points in the x, y and z directions.""" 1671 return np.array(self.dataset.GetDimensions()) 1672 1673 def x_coordinates(self) -> np.ndarray: 1674 """Return the x-coordinates of the grid.""" 1675 return utils.vtk2numpy(self.dataset.GetXCoordinates()) 1676 1677 def y_coordinates(self) -> np.ndarray: 1678 """Return the y-coordinates of the grid.""" 1679 return utils.vtk2numpy(self.dataset.GetYCoordinates()) 1680 1681 def z_coordinates(self) -> np.ndarray: 1682 """Return the z-coordinates of the grid.""" 1683 return utils.vtk2numpy(self.dataset.GetZCoordinates()) 1684 1685 def is_point_visible(self, pid: int) -> bool: 1686 """Return True if point `pid` is visible.""" 1687 return self.dataset.IsPointVisible(pid) 1688 1689 def is_cell_visible(self, cid: int) -> bool: 1690 """Return True if cell `cid` is visible.""" 1691 return self.dataset.IsCellVisible(cid) 1692 1693 def has_blank_points(self) -> bool: 1694 """Return True if the grid has blank points.""" 1695 return self.dataset.HasAnyBlankPoints() 1696 1697 def has_blank_cells(self) -> bool: 1698 """Return True if the grid has blank cells.""" 1699 return self.dataset.HasAnyBlankCells() 1700 1701 def compute_structured_coords(self, x: list) -> dict: 1702 """ 1703 Convenience function computes the structured coordinates for a point `x`. 1704 1705 This method returns a dictionary with keys `ijk`, `pcoords` and `inside`. 1706 The cell is specified by the array `ijk`. 1707 and the parametric coordinates in the cell are specified with `pcoords`. 1708 Value of `inside` is False if the point x is outside of the grid. 1709 """ 1710 ijk = [0, 0, 0] 1711 pcoords = [0., 0., 0.] 1712 inout = self.dataset.ComputeStructuredCoordinates(x, ijk, pcoords) 1713 return {"ijk": np.array(ijk), "pcoords": np.array(pcoords), "inside": bool(inout)} 1714 1715 def compute_pointid(self, ijk: int) -> int: 1716 """Given a location in structured coordinates (i-j-k), return the point id.""" 1717 return self.dataset.ComputePointId(ijk) 1718 1719 def compute_cellid(self, ijk: int) -> int: 1720 """Given a location in structured coordinates (i-j-k), return the cell id.""" 1721 return self.dataset.ComputeCellId(ijk) 1722 1723 def find_point(self, x: list) -> int: 1724 """Given a position `x`, return the id of the closest point.""" 1725 return self.dataset.FindPoint(x) 1726 1727 def find_cell(self, x: list) -> dict: 1728 """Given a position `x`, return the id of the closest cell.""" 1729 cell = vtki.vtkHexagonalPrism() 1730 cellid = vtki.mutable(0) 1731 tol2 = 0.001 # vtki.mutable(0) 1732 subid = vtki.mutable(0) 1733 pcoords = [0.0, 0.0, 0.0] 1734 weights = [0.0, 0.0, 0.0] 1735 res = self.dataset.FindCell(x, cell, cellid, tol2, subid, pcoords, weights) 1736 result = {} 1737 result["cellid"] = cellid 1738 result["subid"] = subid 1739 result["pcoords"] = pcoords 1740 result["weights"] = weights 1741 result["status"] = res 1742 return result 1743 1744 def clone(self, deep=True) -> "RectilinearGrid": 1745 """Return a clone copy of the RectilinearGrid. Alias of `copy()`.""" 1746 if deep: 1747 newrg = vtki.vtkRectilinearGrid() 1748 newrg.CopyStructure(self.dataset) 1749 newrg.CopyAttributes(self.dataset) 1750 newvol = RectilinearGrid(newrg) 1751 else: 1752 newvol = RectilinearGrid(self.dataset) 1753 1754 prop = vtki.vtkProperty() 1755 prop.DeepCopy(self.properties) 1756 newvol.actor.SetProperty(prop) 1757 newvol.properties = prop 1758 newvol.pipeline = utils.OperationNode("clone", parents=[self], c="#bbd0ff", shape="diamond") 1759 return newvol 1760 1761 def bounds(self) -> np.ndarray: 1762 """ 1763 Get the object bounds. 1764 Returns a list in format `[xmin,xmax, ymin,ymax, zmin,zmax]`. 1765 """ 1766 # OVERRIDE CommonAlgorithms.bounds() which is too slow 1767 return np.array(self.dataset.GetBounds()) 1768 1769 def isosurface(self, value=None) -> "vedo.Mesh": 1770 """ 1771 Return a `Mesh` isosurface extracted from the object. 1772 1773 Set `value` as single float or list of values to draw the isosurface(s). 1774 """ 1775 scrange = self.dataset.GetScalarRange() 1776 1777 cf = vtki.new("ContourFilter") 1778 cf.UseScalarTreeOn() 1779 cf.SetInputData(self.dataset) 1780 cf.ComputeNormalsOn() 1781 1782 if value is None: 1783 value = (2 * scrange[0] + scrange[1]) / 3.0 1784 # print("automatic isosurface value =", value) 1785 cf.SetValue(0, value) 1786 else: 1787 if utils.is_sequence(value): 1788 cf.SetNumberOfContours(len(value)) 1789 for i, t in enumerate(value): 1790 cf.SetValue(i, t) 1791 else: 1792 cf.SetValue(0, value) 1793 1794 cf.Update() 1795 poly = cf.GetOutput() 1796 1797 out = vedo.mesh.Mesh(poly, c=None).phong() 1798 out.mapper.SetScalarRange(scrange[0], scrange[1]) 1799 1800 out.pipeline = utils.OperationNode( 1801 "isosurface", 1802 parents=[self], 1803 comment=f"#pts {out.dataset.GetNumberOfPoints()}", 1804 c="#4cc9f0:#e9c46a", 1805 ) 1806 return out 1807 1808 def cut_with_plane(self, origin=(0, 0, 0), normal="x") -> "UnstructuredGrid": 1809 """ 1810 Cut the object with the plane defined by a point and a normal. 1811 1812 Arguments: 1813 origin : (list) 1814 the cutting plane goes through this point 1815 normal : (list, str) 1816 normal vector to the cutting plane 1817 """ 1818 strn = str(normal) 1819 if strn == "x": normal = (1, 0, 0) 1820 elif strn == "y": normal = (0, 1, 0) 1821 elif strn == "z": normal = (0, 0, 1) 1822 elif strn == "-x": normal = (-1, 0, 0) 1823 elif strn == "-y": normal = (0, -1, 0) 1824 elif strn == "-z": normal = (0, 0, -1) 1825 plane = vtki.new("Plane") 1826 plane.SetOrigin(origin) 1827 plane.SetNormal(normal) 1828 clipper = vtki.new("ClipDataSet") 1829 clipper.SetInputData(self.dataset) 1830 clipper.SetClipFunction(plane) 1831 clipper.GenerateClipScalarsOff() 1832 clipper.GenerateClippedOutputOff() 1833 clipper.SetValue(0) 1834 clipper.Update() 1835 cout = clipper.GetOutput() 1836 ug = vedo.UnstructuredGrid(cout) 1837 if isinstance(self, UnstructuredGrid): 1838 self._update(cout) 1839 self.pipeline = utils.OperationNode("cut_with_plane", parents=[self], c="#9e2a2b") 1840 return self 1841 ug.pipeline = utils.OperationNode("cut_with_plane", parents=[self], c="#9e2a2b") 1842 return ug 1843 1844 def cut_with_mesh(self, mesh, invert=False, whole_cells=False, on_boundary=False) -> "UnstructuredGrid": 1845 """ 1846 Cut a `RectilinearGrid` with a `Mesh`. 1847 1848 Use `invert` to return cut off part of the input object. 1849 """ 1850 ug = self.dataset 1851 1852 ippd = vtki.new("ImplicitPolyDataDistance") 1853 ippd.SetInput(mesh.dataset) 1854 1855 if whole_cells or on_boundary: 1856 clipper = vtki.new("ExtractGeometry") 1857 clipper.SetInputData(ug) 1858 clipper.SetImplicitFunction(ippd) 1859 clipper.SetExtractInside(not invert) 1860 clipper.SetExtractBoundaryCells(False) 1861 if on_boundary: 1862 clipper.SetExtractBoundaryCells(True) 1863 clipper.SetExtractOnlyBoundaryCells(True) 1864 else: 1865 signed_dists = vtki.vtkFloatArray() 1866 signed_dists.SetNumberOfComponents(1) 1867 signed_dists.SetName("SignedDistance") 1868 for pointId in range(ug.GetNumberOfPoints()): 1869 p = ug.GetPoint(pointId) 1870 signed_dist = ippd.EvaluateFunction(p) 1871 signed_dists.InsertNextValue(signed_dist) 1872 ug.GetPointData().AddArray(signed_dists) 1873 ug.GetPointData().SetActiveScalars("SignedDistance") # NEEDED 1874 clipper = vtki.new("ClipDataSet") 1875 clipper.SetInputData(ug) 1876 clipper.SetInsideOut(not invert) 1877 clipper.SetValue(0.0) 1878 1879 clipper.Update() 1880 1881 out = UnstructuredGrid(clipper.GetOutput()) 1882 out.pipeline = utils.OperationNode("cut_with_mesh", parents=[self], c="#9e2a2b") 1883 return out
Build a rectilinear grid.
RectilinearGrid(inputobj=None)
1408 def __init__(self, inputobj=None): 1409 """ 1410 A RectilinearGrid is a dataset where edges are parallel to the coordinate axes. 1411 It can be thought of as a tessellation of a box in 3D space, similar to a `Volume` 1412 except that the cells are not necessarily cubes, but they can have different lengths 1413 along each axis. 1414 This can be useful to describe a volume with variable resolution where one needs 1415 to represent a region with higher detail with respect to another region. 1416 1417 Arguments: 1418 inputobj : (vtkRectilinearGrid, list, str) 1419 list of points and tet indices, or filename 1420 1421 Example: 1422 ```python 1423 from vedo import RectilinearGrid, show 1424 1425 xcoords = 7 + np.sqrt(np.arange(0,2500,25)) 1426 ycoords = np.arange(0, 20) 1427 zcoords = np.arange(0, 20) 1428 1429 rgrid = RectilinearGrid([xcoords, ycoords, zcoords]) 1430 1431 print(rgrid) 1432 print(rgrid.x_coordinates().shape) 1433 print(rgrid.compute_structured_coords([20,10,11])) 1434 1435 msh = rgrid.tomesh().lw(1) 1436 1437 show(msh, axes=1, viewup="z") 1438 ``` 1439 """ 1440 1441 super().__init__() 1442 1443 self.dataset = None 1444 1445 self.mapper = vtki.new("PolyDataMapper") 1446 self._actor = vtki.vtkActor() 1447 self._actor.retrieve_object = weak_ref_to(self) 1448 self._actor.SetMapper(self.mapper) 1449 self.properties = self._actor.GetProperty() 1450 1451 self.transform = LinearTransform() 1452 self.point_locator = None 1453 self.cell_locator = None 1454 self.line_locator = None 1455 1456 self.name = "RectilinearGrid" 1457 self.filename = "" 1458 1459 self.info = {} 1460 self.time = time.time() 1461 1462 ############################### 1463 if inputobj is None: 1464 self.dataset = vtki.vtkRectilinearGrid() 1465 1466 elif isinstance(inputobj, vtki.vtkRectilinearGrid): 1467 self.dataset = inputobj 1468 1469 elif isinstance(inputobj, RectilinearGrid): 1470 self.dataset = inputobj.dataset 1471 1472 elif isinstance(inputobj, str): 1473 if "https://" in inputobj: 1474 inputobj = download(inputobj, verbose=False) 1475 if inputobj.endswith(".vtr"): 1476 reader = vtki.new("XMLRectilinearGridReader") 1477 else: 1478 reader = vtki.new("RectilinearGridReader") 1479 self.filename = inputobj 1480 reader.SetFileName(inputobj) 1481 reader.Update() 1482 self.dataset = reader.GetOutput() 1483 1484 elif utils.is_sequence(inputobj): 1485 self.dataset = vtki.vtkRectilinearGrid() 1486 xcoords, ycoords, zcoords = inputobj 1487 nx, ny, nz = len(xcoords), len(ycoords), len(zcoords) 1488 self.dataset.SetDimensions(nx, ny, nz) 1489 self.dataset.SetXCoordinates(utils.numpy2vtk(xcoords)) 1490 self.dataset.SetYCoordinates(utils.numpy2vtk(ycoords)) 1491 self.dataset.SetZCoordinates(utils.numpy2vtk(zcoords)) 1492 1493 ############################### 1494 1495 if not self.dataset: 1496 vedo.logger.error(f"RectilinearGrid: cannot understand input type {type(inputobj)}") 1497 return 1498 1499 self.properties.SetColor(0.352, 0.612, 0.996) # blue7 1500 1501 self.pipeline = utils.OperationNode( 1502 self, comment=f"#tets {self.dataset.GetNumberOfCells()}", c="#9e2a2b" 1503 )
A RectilinearGrid is a dataset where edges are parallel to the coordinate axes.
It can be thought of as a tessellation of a box in 3D space, similar to a Volume
except that the cells are not necessarily cubes, but they can have different lengths
along each axis.
This can be useful to describe a volume with variable resolution where one needs
to represent a region with higher detail with respect to another region.
Arguments:
- inputobj : (vtkRectilinearGrid, list, str) list of points and tet indices, or filename
Example:
from vedo import RectilinearGrid, show xcoords = 7 + np.sqrt(np.arange(0,2500,25)) ycoords = np.arange(0, 20) zcoords = np.arange(0, 20) rgrid = RectilinearGrid([xcoords, ycoords, zcoords]) print(rgrid) print(rgrid.x_coordinates().shape) print(rgrid.compute_structured_coords([20,10,11])) msh = rgrid.tomesh().lw(1) show(msh, axes=1, viewup="z")
actor
1505 @property 1506 def actor(self): 1507 """Return the `vtkActor` of the object.""" 1508 gf = vtki.new("GeometryFilter") 1509 gf.SetInputData(self.dataset) 1510 gf.Update() 1511 self.mapper.SetInputData(gf.GetOutput()) 1512 self.mapper.Modified() 1513 return self._actor
Return the vtkActor of the object.
def
dimensions(self) -> numpy.ndarray:
1669 def dimensions(self) -> np.ndarray: 1670 """Return the number of points in the x, y and z directions.""" 1671 return np.array(self.dataset.GetDimensions())
Return the number of points in the x, y and z directions.
def
x_coordinates(self) -> numpy.ndarray:
1673 def x_coordinates(self) -> np.ndarray: 1674 """Return the x-coordinates of the grid.""" 1675 return utils.vtk2numpy(self.dataset.GetXCoordinates())
Return the x-coordinates of the grid.
def
y_coordinates(self) -> numpy.ndarray:
1677 def y_coordinates(self) -> np.ndarray: 1678 """Return the y-coordinates of the grid.""" 1679 return utils.vtk2numpy(self.dataset.GetYCoordinates())
Return the y-coordinates of the grid.
def
z_coordinates(self) -> numpy.ndarray:
1681 def z_coordinates(self) -> np.ndarray: 1682 """Return the z-coordinates of the grid.""" 1683 return utils.vtk2numpy(self.dataset.GetZCoordinates())
Return the z-coordinates of the grid.
def
is_point_visible(self, pid: int) -> bool:
1685 def is_point_visible(self, pid: int) -> bool: 1686 """Return True if point `pid` is visible.""" 1687 return self.dataset.IsPointVisible(pid)
Return True if point pid is visible.
def
is_cell_visible(self, cid: int) -> bool:
1689 def is_cell_visible(self, cid: int) -> bool: 1690 """Return True if cell `cid` is visible.""" 1691 return self.dataset.IsCellVisible(cid)
Return True if cell cid is visible.
def
has_blank_points(self) -> bool:
1693 def has_blank_points(self) -> bool: 1694 """Return True if the grid has blank points.""" 1695 return self.dataset.HasAnyBlankPoints()
Return True if the grid has blank points.
def
has_blank_cells(self) -> bool:
1697 def has_blank_cells(self) -> bool: 1698 """Return True if the grid has blank cells.""" 1699 return self.dataset.HasAnyBlankCells()
Return True if the grid has blank cells.
def
compute_structured_coords(self, x: list) -> dict:
1701 def compute_structured_coords(self, x: list) -> dict: 1702 """ 1703 Convenience function computes the structured coordinates for a point `x`. 1704 1705 This method returns a dictionary with keys `ijk`, `pcoords` and `inside`. 1706 The cell is specified by the array `ijk`. 1707 and the parametric coordinates in the cell are specified with `pcoords`. 1708 Value of `inside` is False if the point x is outside of the grid. 1709 """ 1710 ijk = [0, 0, 0] 1711 pcoords = [0., 0., 0.] 1712 inout = self.dataset.ComputeStructuredCoordinates(x, ijk, pcoords) 1713 return {"ijk": np.array(ijk), "pcoords": np.array(pcoords), "inside": bool(inout)}
Convenience function computes the structured coordinates for a point x.
This method returns a dictionary with keys ijk, pcoords and inside.
The cell is specified by the array ijk.
and the parametric coordinates in the cell are specified with pcoords.
Value of inside is False if the point x is outside of the grid.
def
compute_pointid(self, ijk: int) -> int:
1715 def compute_pointid(self, ijk: int) -> int: 1716 """Given a location in structured coordinates (i-j-k), return the point id.""" 1717 return self.dataset.ComputePointId(ijk)
Given a location in structured coordinates (i-j-k), return the point id.
def
compute_cellid(self, ijk: int) -> int:
1719 def compute_cellid(self, ijk: int) -> int: 1720 """Given a location in structured coordinates (i-j-k), return the cell id.""" 1721 return self.dataset.ComputeCellId(ijk)
Given a location in structured coordinates (i-j-k), return the cell id.
def
find_point(self, x: list) -> int:
1723 def find_point(self, x: list) -> int: 1724 """Given a position `x`, return the id of the closest point.""" 1725 return self.dataset.FindPoint(x)
Given a position x, return the id of the closest point.
def
find_cell(self, x: list) -> dict:
1727 def find_cell(self, x: list) -> dict: 1728 """Given a position `x`, return the id of the closest cell.""" 1729 cell = vtki.vtkHexagonalPrism() 1730 cellid = vtki.mutable(0) 1731 tol2 = 0.001 # vtki.mutable(0) 1732 subid = vtki.mutable(0) 1733 pcoords = [0.0, 0.0, 0.0] 1734 weights = [0.0, 0.0, 0.0] 1735 res = self.dataset.FindCell(x, cell, cellid, tol2, subid, pcoords, weights) 1736 result = {} 1737 result["cellid"] = cellid 1738 result["subid"] = subid 1739 result["pcoords"] = pcoords 1740 result["weights"] = weights 1741 result["status"] = res 1742 return result
Given a position x, return the id of the closest cell.
1744 def clone(self, deep=True) -> "RectilinearGrid": 1745 """Return a clone copy of the RectilinearGrid. Alias of `copy()`.""" 1746 if deep: 1747 newrg = vtki.vtkRectilinearGrid() 1748 newrg.CopyStructure(self.dataset) 1749 newrg.CopyAttributes(self.dataset) 1750 newvol = RectilinearGrid(newrg) 1751 else: 1752 newvol = RectilinearGrid(self.dataset) 1753 1754 prop = vtki.vtkProperty() 1755 prop.DeepCopy(self.properties) 1756 newvol.actor.SetProperty(prop) 1757 newvol.properties = prop 1758 newvol.pipeline = utils.OperationNode("clone", parents=[self], c="#bbd0ff", shape="diamond") 1759 return newvol
Return a clone copy of the RectilinearGrid. Alias of copy().
def
bounds(self) -> numpy.ndarray:
1761 def bounds(self) -> np.ndarray: 1762 """ 1763 Get the object bounds. 1764 Returns a list in format `[xmin,xmax, ymin,ymax, zmin,zmax]`. 1765 """ 1766 # OVERRIDE CommonAlgorithms.bounds() which is too slow 1767 return np.array(self.dataset.GetBounds())
Get the object bounds.
Returns a list in format [xmin,xmax, ymin,ymax, zmin,zmax].
1769 def isosurface(self, value=None) -> "vedo.Mesh": 1770 """ 1771 Return a `Mesh` isosurface extracted from the object. 1772 1773 Set `value` as single float or list of values to draw the isosurface(s). 1774 """ 1775 scrange = self.dataset.GetScalarRange() 1776 1777 cf = vtki.new("ContourFilter") 1778 cf.UseScalarTreeOn() 1779 cf.SetInputData(self.dataset) 1780 cf.ComputeNormalsOn() 1781 1782 if value is None: 1783 value = (2 * scrange[0] + scrange[1]) / 3.0 1784 # print("automatic isosurface value =", value) 1785 cf.SetValue(0, value) 1786 else: 1787 if utils.is_sequence(value): 1788 cf.SetNumberOfContours(len(value)) 1789 for i, t in enumerate(value): 1790 cf.SetValue(i, t) 1791 else: 1792 cf.SetValue(0, value) 1793 1794 cf.Update() 1795 poly = cf.GetOutput() 1796 1797 out = vedo.mesh.Mesh(poly, c=None).phong() 1798 out.mapper.SetScalarRange(scrange[0], scrange[1]) 1799 1800 out.pipeline = utils.OperationNode( 1801 "isosurface", 1802 parents=[self], 1803 comment=f"#pts {out.dataset.GetNumberOfPoints()}", 1804 c="#4cc9f0:#e9c46a", 1805 ) 1806 return out
Return a Mesh isosurface extracted from the object.
Set value as single float or list of values to draw the isosurface(s).
1808 def cut_with_plane(self, origin=(0, 0, 0), normal="x") -> "UnstructuredGrid": 1809 """ 1810 Cut the object with the plane defined by a point and a normal. 1811 1812 Arguments: 1813 origin : (list) 1814 the cutting plane goes through this point 1815 normal : (list, str) 1816 normal vector to the cutting plane 1817 """ 1818 strn = str(normal) 1819 if strn == "x": normal = (1, 0, 0) 1820 elif strn == "y": normal = (0, 1, 0) 1821 elif strn == "z": normal = (0, 0, 1) 1822 elif strn == "-x": normal = (-1, 0, 0) 1823 elif strn == "-y": normal = (0, -1, 0) 1824 elif strn == "-z": normal = (0, 0, -1) 1825 plane = vtki.new("Plane") 1826 plane.SetOrigin(origin) 1827 plane.SetNormal(normal) 1828 clipper = vtki.new("ClipDataSet") 1829 clipper.SetInputData(self.dataset) 1830 clipper.SetClipFunction(plane) 1831 clipper.GenerateClipScalarsOff() 1832 clipper.GenerateClippedOutputOff() 1833 clipper.SetValue(0) 1834 clipper.Update() 1835 cout = clipper.GetOutput() 1836 ug = vedo.UnstructuredGrid(cout) 1837 if isinstance(self, UnstructuredGrid): 1838 self._update(cout) 1839 self.pipeline = utils.OperationNode("cut_with_plane", parents=[self], c="#9e2a2b") 1840 return self 1841 ug.pipeline = utils.OperationNode("cut_with_plane", parents=[self], c="#9e2a2b") 1842 return ug
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
def
cut_with_mesh( self, mesh, invert=False, whole_cells=False, on_boundary=False) -> UnstructuredGrid:
1844 def cut_with_mesh(self, mesh, invert=False, whole_cells=False, on_boundary=False) -> "UnstructuredGrid": 1845 """ 1846 Cut a `RectilinearGrid` with a `Mesh`. 1847 1848 Use `invert` to return cut off part of the input object. 1849 """ 1850 ug = self.dataset 1851 1852 ippd = vtki.new("ImplicitPolyDataDistance") 1853 ippd.SetInput(mesh.dataset) 1854 1855 if whole_cells or on_boundary: 1856 clipper = vtki.new("ExtractGeometry") 1857 clipper.SetInputData(ug) 1858 clipper.SetImplicitFunction(ippd) 1859 clipper.SetExtractInside(not invert) 1860 clipper.SetExtractBoundaryCells(False) 1861 if on_boundary: 1862 clipper.SetExtractBoundaryCells(True) 1863 clipper.SetExtractOnlyBoundaryCells(True) 1864 else: 1865 signed_dists = vtki.vtkFloatArray() 1866 signed_dists.SetNumberOfComponents(1) 1867 signed_dists.SetName("SignedDistance") 1868 for pointId in range(ug.GetNumberOfPoints()): 1869 p = ug.GetPoint(pointId) 1870 signed_dist = ippd.EvaluateFunction(p) 1871 signed_dists.InsertNextValue(signed_dist) 1872 ug.GetPointData().AddArray(signed_dists) 1873 ug.GetPointData().SetActiveScalars("SignedDistance") # NEEDED 1874 clipper = vtki.new("ClipDataSet") 1875 clipper.SetInputData(ug) 1876 clipper.SetInsideOut(not invert) 1877 clipper.SetValue(0.0) 1878 1879 clipper.Update() 1880 1881 out = UnstructuredGrid(clipper.GetOutput()) 1882 out.pipeline = utils.OperationNode("cut_with_mesh", parents=[self], c="#9e2a2b") 1883 return out
Cut a RectilinearGrid 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
- update_dataset
- 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
- keep_cell_types
- 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
- generate_random_data
- integrate_data
- write
- tomesh
- signed_distance
- unsigned_distance
- smooth_data
- compute_streamlines
- 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
1887class StructuredGrid(PointAlgorithms, MeshVisual): 1888 """ 1889 Build a structured grid. 1890 """ 1891 1892 def __init__(self, inputobj=None): 1893 """ 1894 A StructuredGrid is a dataset where edges of the hexahedrons are 1895 not necessarily parallel to the coordinate axes. 1896 It can be thought of as a tessellation of a block of 3D space, 1897 similar to a `RectilinearGrid` 1898 except that the cells are not necessarily cubes, they can have different 1899 orientations but are connected in the same way as a `RectilinearGrid`. 1900 1901 Arguments: 1902 inputobj : (vtkStructuredGrid, list, str) 1903 list of points and tet indices, or filename 1904 1905 Example: 1906 ```python 1907 from vedo import * 1908 sgrid = StructuredGrid(dataurl+"structgrid.vts") 1909 print(sgrid) 1910 msh = sgrid.tomesh().lw(1) 1911 show(msh, axes=1, viewup="z") 1912 ``` 1913 1914 ```python 1915 from vedo import * 1916 1917 cx = np.sqrt(np.linspace(100, 400, 10)) 1918 cy = np.linspace(30, 40, 20) 1919 cz = np.linspace(40, 50, 30) 1920 x, y, z = np.meshgrid(cx, cy, cz) 1921 1922 sgrid1 = StructuredGrid([x, y, z]) 1923 sgrid1.cmap("viridis", sgrid1.vertices[:, 0]) 1924 print(sgrid1) 1925 1926 sgrid2 = sgrid1.clone().cut_with_plane(normal=(-1,1,1), origin=[14,34,44]) 1927 msh2 = sgrid2.tomesh(shrink=0.9).lw(1).cmap("viridis") 1928 1929 show( 1930 [["StructuredGrid", sgrid1], ["Shrinked Mesh", msh2]], 1931 N=2, axes=1, viewup="z", 1932 ) 1933 ``` 1934 """ 1935 1936 super().__init__() 1937 1938 self.dataset = None 1939 1940 self.mapper = vtki.new("PolyDataMapper") 1941 self._actor = vtki.vtkActor() 1942 self._actor.retrieve_object = weak_ref_to(self) 1943 self._actor.SetMapper(self.mapper) 1944 self.properties = self._actor.GetProperty() 1945 1946 self.transform = LinearTransform() 1947 self.point_locator = None 1948 self.cell_locator = None 1949 self.line_locator = None 1950 1951 self.name = "StructuredGrid" 1952 self.filename = "" 1953 1954 self.info = {} 1955 self.time = time.time() 1956 1957 ############################### 1958 if inputobj is None: 1959 self.dataset = vtki.vtkStructuredGrid() 1960 1961 elif isinstance(inputobj, vtki.vtkStructuredGrid): 1962 self.dataset = inputobj 1963 1964 elif isinstance(inputobj, StructuredGrid): 1965 self.dataset = inputobj.dataset 1966 1967 elif isinstance(inputobj, str): 1968 if "https://" in inputobj: 1969 inputobj = download(inputobj, verbose=False) 1970 if inputobj.endswith(".vts"): 1971 reader = vtki.new("XMLStructuredGridReader") 1972 else: 1973 reader = vtki.new("StructuredGridReader") 1974 self.filename = inputobj 1975 reader.SetFileName(inputobj) 1976 reader.Update() 1977 self.dataset = reader.GetOutput() 1978 1979 elif utils.is_sequence(inputobj): 1980 self.dataset = vtki.vtkStructuredGrid() 1981 x, y, z = inputobj 1982 xyz = np.vstack(( 1983 x.flatten(order="F"), 1984 y.flatten(order="F"), 1985 z.flatten(order="F")) 1986 ).T 1987 dims = x.shape 1988 self.dataset.SetDimensions(dims) 1989 # self.dataset.SetDimensions(dims[1], dims[0], dims[2]) 1990 vpoints = vtki.vtkPoints() 1991 vpoints.SetData(utils.numpy2vtk(xyz)) 1992 self.dataset.SetPoints(vpoints) 1993 1994 1995 ############################### 1996 if not self.dataset: 1997 vedo.logger.error(f"StructuredGrid: cannot understand input type {type(inputobj)}") 1998 return 1999 2000 self.properties.SetColor(0.352, 0.612, 0.996) # blue7 2001 2002 self.pipeline = utils.OperationNode( 2003 self, comment=f"#tets {self.dataset.GetNumberOfCells()}", c="#9e2a2b" 2004 ) 2005 2006 @property 2007 def actor(self): 2008 """Return the `vtkActor` of the object.""" 2009 gf = vtki.new("GeometryFilter") 2010 gf.SetInputData(self.dataset) 2011 gf.Update() 2012 self.mapper.SetInputData(gf.GetOutput()) 2013 self.mapper.Modified() 2014 return self._actor 2015 2016 @actor.setter 2017 def actor(self, _): 2018 pass 2019 2020 def _update(self, data, reset_locators=False): 2021 self.dataset = data 2022 if reset_locators: 2023 self.cell_locator = None 2024 self.point_locator = None 2025 return self 2026 2027 ################################################################## 2028 def __str__(self): 2029 """Print a summary for the `StructuredGrid` object.""" 2030 module = self.__class__.__module__ 2031 name = self.__class__.__name__ 2032 out = vedo.printc( 2033 f"{module}.{name} at ({hex(self.memory_address())})".ljust(75), 2034 c="c", bold=True, invert=True, return_string=True, 2035 ) 2036 out += "\x1b[0m\x1b[36;1m" 2037 2038 out += "name".ljust(14) + ": " + str(self.name) + "\n" 2039 if self.filename: 2040 out += "filename".ljust(14) + ": " + str(self.filename) + "\n" 2041 2042 out += "dimensions".ljust(14) + ": " + str(self.dataset.GetDimensions()) + "\n" 2043 2044 out += "center".ljust(14) + ": " 2045 out += utils.precision(self.dataset.GetCenter(), 6) + "\n" 2046 2047 bnds = self.bounds() 2048 bx1, bx2 = utils.precision(bnds[0], 3), utils.precision(bnds[1], 3) 2049 by1, by2 = utils.precision(bnds[2], 3), utils.precision(bnds[3], 3) 2050 bz1, bz2 = utils.precision(bnds[4], 3), utils.precision(bnds[5], 3) 2051 out += "bounds".ljust(14) + ":" 2052 out += " x=(" + bx1 + ", " + bx2 + ")," 2053 out += " y=(" + by1 + ", " + by2 + ")," 2054 out += " z=(" + bz1 + ", " + bz2 + ")\n" 2055 2056 out += "memory size".ljust(14) + ": " 2057 out += utils.precision(self.dataset.GetActualMemorySize() / 1024, 2) + " MB\n" 2058 2059 for key in self.pointdata.keys(): 2060 arr = self.pointdata[key] 2061 rng = utils.precision(arr.min(), 3) + ", " + utils.precision(arr.max(), 3) 2062 mark_active = "pointdata" 2063 a_scalars = self.dataset.GetPointData().GetScalars() 2064 a_vectors = self.dataset.GetPointData().GetVectors() 2065 a_tensors = self.dataset.GetPointData().GetTensors() 2066 if a_scalars and a_scalars.GetName() == key: 2067 mark_active += " *" 2068 elif a_vectors and a_vectors.GetName() == key: 2069 mark_active += " **" 2070 elif a_tensors and a_tensors.GetName() == key: 2071 mark_active += " ***" 2072 out += mark_active.ljust(14) + f': "{key}" ({arr.dtype}), ndim={arr.ndim}' 2073 out += f", range=({rng})\n" 2074 2075 for key in self.celldata.keys(): 2076 arr = self.celldata[key] 2077 rng = utils.precision(arr.min(), 3) + ", " + utils.precision(arr.max(), 3) 2078 mark_active = "celldata" 2079 a_scalars = self.dataset.GetCellData().GetScalars() 2080 a_vectors = self.dataset.GetCellData().GetVectors() 2081 a_tensors = self.dataset.GetCellData().GetTensors() 2082 if a_scalars and a_scalars.GetName() == key: 2083 mark_active += " *" 2084 elif a_vectors and a_vectors.GetName() == key: 2085 mark_active += " **" 2086 elif a_tensors and a_tensors.GetName() == key: 2087 mark_active += " ***" 2088 out += mark_active.ljust(14) + f': "{key}" ({arr.dtype}), ndim={arr.ndim}' 2089 out += f", range=({rng})\n" 2090 2091 for key in self.metadata.keys(): 2092 arr = self.metadata[key] 2093 out += "metadata".ljust(14) + ": " + f'"{key}" ({len(arr)} values)\n' 2094 2095 return out.rstrip() + "\x1b[0m" 2096 2097 def _repr_html_(self): 2098 """ 2099 HTML representation of the StructuredGrid object for Jupyter Notebooks. 2100 2101 Returns: 2102 HTML text with the image and some properties. 2103 """ 2104 import io 2105 import base64 2106 from PIL import Image 2107 2108 library_name = "vedo.grids.StructuredGrid" 2109 help_url = "https://vedo.embl.es/docs/vedo/grids.html#StructuredGrid" 2110 2111 m = self.tomesh().linewidth(1).lighting("off") 2112 arr= m.thumbnail(zoom=1, elevation=-30, azimuth=-30) 2113 2114 im = Image.fromarray(arr) 2115 buffered = io.BytesIO() 2116 im.save(buffered, format="PNG", quality=100) 2117 encoded = base64.b64encode(buffered.getvalue()).decode("utf-8") 2118 url = "data:image/png;base64," + encoded 2119 image = f"<img src='{url}'></img>" 2120 2121 bounds = "<br/>".join( 2122 [ 2123 utils.precision(min_x, 4) + " ... " + utils.precision(max_x, 4) 2124 for min_x, max_x in zip(self.bounds()[::2], self.bounds()[1::2]) 2125 ] 2126 ) 2127 2128 help_text = "" 2129 if self.name: 2130 help_text += f"<b> {self.name}:   </b>" 2131 help_text += '<b><a href="' + help_url + '" target="_blank">' + library_name + "</a></b>" 2132 if self.filename: 2133 dots = "" 2134 if len(self.filename) > 30: 2135 dots = "..." 2136 help_text += f"<br/><code><i>({dots}{self.filename[-30:]})</i></code>" 2137 2138 pdata = "" 2139 if self.dataset.GetPointData().GetScalars(): 2140 if self.dataset.GetPointData().GetScalars().GetName(): 2141 name = self.dataset.GetPointData().GetScalars().GetName() 2142 pdata = "<tr><td><b> point data array </b></td><td>" + name + "</td></tr>" 2143 2144 cdata = "" 2145 if self.dataset.GetCellData().GetScalars(): 2146 if self.dataset.GetCellData().GetScalars().GetName(): 2147 name = self.dataset.GetCellData().GetScalars().GetName() 2148 cdata = "<tr><td><b> cell data array </b></td><td>" + name + "</td></tr>" 2149 2150 pts = self.vertices 2151 cm = np.mean(pts, axis=0) 2152 2153 all = [ 2154 "<table>", 2155 "<tr>", 2156 "<td>", image, "</td>", 2157 "<td style='text-align: center; vertical-align: center;'><br/>", help_text, 2158 "<table>", 2159 "<tr><td><b> bounds </b> <br/> (x/y/z) </td><td>" + str(bounds) + "</td></tr>", 2160 "<tr><td><b> center of mass </b></td><td>" + utils.precision(cm,3) + "</td></tr>", 2161 "<tr><td><b> nr. points / cells </b></td><td>" 2162 + str(self.npoints) + " / " + str(self.ncells) + "</td></tr>", 2163 pdata, 2164 cdata, 2165 "</table>", 2166 "</table>", 2167 ] 2168 return "\n".join(all) 2169 2170 def dimensions(self) -> np.ndarray: 2171 """Return the number of points in the x, y and z directions.""" 2172 return np.array(self.dataset.GetDimensions()) 2173 2174 def clone(self, deep=True) -> "StructuredGrid": 2175 """Return a clone copy of the StructuredGrid. Alias of `copy()`.""" 2176 if deep: 2177 newrg = vtki.vtkStructuredGrid() 2178 newrg.CopyStructure(self.dataset) 2179 newrg.CopyAttributes(self.dataset) 2180 newvol = StructuredGrid(newrg) 2181 else: 2182 newvol = StructuredGrid(self.dataset) 2183 2184 prop = vtki.vtkProperty() 2185 prop.DeepCopy(self.properties) 2186 newvol.actor.SetProperty(prop) 2187 newvol.properties = prop 2188 newvol.pipeline = utils.OperationNode("clone", parents=[self], c="#bbd0ff", shape="diamond") 2189 return newvol 2190 2191 def find_point(self, x: list) -> int: 2192 """Given a position `x`, return the id of the closest point.""" 2193 return self.dataset.FindPoint(x) 2194 2195 def find_cell(self, x: list) -> dict: 2196 """Given a position `x`, return the id of the closest cell.""" 2197 cell = vtki.vtkHexagonalPrism() 2198 cellid = vtki.mutable(0) 2199 tol2 = 0.001 # vtki.mutable(0) 2200 subid = vtki.mutable(0) 2201 pcoords = [0.0, 0.0, 0.0] 2202 weights = [0.0, 0.0, 0.0] 2203 res = self.dataset.FindCell(x, cell, cellid, tol2, subid, pcoords, weights) 2204 result = {} 2205 result["cellid"] = cellid 2206 result["subid"] = subid 2207 result["pcoords"] = pcoords 2208 result["weights"] = weights 2209 result["status"] = res 2210 return result 2211 2212 def cut_with_plane(self, origin=(0, 0, 0), normal="x") -> "vedo.UnstructuredGrid": 2213 """ 2214 Cut the object with the plane defined by a point and a normal. 2215 2216 Arguments: 2217 origin : (list) 2218 the cutting plane goes through this point 2219 normal : (list, str) 2220 normal vector to the cutting plane 2221 2222 Returns an `UnstructuredGrid` object. 2223 """ 2224 strn = str(normal) 2225 if strn == "x": normal = (1, 0, 0) 2226 elif strn == "y": normal = (0, 1, 0) 2227 elif strn == "z": normal = (0, 0, 1) 2228 elif strn == "-x": normal = (-1, 0, 0) 2229 elif strn == "-y": normal = (0, -1, 0) 2230 elif strn == "-z": normal = (0, 0, -1) 2231 plane = vtki.new("Plane") 2232 plane.SetOrigin(origin) 2233 plane.SetNormal(normal) 2234 clipper = vtki.new("ClipDataSet") 2235 clipper.SetInputData(self.dataset) 2236 clipper.SetClipFunction(plane) 2237 clipper.GenerateClipScalarsOff() 2238 clipper.GenerateClippedOutputOff() 2239 clipper.SetValue(0) 2240 clipper.Update() 2241 cout = clipper.GetOutput() 2242 ug = vedo.UnstructuredGrid(cout) 2243 if isinstance(self, vedo.UnstructuredGrid): 2244 self._update(cout) 2245 self.pipeline = utils.OperationNode("cut_with_plane", parents=[self], c="#9e2a2b") 2246 return self 2247 ug.pipeline = utils.OperationNode("cut_with_plane", parents=[self], c="#9e2a2b") 2248 return ug 2249 2250 def cut_with_mesh(self, mesh: Mesh, invert=False, whole_cells=False, on_boundary=False) -> "UnstructuredGrid": 2251 """ 2252 Cut a `RectilinearGrid` with a `Mesh`. 2253 2254 Use `invert` to return cut off part of the input object. 2255 2256 Returns an `UnstructuredGrid` object. 2257 """ 2258 ug = self.dataset 2259 2260 ippd = vtki.new("ImplicitPolyDataDistance") 2261 ippd.SetInput(mesh.dataset) 2262 2263 if whole_cells or on_boundary: 2264 clipper = vtki.new("ExtractGeometry") 2265 clipper.SetInputData(ug) 2266 clipper.SetImplicitFunction(ippd) 2267 clipper.SetExtractInside(not invert) 2268 clipper.SetExtractBoundaryCells(False) 2269 if on_boundary: 2270 clipper.SetExtractBoundaryCells(True) 2271 clipper.SetExtractOnlyBoundaryCells(True) 2272 else: 2273 signed_dists = vtki.vtkFloatArray() 2274 signed_dists.SetNumberOfComponents(1) 2275 signed_dists.SetName("SignedDistance") 2276 for pointId in range(ug.GetNumberOfPoints()): 2277 p = ug.GetPoint(pointId) 2278 signed_dist = ippd.EvaluateFunction(p) 2279 signed_dists.InsertNextValue(signed_dist) 2280 ug.GetPointData().AddArray(signed_dists) 2281 ug.GetPointData().SetActiveScalars("SignedDistance") # NEEDED 2282 clipper = vtki.new("ClipDataSet") 2283 clipper.SetInputData(ug) 2284 clipper.SetInsideOut(not invert) 2285 clipper.SetValue(0.0) 2286 2287 clipper.Update() 2288 2289 out = UnstructuredGrid(clipper.GetOutput()) 2290 out.pipeline = utils.OperationNode("cut_with_mesh", parents=[self], c="#9e2a2b") 2291 return out 2292 2293 def isosurface(self, value=None) -> "vedo.Mesh": 2294 """ 2295 Return a `Mesh` isosurface extracted from the object. 2296 2297 Set `value` as single float or list of values to draw the isosurface(s). 2298 """ 2299 scrange = self.dataset.GetScalarRange() 2300 2301 cf = vtki.new("ContourFilter") 2302 cf.UseScalarTreeOn() 2303 cf.SetInputData(self.dataset) 2304 cf.ComputeNormalsOn() 2305 2306 if value is None: 2307 value = (2 * scrange[0] + scrange[1]) / 3.0 2308 # print("automatic isosurface value =", value) 2309 cf.SetValue(0, value) 2310 else: 2311 if utils.is_sequence(value): 2312 cf.SetNumberOfContours(len(value)) 2313 for i, t in enumerate(value): 2314 cf.SetValue(i, t) 2315 else: 2316 cf.SetValue(0, value) 2317 2318 cf.Update() 2319 poly = cf.GetOutput() 2320 2321 out = vedo.mesh.Mesh(poly, c=None).phong() 2322 out.mapper.SetScalarRange(scrange[0], scrange[1]) 2323 2324 out.pipeline = utils.OperationNode( 2325 "isosurface", 2326 parents=[self], 2327 comment=f"#pts {out.dataset.GetNumberOfPoints()}", 2328 c="#4cc9f0:#e9c46a", 2329 ) 2330 return out
Build a structured grid.
StructuredGrid(inputobj=None)
1892 def __init__(self, inputobj=None): 1893 """ 1894 A StructuredGrid is a dataset where edges of the hexahedrons are 1895 not necessarily parallel to the coordinate axes. 1896 It can be thought of as a tessellation of a block of 3D space, 1897 similar to a `RectilinearGrid` 1898 except that the cells are not necessarily cubes, they can have different 1899 orientations but are connected in the same way as a `RectilinearGrid`. 1900 1901 Arguments: 1902 inputobj : (vtkStructuredGrid, list, str) 1903 list of points and tet indices, or filename 1904 1905 Example: 1906 ```python 1907 from vedo import * 1908 sgrid = StructuredGrid(dataurl+"structgrid.vts") 1909 print(sgrid) 1910 msh = sgrid.tomesh().lw(1) 1911 show(msh, axes=1, viewup="z") 1912 ``` 1913 1914 ```python 1915 from vedo import * 1916 1917 cx = np.sqrt(np.linspace(100, 400, 10)) 1918 cy = np.linspace(30, 40, 20) 1919 cz = np.linspace(40, 50, 30) 1920 x, y, z = np.meshgrid(cx, cy, cz) 1921 1922 sgrid1 = StructuredGrid([x, y, z]) 1923 sgrid1.cmap("viridis", sgrid1.vertices[:, 0]) 1924 print(sgrid1) 1925 1926 sgrid2 = sgrid1.clone().cut_with_plane(normal=(-1,1,1), origin=[14,34,44]) 1927 msh2 = sgrid2.tomesh(shrink=0.9).lw(1).cmap("viridis") 1928 1929 show( 1930 [["StructuredGrid", sgrid1], ["Shrinked Mesh", msh2]], 1931 N=2, axes=1, viewup="z", 1932 ) 1933 ``` 1934 """ 1935 1936 super().__init__() 1937 1938 self.dataset = None 1939 1940 self.mapper = vtki.new("PolyDataMapper") 1941 self._actor = vtki.vtkActor() 1942 self._actor.retrieve_object = weak_ref_to(self) 1943 self._actor.SetMapper(self.mapper) 1944 self.properties = self._actor.GetProperty() 1945 1946 self.transform = LinearTransform() 1947 self.point_locator = None 1948 self.cell_locator = None 1949 self.line_locator = None 1950 1951 self.name = "StructuredGrid" 1952 self.filename = "" 1953 1954 self.info = {} 1955 self.time = time.time() 1956 1957 ############################### 1958 if inputobj is None: 1959 self.dataset = vtki.vtkStructuredGrid() 1960 1961 elif isinstance(inputobj, vtki.vtkStructuredGrid): 1962 self.dataset = inputobj 1963 1964 elif isinstance(inputobj, StructuredGrid): 1965 self.dataset = inputobj.dataset 1966 1967 elif isinstance(inputobj, str): 1968 if "https://" in inputobj: 1969 inputobj = download(inputobj, verbose=False) 1970 if inputobj.endswith(".vts"): 1971 reader = vtki.new("XMLStructuredGridReader") 1972 else: 1973 reader = vtki.new("StructuredGridReader") 1974 self.filename = inputobj 1975 reader.SetFileName(inputobj) 1976 reader.Update() 1977 self.dataset = reader.GetOutput() 1978 1979 elif utils.is_sequence(inputobj): 1980 self.dataset = vtki.vtkStructuredGrid() 1981 x, y, z = inputobj 1982 xyz = np.vstack(( 1983 x.flatten(order="F"), 1984 y.flatten(order="F"), 1985 z.flatten(order="F")) 1986 ).T 1987 dims = x.shape 1988 self.dataset.SetDimensions(dims) 1989 # self.dataset.SetDimensions(dims[1], dims[0], dims[2]) 1990 vpoints = vtki.vtkPoints() 1991 vpoints.SetData(utils.numpy2vtk(xyz)) 1992 self.dataset.SetPoints(vpoints) 1993 1994 1995 ############################### 1996 if not self.dataset: 1997 vedo.logger.error(f"StructuredGrid: cannot understand input type {type(inputobj)}") 1998 return 1999 2000 self.properties.SetColor(0.352, 0.612, 0.996) # blue7 2001 2002 self.pipeline = utils.OperationNode( 2003 self, comment=f"#tets {self.dataset.GetNumberOfCells()}", c="#9e2a2b" 2004 )
A StructuredGrid is a dataset where edges of the hexahedrons are
not necessarily parallel to the coordinate axes.
It can be thought of as a tessellation of a block of 3D space,
similar to a RectilinearGrid
except that the cells are not necessarily cubes, they can have different
orientations but are connected in the same way as a RectilinearGrid.
Arguments:
- inputobj : (vtkStructuredGrid, list, str) list of points and tet indices, or filename
Example:
from vedo import * sgrid = StructuredGrid(dataurl+"structgrid.vts") print(sgrid) msh = sgrid.tomesh().lw(1) show(msh, axes=1, viewup="z")from vedo import * cx = np.sqrt(np.linspace(100, 400, 10)) cy = np.linspace(30, 40, 20) cz = np.linspace(40, 50, 30) x, y, z = np.meshgrid(cx, cy, cz) sgrid1 = StructuredGrid([x, y, z]) sgrid1.cmap("viridis", sgrid1.vertices[:, 0]) print(sgrid1) sgrid2 = sgrid1.clone().cut_with_plane(normal=(-1,1,1), origin=[14,34,44]) msh2 = sgrid2.tomesh(shrink=0.9).lw(1).cmap("viridis") show( [["StructuredGrid", sgrid1], ["Shrinked Mesh", msh2]], N=2, axes=1, viewup="z", )
actor
2006 @property 2007 def actor(self): 2008 """Return the `vtkActor` of the object.""" 2009 gf = vtki.new("GeometryFilter") 2010 gf.SetInputData(self.dataset) 2011 gf.Update() 2012 self.mapper.SetInputData(gf.GetOutput()) 2013 self.mapper.Modified() 2014 return self._actor
Return the vtkActor of the object.
def
dimensions(self) -> numpy.ndarray:
2170 def dimensions(self) -> np.ndarray: 2171 """Return the number of points in the x, y and z directions.""" 2172 return np.array(self.dataset.GetDimensions())
Return the number of points in the x, y and z directions.
2174 def clone(self, deep=True) -> "StructuredGrid": 2175 """Return a clone copy of the StructuredGrid. Alias of `copy()`.""" 2176 if deep: 2177 newrg = vtki.vtkStructuredGrid() 2178 newrg.CopyStructure(self.dataset) 2179 newrg.CopyAttributes(self.dataset) 2180 newvol = StructuredGrid(newrg) 2181 else: 2182 newvol = StructuredGrid(self.dataset) 2183 2184 prop = vtki.vtkProperty() 2185 prop.DeepCopy(self.properties) 2186 newvol.actor.SetProperty(prop) 2187 newvol.properties = prop 2188 newvol.pipeline = utils.OperationNode("clone", parents=[self], c="#bbd0ff", shape="diamond") 2189 return newvol
Return a clone copy of the StructuredGrid. Alias of copy().
def
find_point(self, x: list) -> int:
2191 def find_point(self, x: list) -> int: 2192 """Given a position `x`, return the id of the closest point.""" 2193 return self.dataset.FindPoint(x)
Given a position x, return the id of the closest point.
def
find_cell(self, x: list) -> dict:
2195 def find_cell(self, x: list) -> dict: 2196 """Given a position `x`, return the id of the closest cell.""" 2197 cell = vtki.vtkHexagonalPrism() 2198 cellid = vtki.mutable(0) 2199 tol2 = 0.001 # vtki.mutable(0) 2200 subid = vtki.mutable(0) 2201 pcoords = [0.0, 0.0, 0.0] 2202 weights = [0.0, 0.0, 0.0] 2203 res = self.dataset.FindCell(x, cell, cellid, tol2, subid, pcoords, weights) 2204 result = {} 2205 result["cellid"] = cellid 2206 result["subid"] = subid 2207 result["pcoords"] = pcoords 2208 result["weights"] = weights 2209 result["status"] = res 2210 return result
Given a position x, return the id of the closest cell.
2212 def cut_with_plane(self, origin=(0, 0, 0), normal="x") -> "vedo.UnstructuredGrid": 2213 """ 2214 Cut the object with the plane defined by a point and a normal. 2215 2216 Arguments: 2217 origin : (list) 2218 the cutting plane goes through this point 2219 normal : (list, str) 2220 normal vector to the cutting plane 2221 2222 Returns an `UnstructuredGrid` object. 2223 """ 2224 strn = str(normal) 2225 if strn == "x": normal = (1, 0, 0) 2226 elif strn == "y": normal = (0, 1, 0) 2227 elif strn == "z": normal = (0, 0, 1) 2228 elif strn == "-x": normal = (-1, 0, 0) 2229 elif strn == "-y": normal = (0, -1, 0) 2230 elif strn == "-z": normal = (0, 0, -1) 2231 plane = vtki.new("Plane") 2232 plane.SetOrigin(origin) 2233 plane.SetNormal(normal) 2234 clipper = vtki.new("ClipDataSet") 2235 clipper.SetInputData(self.dataset) 2236 clipper.SetClipFunction(plane) 2237 clipper.GenerateClipScalarsOff() 2238 clipper.GenerateClippedOutputOff() 2239 clipper.SetValue(0) 2240 clipper.Update() 2241 cout = clipper.GetOutput() 2242 ug = vedo.UnstructuredGrid(cout) 2243 if isinstance(self, vedo.UnstructuredGrid): 2244 self._update(cout) 2245 self.pipeline = utils.OperationNode("cut_with_plane", parents=[self], c="#9e2a2b") 2246 return self 2247 ug.pipeline = utils.OperationNode("cut_with_plane", parents=[self], c="#9e2a2b") 2248 return ug
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
Returns an UnstructuredGrid object.
def
cut_with_mesh( self, mesh: vedo.mesh.Mesh, invert=False, whole_cells=False, on_boundary=False) -> UnstructuredGrid:
2250 def cut_with_mesh(self, mesh: Mesh, invert=False, whole_cells=False, on_boundary=False) -> "UnstructuredGrid": 2251 """ 2252 Cut a `RectilinearGrid` with a `Mesh`. 2253 2254 Use `invert` to return cut off part of the input object. 2255 2256 Returns an `UnstructuredGrid` object. 2257 """ 2258 ug = self.dataset 2259 2260 ippd = vtki.new("ImplicitPolyDataDistance") 2261 ippd.SetInput(mesh.dataset) 2262 2263 if whole_cells or on_boundary: 2264 clipper = vtki.new("ExtractGeometry") 2265 clipper.SetInputData(ug) 2266 clipper.SetImplicitFunction(ippd) 2267 clipper.SetExtractInside(not invert) 2268 clipper.SetExtractBoundaryCells(False) 2269 if on_boundary: 2270 clipper.SetExtractBoundaryCells(True) 2271 clipper.SetExtractOnlyBoundaryCells(True) 2272 else: 2273 signed_dists = vtki.vtkFloatArray() 2274 signed_dists.SetNumberOfComponents(1) 2275 signed_dists.SetName("SignedDistance") 2276 for pointId in range(ug.GetNumberOfPoints()): 2277 p = ug.GetPoint(pointId) 2278 signed_dist = ippd.EvaluateFunction(p) 2279 signed_dists.InsertNextValue(signed_dist) 2280 ug.GetPointData().AddArray(signed_dists) 2281 ug.GetPointData().SetActiveScalars("SignedDistance") # NEEDED 2282 clipper = vtki.new("ClipDataSet") 2283 clipper.SetInputData(ug) 2284 clipper.SetInsideOut(not invert) 2285 clipper.SetValue(0.0) 2286 2287 clipper.Update() 2288 2289 out = UnstructuredGrid(clipper.GetOutput()) 2290 out.pipeline = utils.OperationNode("cut_with_mesh", parents=[self], c="#9e2a2b") 2291 return out
Cut a RectilinearGrid with a Mesh.
Use invert to return cut off part of the input object.
Returns an UnstructuredGrid object.
2293 def isosurface(self, value=None) -> "vedo.Mesh": 2294 """ 2295 Return a `Mesh` isosurface extracted from the object. 2296 2297 Set `value` as single float or list of values to draw the isosurface(s). 2298 """ 2299 scrange = self.dataset.GetScalarRange() 2300 2301 cf = vtki.new("ContourFilter") 2302 cf.UseScalarTreeOn() 2303 cf.SetInputData(self.dataset) 2304 cf.ComputeNormalsOn() 2305 2306 if value is None: 2307 value = (2 * scrange[0] + scrange[1]) / 3.0 2308 # print("automatic isosurface value =", value) 2309 cf.SetValue(0, value) 2310 else: 2311 if utils.is_sequence(value): 2312 cf.SetNumberOfContours(len(value)) 2313 for i, t in enumerate(value): 2314 cf.SetValue(i, t) 2315 else: 2316 cf.SetValue(0, value) 2317 2318 cf.Update() 2319 poly = cf.GetOutput() 2320 2321 out = vedo.mesh.Mesh(poly, c=None).phong() 2322 out.mapper.SetScalarRange(scrange[0], scrange[1]) 2323 2324 out.pipeline = utils.OperationNode( 2325 "isosurface", 2326 parents=[self], 2327 comment=f"#pts {out.dataset.GetNumberOfPoints()}", 2328 c="#4cc9f0:#e9c46a", 2329 ) 2330 return out
Return a Mesh isosurface extracted from the object.
Set value as single float or list of values to draw the isosurface(s).
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
- update_dataset
- bounds
- 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
- keep_cell_types
- 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
- generate_random_data
- integrate_data
- write
- tomesh
- signed_distance
- unsigned_distance
- smooth_data
- compute_streamlines
- 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