vedo.grids

Work with tetrahedral meshes.

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

Support for UnstructuredGrid objects.

UnstructuredGrid(inputobj=None)
 41    def __init__(self, inputobj=None):
 42        """
 43        Support for UnstructuredGrid objects.
 44
 45        Arguments:
 46            inputobj : (list, vtkUnstructuredGrid, str)
 47                A list in the form `[points, cells, celltypes]`,
 48                or a vtkUnstructuredGrid object, or a filename
 49
 50        Celltypes are identified by the following 
 51        [convention](https://vtk.org/doc/nightly/html/vtkCellType_8h_source.html).
 52        """
 53        super().__init__()
 54
 55        self.dataset = None
 56
 57        self.mapper = vtki.new("PolyDataMapper")
 58        self._actor = vtki.vtkActor()
 59        self._actor.retrieve_object = weak_ref_to(self)
 60        self._actor.SetMapper(self.mapper)
 61        self.properties = self._actor.GetProperty()
 62
 63        self.transform = LinearTransform()
 64        self.point_locator = None
 65        self.cell_locator = None
 66        self.line_locator = None
 67
 68        self.name = "UnstructuredGrid"
 69        self.filename = ""
 70        self.file_size = ""
 71
 72        self.info = {}
 73        self.time = time.time()
 74        self.rendered_at = set()
 75
 76        ######################################
 77        inputtype = str(type(inputobj))
 78
 79        if inputobj is None:
 80            self.dataset = vtki.vtkUnstructuredGrid()
 81
 82        elif utils.is_sequence(inputobj):
 83
 84            pts, cells, celltypes = inputobj
 85            assert len(cells) == len(celltypes)
 86
 87            self.dataset = vtki.vtkUnstructuredGrid()
 88
 89            if not utils.is_sequence(cells[0]):
 90                tets = []
 91                nf = cells[0] + 1
 92                for i, cl in enumerate(cells):
 93                    if i in (nf, 0):
 94                        k = i + 1
 95                        nf = cl + k
 96                        cell = [cells[j + k] for j in range(cl)]
 97                        tets.append(cell)
 98                cells = tets
 99
100            # This would fill the points and use those to define orientation
101            vpts = utils.numpy2vtk(pts, dtype=np.float32)
102            points = vtki.vtkPoints()
103            points.SetData(vpts)
104            self.dataset.SetPoints(points)
105
106            # Fill cells
107            # https://vtk.org/doc/nightly/html/vtkCellType_8h_source.html
108            for i, ct in enumerate(celltypes):
109                if   ct == vtki.cell_types["VERTEX"]:
110                    cell = vtki.vtkVertex()
111                elif ct == vtki.cell_types["POLY_VERTEX"]:
112                    cell = vtki.vtkPolyVertex()
113                elif ct == vtki.cell_types["TETRA"]:
114                    cell = vtki.vtkTetra()
115                elif ct == vtki.cell_types["WEDGE"]:
116                    cell = vtki.vtkWedge()
117                elif ct == vtki.cell_types["LINE"]:
118                    cell = vtki.vtkLine()
119                elif ct == vtki.cell_types["POLY_LINE"]:
120                    cell = vtki.vtkPolyLine()
121                elif ct == vtki.cell_types["TRIANGLE"]:
122                    cell = vtki.vtkTriangle()
123                elif ct == vtki.cell_types["TRIANGLE_STRIP"]:
124                    cell = vtki.vtkTriangleStrip()
125                elif ct == vtki.cell_types["POLYGON"]:
126                    cell = vtki.vtkPolygon()
127                elif ct == vtki.cell_types["PIXEL"]:
128                    cell = vtki.vtkPixel()
129                elif ct == vtki.cell_types["QUAD"]:
130                    cell = vtki.vtkQuad()
131                elif ct == vtki.cell_types["VOXEL"]:
132                    cell = vtki.vtkVoxel()
133                elif ct == vtki.cell_types["PYRAMID"]:
134                    cell = vtki.vtkPyramid()
135                elif ct == vtki.cell_types["HEXAHEDRON"]:
136                    cell = vtki.vtkHexahedron()
137                elif ct == vtki.cell_types["HEXAGONAL_PRISM"]:
138                    cell = vtki.vtkHexagonalPrism()
139                elif ct == vtki.cell_types["PENTAGONAL_PRISM"]:
140                    cell = vtki.vtkPentagonalPrism()
141                elif ct == vtki.cell_types["QUADRATIC_TETRA"]:
142                    from vtkmodules.vtkCommonDataModel import vtkQuadraticTetra
143                    cell = vtkQuadraticTetra()
144                elif ct == vtki.cell_types["QUADRATIC_HEXAHEDRON"]:
145                    from vtkmodules.vtkCommonDataModel import vtkQuadraticHexahedron
146                    cell = vtkQuadraticHexahedron()
147                elif ct == vtki.cell_types["QUADRATIC_WEDGE"]:
148                    from vtkmodules.vtkCommonDataModel import vtkQuadraticWedge
149                    cell = vtkQuadraticWedge()
150                elif ct == vtki.cell_types["QUADRATIC_PYRAMID"]:
151                    from vtkmodules.vtkCommonDataModel import vtkQuadraticPyramid
152                    cell = vtkQuadraticPyramid()
153                elif ct == vtki.cell_types["QUADRATIC_LINEAR_QUAD"]:
154                    from vtkmodules.vtkCommonDataModel import vtkQuadraticLinearQuad
155                    cell = vtkQuadraticLinearQuad()
156                elif ct == vtki.cell_types["QUADRATIC_LINEAR_WEDGE"]:
157                    from vtkmodules.vtkCommonDataModel import vtkQuadraticLinearWedge
158                    cell = vtkQuadraticLinearWedge()
159                elif ct == vtki.cell_types["BIQUADRATIC_QUADRATIC_WEDGE"]:
160                    from vtkmodules.vtkCommonDataModel import vtkBiQuadraticQuadraticWedge
161                    cell = vtkBiQuadraticQuadraticWedge()
162                elif ct == vtki.cell_types["BIQUADRATIC_QUADRATIC_HEXAHEDRON"]:
163                    from vtkmodules.vtkCommonDataModel import vtkBiQuadraticQuadraticHexahedron
164                    cell = vtkBiQuadraticQuadraticHexahedron()
165                elif ct == vtki.cell_types["BIQUADRATIC_TRIANGLE"]:
166                    from vtkmodules.vtkCommonDataModel import vtkBiQuadraticTriangle
167                    cell = vtkBiQuadraticTriangle()
168                elif ct == vtki.cell_types["CUBIC_LINE"]:
169                    from vtkmodules.vtkCommonDataModel import vtkCubicLine
170                    cell = vtkCubicLine()
171                elif ct == vtki.cell_types["CONVEX_POINT_SET"]:
172                    from vtkmodules.vtkCommonDataModel import vtkConvexPointSet
173                    cell = vtkConvexPointSet()
174                elif ct == vtki.cell_types["POLYHEDRON"]:
175                    from vtkmodules.vtkCommonDataModel import vtkPolyhedron
176                    cell = vtkPolyhedron()
177                elif ct == vtki.cell_types["HIGHER_ORDER_TRIANGLE"]:
178                    from vtkmodules.vtkCommonDataModel import vtkHigherOrderTriangle
179                    cell = vtkHigherOrderTriangle()
180                elif ct == vtki.cell_types["HIGHER_ORDER_QUAD"]:
181                    from vtkmodules.vtkCommonDataModel import vtkHigherOrderQuadrilateral
182                    cell = vtkHigherOrderQuadrilateral()
183                elif ct == vtki.cell_types["HIGHER_ORDER_TETRAHEDRON"]:
184                    from vtkmodules.vtkCommonDataModel import vtkHigherOrderTetra
185                    cell = vtkHigherOrderTetra()
186                elif ct == vtki.cell_types["HIGHER_ORDER_WEDGE"]:
187                    from vtkmodules.vtkCommonDataModel import vtkHigherOrderWedge
188                    cell = vtkHigherOrderWedge()
189                elif ct == vtki.cell_types["HIGHER_ORDER_HEXAHEDRON"]:
190                    from vtkmodules.vtkCommonDataModel import vtkHigherOrderHexahedron
191                    cell = vtkHigherOrderHexahedron()
192                elif ct == vtki.cell_types["LAGRANGE_CURVE"]:
193                    from vtkmodules.vtkCommonDataModel import vtkLagrangeCurve
194                    cell = vtkLagrangeCurve()
195                elif ct == vtki.cell_types["LAGRANGE_TRIANGLE"]:
196                    from vtkmodules.vtkCommonDataModel import vtkLagrangeTriangle
197                    cell = vtkLagrangeTriangle()
198                elif ct == vtki.cell_types["LAGRANGE_QUADRILATERAL"]:
199                    from vtkmodules.vtkCommonDataModel import vtkLagrangeQuadrilateral
200                    cell = vtkLagrangeQuadrilateral()
201                elif ct == vtki.cell_types["LAGRANGE_TETRAHEDRON"]:
202                    from vtkmodules.vtkCommonDataModel import vtkLagrangeTetra
203                    cell = vtkLagrangeTetra()
204                elif ct == vtki.cell_types["LAGRANGE_HEXAHEDRON"]:
205                    from vtkmodules.vtkCommonDataModel import vtkLagrangeHexahedron
206                    cell = vtkLagrangeHexahedron()
207                elif ct == vtki.cell_types["LAGRANGE_WEDGE"]:
208                    from vtkmodules.vtkCommonDataModel import vtkLagrangeWedge
209                    cell = vtkLagrangeWedge()
210                elif ct == vtki.cell_types["BEZIER_CURVE"]:
211                    from vtkmodules.vtkCommonDataModel import vtkBezierCurve
212                    cell = vtkBezierCurve()
213                elif ct == vtki.cell_types["BEZIER_TRIANGLE"]:
214                    from vtkmodules.vtkCommonDataModel import vtkBezierTriangle
215                    cell = vtkBezierTriangle()
216                elif ct == vtki.cell_types["BEZIER_QUADRILATERAL"]:
217                    from vtkmodules.vtkCommonDataModel import vtkBezierQuadrilateral
218                    cell = vtkBezierQuadrilateral()
219                elif ct == vtki.cell_types["BEZIER_TETRAHEDRON"]:
220                    from vtkmodules.vtkCommonDataModel import vtkBezierTetra
221                    cell = vtkBezierTetra()
222                elif ct == vtki.cell_types["BEZIER_HEXAHEDRON"]:
223                    from vtkmodules.vtkCommonDataModel import vtkBezierHexahedron
224                    cell = vtkBezierHexahedron()
225                elif ct == vtki.cell_types["BEZIER_WEDGE"]:
226                    from vtkmodules.vtkCommonDataModel import vtkBezierWedge
227                    cell = vtkBezierWedge()
228                else:
229                    vedo.logger.error(
230                        f"UnstructuredGrid: cell type {ct} not supported. Skip.")
231                    continue
232
233                cpids = cell.GetPointIds()
234                cell_conn = cells[i]
235                for j, pid in enumerate(cell_conn):
236                    cpids.SetId(j, pid)
237                self.dataset.InsertNextCell(ct, cpids)
238
239        elif "UnstructuredGrid" in inputtype:
240            self.dataset = inputobj
241
242        elif isinstance(inputobj, str):
243            if "https://" in inputobj:
244                inputobj = download(inputobj, verbose=False)
245            self.filename = inputobj
246            if inputobj.endswith(".vtu"):
247                reader = vtki.new("XMLUnstructuredGridReader")
248            else:
249                reader = vtki.new("UnstructuredGridReader")
250            self.filename = inputobj
251            reader.SetFileName(inputobj)
252            reader.Update()
253            self.dataset = reader.GetOutput()
254
255        else:
256            # this converts other types of vtk objects to UnstructuredGrid
257            apf = vtki.new("AppendFilter")
258            try:
259                apf.AddInputData(inputobj)
260            except TypeError:
261                apf.AddInputData(inputobj.dataset)
262            apf.Update()
263            self.dataset = apf.GetOutput()
264
265        self.properties.SetColor(0.89, 0.455, 0.671)  # pink7
266
267        self.pipeline = utils.OperationNode(
268            self, comment=f"#cells {self.dataset.GetNumberOfCells()}", c="#4cc9f0"
269        )

Support for UnstructuredGrid objects.

Arguments:
  • inputobj : (list, vtkUnstructuredGrid, str) A list in the form [points, cells, celltypes], or a vtkUnstructuredGrid object, or a filename

Celltypes are identified by the following convention.

actor
412    @property
413    def actor(self):
414        """Return the `vtkActor` of the object."""
415        # print("building actor")
416        gf = vtki.new("GeometryFilter")
417        gf.SetInputData(self.dataset)
418        gf.Update()
419        out = gf.GetOutput()
420        self.mapper.SetInputData(out)
421        self.mapper.Modified()
422        return self._actor

Return the vtkActor of the object.

def merge(self, *others) -> Self:
435    def merge(self, *others) -> Self:
436        """
437        Merge multiple datasets into one single `UnstrcturedGrid`.
438        """
439        apf = vtki.new("AppendFilter")
440        for o in others:
441            if isinstance(o, UnstructuredGrid):
442                apf.AddInputData(o.dataset)
443            elif isinstance(o, vtki.vtkUnstructuredGrid):
444                apf.AddInputData(o)
445            else:
446                vedo.printc("Error: cannot merge type", type(o), c="r")
447        apf.Update()
448        self._update(apf.GetOutput())
449        self.pipeline = utils.OperationNode(
450            "merge", parents=[self, *others], c="#9e2a2b"
451        )
452        return self

Merge multiple datasets into one single UnstrcturedGrid.

def copy(self, deep=True) -> UnstructuredGrid:
454    def copy(self, deep=True) -> "UnstructuredGrid":
455        """Return a copy of the object. Alias of `clone()`."""
456        return self.clone(deep=deep)

Return a copy of the object. Alias of clone().

def clone(self, deep=True) -> UnstructuredGrid:
458    def clone(self, deep=True) -> "UnstructuredGrid":
459        """Clone the UnstructuredGrid object to yield an exact copy."""
460        ug = vtki.vtkUnstructuredGrid()
461        if deep:
462            ug.DeepCopy(self.dataset)
463        else:
464            ug.ShallowCopy(self.dataset)
465        if isinstance(self, vedo.UnstructuredGrid):
466            cloned = vedo.UnstructuredGrid(ug)
467        else:
468            cloned = vedo.TetMesh(ug)
469
470        cloned.copy_properties_from(self)
471
472        cloned.pipeline = utils.OperationNode(
473            "clone", parents=[self], shape="diamond", c="#bbe1ed"
474        )
475        return cloned

Clone the UnstructuredGrid object to yield an exact copy.

def bounds(self) -> numpy.ndarray:
477    def bounds(self) -> np.ndarray:
478        """
479        Get the object bounds.
480        Returns a list in format `[xmin,xmax, ymin,ymax, zmin,zmax]`.
481        """
482        # OVERRIDE CommonAlgorithms.bounds() which is too slow
483        return np.array(self.dataset.GetBounds())

Get the object bounds. Returns a list in format [xmin,xmax, ymin,ymax, zmin,zmax].

def threshold(self, name=None, above=None, below=None, on='cells') -> Self:
485    def threshold(self, name=None, above=None, below=None, on="cells") -> Self:
486        """
487        Threshold the tetrahedral mesh by a cell scalar value.
488        Reduce to only tets which satisfy the threshold limits.
489
490        - if `above = below` will only select tets with that specific value.
491        - if `above > below` selection range is flipped.
492
493        Set keyword "on" to either "cells" or "points".
494        """
495        th = vtki.new("Threshold")
496        th.SetInputData(self.dataset)
497
498        if name is None:
499            if self.celldata.keys():
500                name = self.celldata.keys()[0]
501                th.SetInputArrayToProcess(0, 0, 0, 1, name)
502            elif self.pointdata.keys():
503                name = self.pointdata.keys()[0]
504                th.SetInputArrayToProcess(0, 0, 0, 0, name)
505            if name is None:
506                vedo.logger.warning("cannot find active array. Skip.")
507                return self
508        else:
509            if on.startswith("c"):
510                th.SetInputArrayToProcess(0, 0, 0, 1, name)
511            else:
512                th.SetInputArrayToProcess(0, 0, 0, 0, name)
513
514        if above is not None:
515            th.SetLowerThreshold(above)
516
517        if below is not None:
518            th.SetUpperThreshold(below)
519
520        th.Update()
521        return self._update(th.GetOutput())

Threshold the tetrahedral mesh by a cell scalar value. Reduce to only tets which satisfy the threshold limits.

  • if above = below will only select tets with that specific value.
  • if above > below selection range is flipped.

Set keyword "on" to either "cells" or "points".

def isosurface(self, value=None, flying_edges=False) -> vedo.mesh.Mesh:
523    def isosurface(self, value=None, flying_edges=False) -> "vedo.mesh.Mesh":
524        """
525        Return an `Mesh` isosurface extracted from the `Volume` object.
526
527        Set `value` as single float or list of values to draw the isosurface(s).
528        Use flying_edges for faster results (but sometimes can interfere with `smooth()`).
529
530        Examples:
531            - [isosurfaces1.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/isosurfaces1.py)
532
533                ![](https://vedo.embl.es/images/volumetric/isosurfaces.png)
534        """
535        scrange = self.dataset.GetScalarRange()
536
537        if flying_edges:
538            cf = vtki.new("FlyingEdges3D")
539            cf.InterpolateAttributesOn()
540        else:
541            cf = vtki.new("ContourFilter")
542            cf.UseScalarTreeOn()
543
544        cf.SetInputData(self.dataset)
545        cf.ComputeNormalsOn()
546
547        if utils.is_sequence(value):
548            cf.SetNumberOfContours(len(value))
549            for i, t in enumerate(value):
550                cf.SetValue(i, t)
551        else:
552            if value is None:
553                value = (2 * scrange[0] + scrange[1]) / 3.0
554                # print("automatic isosurface value =", value)
555            cf.SetValue(0, value)
556
557        cf.Update()
558        poly = cf.GetOutput()
559
560        out = vedo.mesh.Mesh(poly, c=None).flat()
561        out.mapper.SetScalarRange(scrange[0], scrange[1])
562
563        out.pipeline = utils.OperationNode(
564            "isosurface",
565            parents=[self],
566            comment=f"#pts {out.dataset.GetNumberOfPoints()}",
567            c="#4cc9f0:#e9c46a",
568        )
569        return out

Return an Mesh isosurface extracted from the Volume object.

Set value as single float or list of values to draw the isosurface(s). Use flying_edges for faster results (but sometimes can interfere with smooth()).

Examples:
def shrink(self, fraction=0.8) -> Self:
571    def shrink(self, fraction=0.8) -> Self:
572        """
573        Shrink the individual cells.
574
575        ![](https://vedo.embl.es/images/feats/shrink_hex.png)
576        """
577        sf = vtki.new("ShrinkFilter")
578        sf.SetInputData(self.dataset)
579        sf.SetShrinkFactor(fraction)
580        sf.Update()
581        out = sf.GetOutput()
582        self._update(out)
583        self.pipeline = utils.OperationNode(
584            "shrink", comment=f"by {fraction}", parents=[self], c="#9e2a2b"
585        )
586        return self

Shrink the individual cells.

def tomesh(self, fill=False, shrink=1.0) -> vedo.mesh.Mesh:
588    def tomesh(self, fill=False, shrink=1.0) -> "vedo.mesh.Mesh":
589        """
590        Build a polygonal `Mesh` from the current object.
591
592        If `fill=True`, the interior faces of all the cells are created.
593        (setting a `shrink` value slightly smaller than the default 1.0
594        can avoid flickering due to internal adjacent faces).
595
596        If `fill=False`, only the boundary faces will be generated.
597        """
598        gf = vtki.new("GeometryFilter")
599        if fill:
600            sf = vtki.new("ShrinkFilter")
601            sf.SetInputData(self.dataset)
602            sf.SetShrinkFactor(shrink)
603            sf.Update()
604            gf.SetInputData(sf.GetOutput())
605            gf.Update()
606            poly = gf.GetOutput()
607        else:
608            gf.SetInputData(self.dataset)
609            gf.Update()
610            poly = gf.GetOutput()
611
612        msh = vedo.mesh.Mesh(poly)
613        msh.copy_properties_from(self)
614        msh.pipeline = utils.OperationNode(
615            "tomesh", parents=[self], comment=f"fill={fill}", c="#9e2a2b:#e9c46a"
616        )
617        return msh

Build a polygonal Mesh from the current object.

If fill=True, the interior faces of all the cells are created. (setting a shrink value slightly smaller than the default 1.0 can avoid flickering due to internal adjacent faces).

If fill=False, only the boundary faces will be generated.

cell_types_array
619    @property
620    def cell_types_array(self):
621        """Return the list of cell types in the dataset."""
622        uarr = self.dataset.GetCellTypesArray()
623        return utils.vtk2numpy(uarr)

Return the list of cell types in the dataset.

def extract_cells_by_type(self, ctype) -> UnstructuredGrid:
625    def extract_cells_by_type(self, ctype) -> "UnstructuredGrid":
626        """Extract a specific cell type and return a new `UnstructuredGrid`."""
627        if isinstance(ctype, str):
628            try:
629                ctype = vtki.cell_types[ctype.upper()]
630            except KeyError:
631                vedo.logger.error(f"extract_cells_by_type: cell type {ctype} does not exist. Skip.")
632                return self
633        uarr = self.dataset.GetCellTypesArray()
634        ctarrtyp = np.where(utils.vtk2numpy(uarr) == ctype)[0]
635        uarrtyp = utils.numpy2vtk(ctarrtyp, deep=False, dtype="id")
636        selection_node = vtki.new("SelectionNode")
637        selection_node.SetFieldType(vtki.get_class("SelectionNode").CELL)
638        selection_node.SetContentType(vtki.get_class("SelectionNode").INDICES)
639        selection_node.SetSelectionList(uarrtyp)
640        selection = vtki.new("Selection")
641        selection.AddNode(selection_node)
642        es = vtki.new("ExtractSelection")
643        es.SetInputData(0, self.dataset)
644        es.SetInputData(1, selection)
645        es.Update()
646
647        ug = UnstructuredGrid(es.GetOutput())
648        ug.pipeline = utils.OperationNode(
649            "extract_cell_type", comment=f"type {ctype}", c="#edabab", parents=[self]
650        )
651        return ug

Extract a specific cell type and return a new UnstructuredGrid.

def extract_cells_by_id(self, idlist, use_point_ids=False) -> UnstructuredGrid:
653    def extract_cells_by_id(self, idlist, use_point_ids=False) -> "UnstructuredGrid":
654        """Return a new `UnstructuredGrid` composed of the specified subset of indices."""
655        selection_node = vtki.new("SelectionNode")
656        if use_point_ids:
657            selection_node.SetFieldType(vtki.get_class("SelectionNode").POINT)
658            contcells = vtki.get_class("SelectionNode").CONTAINING_CELLS()
659            selection_node.GetProperties().Set(contcells, 1)
660        else:
661            selection_node.SetFieldType(vtki.get_class("SelectionNode").CELL)
662        selection_node.SetContentType(vtki.get_class("SelectionNode").INDICES)
663        vidlist = utils.numpy2vtk(idlist, dtype="id")
664        selection_node.SetSelectionList(vidlist)
665        selection = vtki.new("Selection")
666        selection.AddNode(selection_node)
667        es = vtki.new("ExtractSelection")
668        es.SetInputData(0, self)
669        es.SetInputData(1, selection)
670        es.Update()
671
672        ug = UnstructuredGrid(es.GetOutput())
673        pr = vtki.vtkProperty()
674        pr.DeepCopy(self.properties)
675        ug.actor.SetProperty(pr)
676        ug.properties = pr
677
678        ug.mapper.SetLookupTable(utils.ctf2lut(self))
679        ug.pipeline = utils.OperationNode(
680            "extract_cells_by_id",
681            parents=[self],
682            comment=f"#cells {self.dataset.GetNumberOfCells()}",
683            c="#9e2a2b",
684        )
685        return ug

Return a new UnstructuredGrid composed of the specified subset of indices.

def find_cell(self, p: list) -> int:
687    def find_cell(self, p: list) -> int:
688        """Locate the cell that contains a point and return the cell ID."""
689        cell = vtki.vtkTetra()
690        cell_id = vtki.mutable(0)
691        tol2 = vtki.mutable(0)
692        sub_id = vtki.mutable(0)
693        pcoords = [0, 0, 0]
694        weights = [0, 0, 0