vedo.tetmesh

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