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