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

Support for UnstructuredGrid objects.

UnstructuredGrid(inputobj=None)
122    def __init__(self, inputobj=None):
123        """
124        Support for UnstructuredGrid objects.
125
126        Arguments:
127            inputobj : (list, vtkUnstructuredGrid, str)
128                A list in the form `[points, cells, celltypes]`,
129                or a vtkUnstructuredGrid object, or a filename
130
131        Celltypes are identified by the following 
132        [convention](https://vtk.org/doc/nightly/html/vtkCellType_8h_source.html).
133        """
134        super().__init__()
135
136        self.dataset = None
137
138        self.mapper = vtki.new("PolyDataMapper")
139        self._actor = vtki.vtkActor()
140        self._actor.retrieve_object = weak_ref_to(self)
141        self._actor.SetMapper(self.mapper)
142        self.properties = self._actor.GetProperty()
143
144        self.transform = LinearTransform()
145        self.point_locator = None
146        self.cell_locator = None
147        self.line_locator = None
148
149        self.name = "UnstructuredGrid"
150        self.filename = ""
151        self.file_size = ""
152
153        self.info = {}
154        self.time = time.time()
155        self.rendered_at = set()
156
157        ######################################
158        inputtype = str(type(inputobj))
159
160        if inputobj is None:
161            self.dataset = vtki.vtkUnstructuredGrid()
162
163        elif utils.is_sequence(inputobj):
164
165            pts, cells, celltypes = inputobj
166            assert len(cells) == len(celltypes)
167
168            self.dataset = vtki.vtkUnstructuredGrid()
169
170            if not utils.is_sequence(cells[0]):
171                tets = []
172                nf = cells[0] + 1
173                for i, cl in enumerate(cells):
174                    if i in (nf, 0):
175                        k = i + 1
176                        nf = cl + k
177                        cell = [cells[j + k] for j in range(cl)]
178                        tets.append(cell)
179                cells = tets
180
181            # This would fill the points and use those to define orientation
182            vpts = utils.numpy2vtk(pts, dtype=np.float32)
183            points = vtki.vtkPoints()
184            points.SetData(vpts)
185            self.dataset.SetPoints(points)
186
187            # Fill cells
188            # https://vtk.org/doc/nightly/html/vtkCellType_8h_source.html
189            for i, ct in enumerate(celltypes):
190                if   ct == cell_types["VERTEX"]:
191                    cell = vtki.vtkVertex()
192                elif ct == cell_types["POLY_VERTEX"]:
193                    cell = vtki.vtkPolyVertex()
194                elif ct == cell_types["TETRA"]:
195                    cell = vtki.vtkTetra()
196                elif ct == cell_types["WEDGE"]:
197                    cell = vtki.vtkWedge()
198                elif ct == cell_types["LINE"]:
199                    cell = vtki.vtkLine()
200                elif ct == cell_types["POLY_LINE"]:
201                    cell = vtki.vtkPolyLine()
202                elif ct == cell_types["TRIANGLE"]:
203                    cell = vtki.vtkTriangle()
204                elif ct == cell_types["TRIANGLE_STRIP"]:
205                    cell = vtki.vtkTriangleStrip()
206                elif ct == cell_types["POLYGON"]:
207                    cell = vtki.vtkPolygon()
208                elif ct == cell_types["PIXEL"]:
209                    cell = vtki.vtkPixel()
210                elif ct == cell_types["QUAD"]:
211                    cell = vtki.vtkQuad()
212                elif ct == cell_types["VOXEL"]:
213                    cell = vtki.vtkVoxel()
214                elif ct == cell_types["PYRAMID"]:
215                    cell = vtki.vtkPyramid()
216                elif ct == cell_types["HEXAHEDRON"]:
217                    cell = vtki.vtkHexahedron()
218                elif ct == cell_types["HEXAGONAL_PRISM"]:
219                    cell = vtki.vtkHexagonalPrism()
220                elif ct == cell_types["PENTAGONAL_PRISM"]:
221                    cell = vtki.vtkPentagonalPrism()
222                elif ct == cell_types["QUADRATIC_TETRA"]:
223                    from vtkmodules.vtkCommonDataModel import vtkQuadraticTetra
224                    cell = vtkQuadraticTetra()
225                elif ct == cell_types["QUADRATIC_HEXAHEDRON"]:
226                    from vtkmodules.vtkCommonDataModel import vtkQuadraticHexahedron
227                    cell = vtkQuadraticHexahedron()
228                elif ct == cell_types["QUADRATIC_WEDGE"]:
229                    from vtkmodules.vtkCommonDataModel import vtkQuadraticWedge
230                    cell = vtkQuadraticWedge()
231                elif ct == cell_types["QUADRATIC_PYRAMID"]:
232                    from vtkmodules.vtkCommonDataModel import vtkQuadraticPyramid
233                    cell = vtkQuadraticPyramid()
234                elif ct == cell_types["QUADRATIC_LINEAR_QUAD"]:
235                    from vtkmodules.vtkCommonDataModel import vtkQuadraticLinearQuad
236                    cell = vtkQuadraticLinearQuad()
237                elif ct == cell_types["QUADRATIC_LINEAR_WEDGE"]:
238                    from vtkmodules.vtkCommonDataModel import vtkQuadraticLinearWedge
239                    cell = vtkQuadraticLinearWedge()
240                elif ct == cell_types["BIQUADRATIC_QUADRATIC_WEDGE"]:
241                    from vtkmodules.vtkCommonDataModel import vtkBiQuadraticQuadraticWedge
242                    cell = vtkBiQuadraticQuadraticWedge()
243                elif ct == cell_types["BIQUADRATIC_QUADRATIC_HEXAHEDRON"]:
244                    from vtkmodules.vtkCommonDataModel import vtkBiQuadraticQuadraticHexahedron
245                    cell = vtkBiQuadraticQuadraticHexahedron()
246                elif ct == cell_types["BIQUADRATIC_TRIANGLE"]:
247                    from vtkmodules.vtkCommonDataModel import vtkBiQuadraticTriangle
248                    cell = vtkBiQuadraticTriangle()
249                elif ct == cell_types["CUBIC_LINE"]:
250                    from vtkmodules.vtkCommonDataModel import vtkCubicLine
251                    cell = vtkCubicLine()
252                elif ct == cell_types["CONVEX_POINT_SET"]:
253                    from vtkmodules.vtkCommonDataModel import vtkConvexPointSet
254                    cell = vtkConvexPointSet()
255                elif ct == cell_types["POLYHEDRON"]:
256                    from vtkmodules.vtkCommonDataModel import vtkPolyhedron
257                    cell = vtkPolyhedron()
258                elif ct == cell_types["HIGHER_ORDER_TRIANGLE"]:
259                    from vtkmodules.vtkCommonDataModel import vtkHigherOrderTriangle
260                    cell = vtkHigherOrderTriangle()
261                elif ct == cell_types["HIGHER_ORDER_QUAD"]:
262                    from vtkmodules.vtkCommonDataModel import vtkHigherOrderQuadrilateral
263                    cell = vtkHigherOrderQuadrilateral()
264                elif ct == cell_types["HIGHER_ORDER_TETRAHEDRON"]:
265                    from vtkmodules.vtkCommonDataModel import vtkHigherOrderTetra
266                    cell = vtkHigherOrderTetra()
267                elif ct == cell_types["HIGHER_ORDER_WEDGE"]:
268                    from vtkmodules.vtkCommonDataModel import vtkHigherOrderWedge
269                    cell = vtkHigherOrderWedge()
270                elif ct == cell_types["HIGHER_ORDER_HEXAHEDRON"]:
271                    from vtkmodules.vtkCommonDataModel import vtkHigherOrderHexahedron
272                    cell = vtkHigherOrderHexahedron()
273                elif ct == cell_types["LAGRANGE_CURVE"]:
274                    from vtkmodules.vtkCommonDataModel import vtkLagrangeCurve
275                    cell = vtkLagrangeCurve()
276                elif ct == cell_types["LAGRANGE_TRIANGLE"]:
277                    from vtkmodules.vtkCommonDataModel import vtkLagrangeTriangle
278                    cell = vtkLagrangeTriangle()
279                elif ct == cell_types["LAGRANGE_QUADRILATERAL"]:
280                    from vtkmodules.vtkCommonDataModel import vtkLagrangeQuadrilateral
281                    cell = vtkLagrangeQuadrilateral()
282                elif ct == cell_types["LAGRANGE_TETRAHEDRON"]:
283                    from vtkmodules.vtkCommonDataModel import vtkLagrangeTetra
284                    cell = vtkLagrangeTetra()
285                elif ct == cell_types["LAGRANGE_HEXAHEDRON"]:
286                    from vtkmodules.vtkCommonDataModel import vtkLagrangeHexahedron
287                    cell = vtkLagrangeHexahedron()
288                elif ct == cell_types["LAGRANGE_WEDGE"]:
289                    from vtkmodules.vtkCommonDataModel import vtkLagrangeWedge
290                    cell = vtkLagrangeWedge()
291                elif ct == cell_types["BEZIER_CURVE"]:
292                    from vtkmodules.vtkCommonDataModel import vtkBezierCurve
293                    cell = vtkBezierCurve()
294                elif ct == cell_types["BEZIER_TRIANGLE"]:
295                    from vtkmodules.vtkCommonDataModel import vtkBezierTriangle
296                    cell = vtkBezierTriangle()
297                elif ct == cell_types["BEZIER_QUADRILATERAL"]:
298                    from vtkmodules.vtkCommonDataModel import vtkBezierQuadrilateral
299                    cell = vtkBezierQuadrilateral()
300                elif ct == cell_types["BEZIER_TETRAHEDRON"]:
301                    from vtkmodules.vtkCommonDataModel import vtkBezierTetra
302                    cell = vtkBezierTetra()
303                elif ct == cell_types["BEZIER_HEXAHEDRON"]:
304                    from vtkmodules.vtkCommonDataModel import vtkBezierHexahedron
305                    cell = vtkBezierHexahedron()
306                elif ct == cell_types["BEZIER_WEDGE"]:
307                    from vtkmodules.vtkCommonDataModel import vtkBezierWedge
308                    cell = vtkBezierWedge()
309                else:
310                    vedo.logger.error(
311                        f"UnstructuredGrid: cell type {ct} not supported. Skip.")
312                    continue
313
314                cpids = cell.GetPointIds()
315                cell_conn = cells[i]
316                for j, pid in enumerate(cell_conn):
317                    cpids.SetId(j, pid)
318                self.dataset.InsertNextCell(ct, cpids)
319
320        elif "UnstructuredGrid" in inputtype:
321            self.dataset = inputobj
322
323        elif isinstance(inputobj, str):
324            if "https://" in inputobj:
325                inputobj = download(inputobj, verbose=False)
326            self.filename = inputobj
327            if inputobj.endswith(".vtu"):
328                reader = vtki.new("XMLUnstructuredGridReader")
329            else:
330                reader = vtki.new("UnstructuredGridReader")
331            self.filename = inputobj
332            reader.SetFileName(inputobj)
333            reader.Update()
334            self.dataset = reader.GetOutput()
335
336        else:
337            # this converts other types of vtk objects to UnstructuredGrid
338            apf = vtki.new("AppendFilter")
339            try:
340                apf.AddInputData(inputobj)
341            except TypeError:
342                apf.AddInputData(inputobj.dataset)
343            apf.Update()
344            self.dataset = apf.GetOutput()
345
346        self.properties.SetColor(0.89, 0.455, 0.671)  # pink7
347
348        self.pipeline = utils.OperationNode(
349            self, comment=f"#cells {self.dataset.GetNumberOfCells()}", c="#4cc9f0"
350        )

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

Return the vtkActor of the object.

def merge(self, *others):
516    def merge(self, *others):
517        """
518        Merge multiple datasets into one single `UnstrcturedGrid`.
519        """
520        apf = vtki.new("AppendFilter")
521        for o in others:
522            if isinstance(o, UnstructuredGrid):
523                apf.AddInputData(o.dataset)
524            elif isinstance(o, vtki.vtkUnstructuredGrid):
525                apf.AddInputData(o)
526            else:
527                vedo.printc("Error: cannot merge type", type(o), c="r")
528        apf.Update()
529        self._update(apf.GetOutput())
530        self.pipeline = utils.OperationNode(
531            "merge", parents=[self, *others], c="#9e2a2b"
532        )
533        return self

Merge multiple datasets into one single UnstrcturedGrid.

def copy(self, deep=True):
535    def copy(self, deep=True):
536        """Return a copy of the object. Alias of `clone()`."""
537        return self.clone(deep=deep)

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

def clone(self, deep=True):
539    def clone(self, deep=True):
540        """Clone the UnstructuredGrid object to yield an exact copy."""
541        ug = vtki.vtkUnstructuredGrid()
542        if deep:
543            ug.DeepCopy(self.dataset)
544        else:
545            ug.ShallowCopy(self.dataset)
546        if isinstance(self, vedo.UnstructuredGrid):
547            cloned = vedo.UnstructuredGrid(ug)
548        else:
549            cloned = vedo.TetMesh(ug)
550
551        cloned.copy_properties_from(self)
552
553        cloned.pipeline = utils.OperationNode(
554            "clone", parents=[self], shape="diamond", c="#bbe1ed"
555        )
556        return cloned

Clone the UnstructuredGrid object to yield an exact copy.

def bounds(self):
558    def bounds(self):
559        """
560        Get the object bounds.
561        Returns a list in format `[xmin,xmax, ymin,ymax, zmin,zmax]`.
562        """
563        # OVERRIDE CommonAlgorithms.bounds() which is too slow
564        return 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'):
566    def threshold(self, name=None, above=None, below=None, on="cells"):
567        """
568        Threshold the tetrahedral mesh by a cell scalar value.
569        Reduce to only tets which satisfy the threshold limits.
570
571        - if `above = below` will only select tets with that specific value.
572        - if `above > below` selection range is flipped.
573
574        Set keyword "on" to either "cells" or "points".
575        """
576        th = vtki.new("Threshold")
577        th.SetInputData(self.dataset)
578
579        if name is None:
580            if self.celldata.keys():
581                name = self.celldata.keys()[0]
582                th.SetInputArrayToProcess(0, 0, 0, 1, name)
583            elif self.pointdata.keys():
584                name = self.pointdata.keys()[0]
585                th.SetInputArrayToProcess(0, 0, 0, 0, name)
586            if name is None:
587                vedo.logger.warning("cannot find active array. Skip.")
588                return self
589        else:
590            if on.startswith("c"):
591                th.SetInputArrayToProcess(0, 0, 0, 1, name)
592            else:
593                th.SetInputArrayToProcess(0, 0, 0, 0, name)
594
595        if above is not None:
596            th.SetLowerThreshold(above)
597
598        if below is not None:
599            th.SetUpperThreshold(below)
600
601        th.Update()
602        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):
604    def isosurface(self, value=None, flying_edges=False):
605        """
606        Return an `Mesh` isosurface extracted from the `Volume` object.
607
608        Set `value` as single float or list of values to draw the isosurface(s).
609        Use flying_edges for faster results (but sometimes can interfere with `smooth()`).
610
611        Examples:
612            - [isosurfaces1.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/isosurfaces1.py)
613
614                ![](https://vedo.embl.es/images/volumetric/isosurfaces.png)
615        """
616        scrange = self.dataset.GetScalarRange()
617
618        if flying_edges:
619            cf = vtki.new("FlyingEdges3D")
620            cf.InterpolateAttributesOn()
621        else:
622            cf = vtki.new("ContourFilter")
623            cf.UseScalarTreeOn()
624
625        cf.SetInputData(self.dataset)
626        cf.ComputeNormalsOn()
627
628        if utils.is_sequence(value):
629            cf.SetNumberOfContours(len(value))
630            for i, t in enumerate(value):
631                cf.SetValue(i, t)
632        else:
633            if value is None:
634                value = (2 * scrange[0] + scrange[1]) / 3.0
635                # print("automatic isosurface value =", value)
636            cf.SetValue(0, value)
637
638        cf.Update()
639        poly = cf.GetOutput()
640
641        out = vedo.mesh.Mesh(poly, c=None).flat()
642        out.mapper.SetScalarRange(scrange[0], scrange[1])
643
644        out.pipeline = utils.OperationNode(
645            "isosurface",
646            parents=[self],
647            comment=f"#pts {out.dataset.GetNumberOfPoints()}",
648            c="#4cc9f0:#e9c46a",
649        )
650        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):
652    def shrink(self, fraction=0.8):
653        """
654        Shrink the individual cells.
655
656        ![](https://vedo.embl.es/images/feats/shrink_hex.png)
657        """
658        sf = vtki.new("ShrinkFilter")
659        sf.SetInputData(self.dataset)
660        sf.SetShrinkFactor(fraction)
661        sf.Update()
662        out = sf.GetOutput()
663        self._update(out)
664        self.pipeline = utils.OperationNode(
665            "shrink", comment=f"by {fraction}", parents=[self], c="#9e2a2b"
666        )
667        return self

Shrink the individual cells.

def tomesh(self, fill=False, shrink=1.0):
669    def tomesh(self, fill=False, shrink=1.0):
670        """
671        Build a polygonal `Mesh` from the current object.
672
673        If `fill=True`, the interior faces of all the cells are created.
674        (setting a `shrink` value slightly smaller than the default 1.0
675        can avoid flickering due to internal adjacent faces).
676
677        If `fill=False`, only the boundary faces will be generated.
678        """
679        gf = vtki.new("GeometryFilter")
680        if fill:
681            sf = vtki.new("ShrinkFilter")
682            sf.SetInputData(self.dataset)
683            sf.SetShrinkFactor(shrink)
684            sf.Update()
685            gf.SetInputData(sf.GetOutput())
686            gf.Update()
687            poly = gf.GetOutput()
688        else:
689            gf.SetInputData(self.dataset)
690            gf.Update()
691            poly = gf.GetOutput()
692
693        msh = vedo.mesh.Mesh(poly)
694        msh.copy_properties_from(self)
695        msh.pipeline = utils.OperationNode(
696            "tomesh", parents=[self], comment=f"fill={fill}", c="#9e2a2b:#e9c46a"
697        )
698        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

Return the list of cell types in the dataset.

def extract_cells_by_type(self, ctype):
706    def extract_cells_by_type(self, ctype):
707        """Extract a specific cell type and return a new `UnstructuredGrid`."""
708        if isinstance(ctype, str):
709            try:
710                ctype = cell_types[ctype.upper()]
711            except KeyError:
712                vedo.logger.error(f"extract_cells_by_type: cell type {ctype} does not exist. Skip.")
713                return self
714        uarr = self.dataset.GetCellTypesArray()
715        ctarrtyp = np.where(utils.vtk2numpy(uarr) == ctype)[0]
716        uarrtyp = utils.numpy2vtk(ctarrtyp, deep=False, dtype="id")
717        selection_node = vtki.new("SelectionNode")
718        selection_node.SetFieldType(vtki.get_class("SelectionNode").CELL)
719        selection_node.SetContentType(vtki.get_class("SelectionNode").INDICES)
720        selection_node.SetSelectionList(uarrtyp)
721        selection = vtki.new("Selection")
722        selection.AddNode(selection_node)
723        es = vtki.new("ExtractSelection")
724        es.SetInputData(0, self.dataset)
725        es.SetInputData(1, selection)
726        es.Update()
727
728        ug = UnstructuredGrid(es.GetOutput())
729        ug.pipeline = utils.OperationNode(
730            "extract_cell_type", comment=f"type {ctype}", c="#edabab", parents=[self]
731        )
732        return ug

Extract a specific cell type and return a new UnstructuredGrid.

def extract_cells_by_id(self, idlist, use_point_ids=False):
734    def extract_cells_by_id(self, idlist, use_point_ids=False):
735        """Return a new `UnstructuredGrid` composed of the specified subset of indices."""
736        selection_node = vtki.new("SelectionNode")
737        if use_point_ids:
738            selection_node.SetFieldType(vtki.get_class("SelectionNode").POINT)
739            contcells = vtki.get_class("SelectionNode").CONTAINING_CELLS()
740            selection_node.GetProperties().Set(contcells, 1)
741        else:
742            selection_node.SetFieldType(vtki.get_class("SelectionNode").CELL)
743        selection_node.SetContentType(vtki.get_class("SelectionNode").INDICES)
744        vidlist = utils.numpy2vtk(idlist, dtype="id")
745        selection_node.SetSelectionList(vidlist)
746        selection = vtki.new("Selection")
747        selection.AddNode(selection_node)
748        es = vtki.new("ExtractSelection")
749        es.SetInputData(0, self)
750        es.SetInputData(1, selection)
751        es.Update()
752
753        ug = UnstructuredGrid(es.GetOutput())
754        pr = vtki.vtkProperty()
755        pr.DeepCopy(self.properties)
756        ug.SetProperty(pr)
757        ug.properties = pr
758
759        ug.mapper.SetLookupTable(utils.ctf2lut(self))
760        ug.pipeline = utils.OperationNode(
761            "extract_cells_by_id",
762            parents=[self],
763            comment=f"#cells {self.dataset.GetNumberOfCells()}",
764            c="#9e2a2b",
765        )
766        return ug

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

def find_cell(self, p):
768    def find_cell(self, p):
769        """Locate the cell that contains a point and return the cell ID."""
770        cell = vtki.vtkTetra()
771        cell_id = vtki.mutable(0)
772        tol2 = vtki.mutable(0)
773        sub_id = vtki.mutable(0)
774        pcoords = [0, 0, 0]
775        weights = [0, 0, 0]
776        cid = self.dataset.FindCell(p, cell, cell_id, tol2, sub_id, pcoords, weights)
777        return cid

Locate the cell that contains a point and return the cell ID.

def clean(self):
779    def clean(self):
780        """
781        Cleanup unused points and empty cells
782        """
783        cl = vtki.new("StaticCleanUnstructuredGrid")
784        cl.SetInputData(self.dataset)
785        cl.RemoveUnusedPointsOn()
786        cl.ProduceMergeMapOff()
787        cl.AveragePointDataOff()
788        cl.Update()