vedo.mesh

Submodule to work with polygonal meshes

   1#!/usr/bin/env python3
   2# -*- coding: utf-8 -*-
   3import numpy as np
   4from typing import List, Tuple, Union, MutableSequence, Any
   5from typing_extensions import Self
   6
   7import vedo.vtkclasses as vtki  # a wrapper for lazy imports
   8
   9import vedo
  10from vedo.colors import get_color
  11from vedo.pointcloud import Points
  12from vedo.utils import buildPolyData, is_sequence, mag, precision
  13from vedo.utils import numpy2vtk, vtk2numpy, OperationNode
  14from vedo.visual import MeshVisual
  15
  16__docformat__ = "google"
  17
  18__doc__ = """
  19Submodule to work with polygonal meshes
  20
  21![](https://vedo.embl.es/images/advanced/mesh_smoother2.png)
  22"""
  23
  24__all__ = ["Mesh"]
  25
  26
  27####################################################
  28class Mesh(MeshVisual, Points):
  29    """
  30    Build an instance of object `Mesh` derived from `vedo.PointCloud`.
  31    """
  32
  33    def __init__(self, inputobj=None, c="gold", alpha=1):
  34        """
  35        Initialize a ``Mesh`` object.
  36
  37        Arguments:
  38            inputobj : (str, vtkPolyData, vtkActor, vedo.Mesh)
  39                If inputobj is `None` an empty mesh is created.
  40                If inputobj is a `str` then it is interpreted as the name of a file to load as mesh.
  41                If inputobj is an `vtkPolyData` or `vtkActor` or `vedo.Mesh`
  42                then a shallow copy of it is created.
  43                If inputobj is a `vedo.Mesh` then a shallow copy of it is created.
  44
  45        Examples:
  46            - [buildmesh.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/buildmesh.py)
  47            (and many others!)
  48
  49            ![](https://vedo.embl.es/images/basic/buildmesh.png)
  50        """
  51        # print("INIT MESH", super())
  52        super().__init__()
  53
  54        self.name = "Mesh"
  55
  56        if inputobj is None:
  57            # self.dataset = vtki.vtkPolyData()
  58            pass
  59
  60        elif isinstance(inputobj, str):
  61            self.dataset = vedo.file_io.load(inputobj).dataset
  62            self.filename = inputobj
  63
  64        elif isinstance(inputobj, vtki.vtkPolyData):
  65            # self.dataset.DeepCopy(inputobj) # NO
  66            self.dataset = inputobj
  67            if self.dataset.GetNumberOfCells() == 0:
  68                carr = vtki.vtkCellArray()
  69                for i in range(inputobj.GetNumberOfPoints()):
  70                    carr.InsertNextCell(1)
  71                    carr.InsertCellPoint(i)
  72                self.dataset.SetVerts(carr)
  73
  74        elif isinstance(inputobj, Mesh):
  75            self.dataset = inputobj.dataset
  76
  77        elif is_sequence(inputobj):
  78            ninp = len(inputobj)
  79            if   ninp == 4:  # assume input is [vertices, faces, lines, strips]
  80                self.dataset = buildPolyData(inputobj[0], inputobj[1], inputobj[2], inputobj[3])
  81            elif ninp == 3:  # assume input is [vertices, faces, lines]
  82                self.dataset = buildPolyData(inputobj[0], inputobj[1], inputobj[2])
  83            elif ninp == 2:  # assume input is [vertices, faces]
  84                self.dataset = buildPolyData(inputobj[0], inputobj[1])
  85            elif ninp == 1:  # assume input is [vertices]
  86                self.dataset = buildPolyData(inputobj[0])
  87            else:
  88                vedo.logger.error("input must be a list of max 4 elements.")
  89                raise ValueError()
  90
  91        elif isinstance(inputobj, vtki.vtkActor):
  92            self.dataset.DeepCopy(inputobj.GetMapper().GetInput())
  93            v = inputobj.GetMapper().GetScalarVisibility()
  94            self.mapper.SetScalarVisibility(v)
  95            pr = vtki.vtkProperty()
  96            pr.DeepCopy(inputobj.GetProperty())
  97            self.actor.SetProperty(pr)
  98            self.properties = pr
  99
 100        elif isinstance(inputobj, (vtki.vtkStructuredGrid, vtki.vtkRectilinearGrid)):
 101            gf = vtki.new("GeometryFilter")
 102            gf.SetInputData(inputobj)
 103            gf.Update()
 104            self.dataset = gf.GetOutput()
 105
 106        elif "meshlab" in str(type(inputobj)):
 107            self.dataset = vedo.utils.meshlab2vedo(inputobj).dataset
 108
 109        elif "meshlib" in str(type(inputobj)):
 110            import meshlib.mrmeshnumpy as mrmeshnumpy
 111            self.dataset = buildPolyData(
 112                mrmeshnumpy.getNumpyVerts(inputobj),
 113                mrmeshnumpy.getNumpyFaces(inputobj.topology),
 114            )
 115
 116        elif "trimesh" in str(type(inputobj)):
 117            self.dataset = vedo.utils.trimesh2vedo(inputobj).dataset
 118
 119        elif "meshio" in str(type(inputobj)):
 120            # self.dataset = vedo.utils.meshio2vedo(inputobj) ##TODO
 121            if len(inputobj.cells) > 0:
 122                mcells = []
 123                for cellblock in inputobj.cells:
 124                    if cellblock.type in ("triangle", "quad"):
 125                        mcells += cellblock.data.tolist()
 126                self.dataset = buildPolyData(inputobj.points, mcells)
 127            else:
 128                self.dataset = buildPolyData(inputobj.points, None)
 129            # add arrays:
 130            try:
 131                if len(inputobj.point_data) > 0:
 132                    for k in inputobj.point_data.keys():
 133                        vdata = numpy2vtk(inputobj.point_data[k])
 134                        vdata.SetName(str(k))
 135                        self.dataset.GetPointData().AddArray(vdata)
 136            except AssertionError:
 137                print("Could not add meshio point data, skip.")
 138
 139        else:
 140            try:
 141                gf = vtki.new("GeometryFilter")
 142                gf.SetInputData(inputobj)
 143                gf.Update()
 144                self.dataset = gf.GetOutput()
 145            except:
 146                vedo.logger.error(f"cannot build mesh from type {type(inputobj)}")
 147                raise RuntimeError()
 148
 149        self.mapper.SetInputData(self.dataset)
 150        self.actor.SetMapper(self.mapper)
 151
 152        self.properties.SetInterpolationToPhong()
 153        self.properties.SetColor(get_color(c))
 154
 155        if alpha is not None:
 156            self.properties.SetOpacity(alpha)
 157
 158        self.mapper.SetInterpolateScalarsBeforeMapping(
 159            vedo.settings.interpolate_scalars_before_mapping
 160        )
 161
 162        if vedo.settings.use_polygon_offset:
 163            self.mapper.SetResolveCoincidentTopologyToPolygonOffset()
 164            pof = vedo.settings.polygon_offset_factor
 165            pou = vedo.settings.polygon_offset_units
 166            self.mapper.SetResolveCoincidentTopologyPolygonOffsetParameters(pof, pou)
 167
 168        n = self.dataset.GetNumberOfPoints()
 169        self.pipeline = OperationNode(self, comment=f"#pts {n}")
 170
 171    def _repr_html_(self):
 172        """
 173        HTML representation of the Mesh object for Jupyter Notebooks.
 174
 175        Returns:
 176            HTML text with the image and some properties.
 177        """
 178        import io
 179        import base64
 180        from PIL import Image
 181
 182        library_name = "vedo.mesh.Mesh"
 183        help_url = "https://vedo.embl.es/docs/vedo/mesh.html#Mesh"
 184
 185        arr = self.thumbnail()
 186        im = Image.fromarray(arr)
 187        buffered = io.BytesIO()
 188        im.save(buffered, format="PNG", quality=100)
 189        encoded = base64.b64encode(buffered.getvalue()).decode("utf-8")
 190        url = "data:image/png;base64," + encoded
 191        image = f"<img src='{url}'></img>"
 192
 193        bounds = "<br/>".join(
 194            [
 195                precision(min_x, 4) + " ... " + precision(max_x, 4)
 196                for min_x, max_x in zip(self.bounds()[::2], self.bounds()[1::2])
 197            ]
 198        )
 199        average_size = "{size:.3f}".format(size=self.average_size())
 200
 201        help_text = ""
 202        if self.name:
 203            help_text += f"<b> {self.name}: &nbsp&nbsp</b>"
 204        help_text += '<b><a href="' + help_url + '" target="_blank">' + library_name + "</a></b>"
 205        if self.filename:
 206            dots = ""
 207            if len(self.filename) > 30:
 208                dots = "..."
 209            help_text += f"<br/><code><i>({dots}{self.filename[-30:]})</i></code>"
 210
 211        pdata = ""
 212        if self.dataset.GetPointData().GetScalars():
 213            if self.dataset.GetPointData().GetScalars().GetName():
 214                name = self.dataset.GetPointData().GetScalars().GetName()
 215                pdata = "<tr><td><b> point data array </b></td><td>" + name + "</td></tr>"
 216
 217        cdata = ""
 218        if self.dataset.GetCellData().GetScalars():
 219            if self.dataset.GetCellData().GetScalars().GetName():
 220                name = self.dataset.GetCellData().GetScalars().GetName()
 221                cdata = "<tr><td><b> cell data array </b></td><td>" + name + "</td></tr>"
 222
 223        allt = [
 224            "<table>",
 225            "<tr>",
 226            "<td>",
 227            image,
 228            "</td>",
 229            "<td style='text-align: center; vertical-align: center;'><br/>",
 230            help_text,
 231            "<table>",
 232            "<tr><td><b> bounds </b> <br/> (x/y/z) </td><td>" + str(bounds) + "</td></tr>",
 233            "<tr><td><b> center of mass </b></td><td>"
 234            + precision(self.center_of_mass(), 3)
 235            + "</td></tr>",
 236            "<tr><td><b> average size </b></td><td>" + str(average_size) + "</td></tr>",
 237            "<tr><td><b> nr. points&nbsp/&nbspfaces </b></td><td>"
 238            + str(self.npoints)
 239            + "&nbsp/&nbsp"
 240            + str(self.ncells)
 241            + "</td></tr>",
 242            pdata,
 243            cdata,
 244            "</table>",
 245            "</table>",
 246        ]
 247        return "\n".join(allt)
 248
 249    def faces(self, ids=()):
 250        """DEPRECATED. Use property `mesh.cells` instead."""
 251        vedo.printc("WARNING: use property mesh.cells instead of mesh.faces()",c='y')
 252        return self.cells
 253    
 254    @property
 255    def edges(self):
 256        """Return an array containing the edges connectivity."""
 257        extractEdges = vtki.new("ExtractEdges")
 258        extractEdges.SetInputData(self.dataset)
 259        # eed.UseAllPointsOn()
 260        extractEdges.Update()
 261        lpoly = extractEdges.GetOutput()
 262
 263        arr1d = vtk2numpy(lpoly.GetLines().GetData())
 264        # [nids1, id0 ... idn, niids2, id0 ... idm,  etc].
 265
 266        i = 0
 267        conn = []
 268        n = len(arr1d)
 269        for _ in range(n):
 270            cell = [arr1d[i + k + 1] for k in range(arr1d[i])]
 271            conn.append(cell)
 272            i += arr1d[i] + 1
 273            if i >= n:
 274                break
 275        return conn  # cannot always make a numpy array of it!
 276
 277    @property
 278    def cell_normals(self):
 279        """
 280        Retrieve face normals as a numpy array.
 281        Check out also `compute_normals(cells=True)` and `compute_normals_with_pca()`.
 282        """
 283        vtknormals = self.dataset.GetCellData().GetNormals()
 284        numpy_normals = vtk2numpy(vtknormals)
 285        if len(numpy_normals) == 0 and len(self.cells) != 0:
 286            raise ValueError("VTK failed to return any normal vectors. You may need to call `Mesh.compute_normals()` before accessing `Mesh.cell_normals`.")
 287        return numpy_normals
 288
 289    def compute_normals(self, points=True, cells=True, feature_angle=None, consistency=True) -> Self:
 290        """
 291        Compute cell and vertex normals for the mesh.
 292
 293        Arguments:
 294            points : (bool)
 295                do the computation for the vertices too
 296            cells : (bool)
 297                do the computation for the cells too
 298            feature_angle : (float)
 299                specify the angle that defines a sharp edge.
 300                If the difference in angle across neighboring polygons is greater than this value,
 301                the shared edge is considered "sharp" and it is split.
 302            consistency : (bool)
 303                turn on/off the enforcement of consistent polygon ordering.
 304
 305        .. warning::
 306            If `feature_angle` is set then the Mesh can be modified, and it
 307            can have a different nr. of vertices from the original.
 308
 309            Note that the appearance of the mesh may change if the normals are computed,
 310            as shading is automatically enabled when such information is present.
 311            Use `mesh.flat()` to avoid smoothing effects.
 312        """
 313        pdnorm = vtki.new("PolyDataNormals")
 314        pdnorm.SetInputData(self.dataset)
 315        pdnorm.SetComputePointNormals(points)
 316        pdnorm.SetComputeCellNormals(cells)
 317        pdnorm.SetConsistency(consistency)
 318        pdnorm.FlipNormalsOff()
 319        if feature_angle:
 320            pdnorm.SetSplitting(True)
 321            pdnorm.SetFeatureAngle(feature_angle)
 322        else:
 323            pdnorm.SetSplitting(False)
 324        pdnorm.Update()
 325        out = pdnorm.GetOutput()
 326        self._update(out, reset_locators=False)
 327        return self
 328
 329    def reverse(self, cells=True, normals=False) -> Self:
 330        """
 331        Reverse the order of polygonal cells
 332        and/or reverse the direction of point and cell normals.
 333
 334        Two flags are used to control these operations:
 335            - `cells=True` reverses the order of the indices in the cell connectivity list.
 336                If cell is a list of IDs only those cells will be reversed.
 337            - `normals=True` reverses the normals by multiplying the normal vector by -1
 338                (both point and cell normals, if present).
 339        """
 340        poly = self.dataset
 341
 342        if is_sequence(cells):
 343            for cell in cells:
 344                poly.ReverseCell(cell)
 345            poly.GetCellData().Modified()
 346            return self  ##############
 347
 348        rev = vtki.new("ReverseSense")
 349        if cells:
 350            rev.ReverseCellsOn()
 351        else:
 352            rev.ReverseCellsOff()
 353        if normals:
 354            rev.ReverseNormalsOn()
 355        else:
 356            rev.ReverseNormalsOff()
 357        rev.SetInputData(poly)
 358        rev.Update()
 359        self._update(rev.GetOutput(), reset_locators=False)
 360        self.pipeline = OperationNode("reverse", parents=[self])
 361        return self
 362
 363    def volume(self) -> float:
 364        """
 365        Compute the volume occupied by mesh.
 366        The mesh must be triangular for this to work.
 367        To triangulate a mesh use `mesh.triangulate()`.
 368        """
 369        mass = vtki.new("MassProperties")
 370        mass.SetGlobalWarningDisplay(0)
 371        mass.SetInputData(self.dataset)
 372        mass.Update()
 373        mass.SetGlobalWarningDisplay(1)
 374        return mass.GetVolume()
 375
 376    def area(self) -> float:
 377        """
 378        Compute the surface area of the mesh.
 379        The mesh must be triangular for this to work.
 380        To triangulate a mesh use `mesh.triangulate()`.
 381        """
 382        mass = vtki.new("MassProperties")
 383        mass.SetGlobalWarningDisplay(0)
 384        mass.SetInputData(self.dataset)
 385        mass.Update()
 386        mass.SetGlobalWarningDisplay(1)
 387        return mass.GetSurfaceArea()
 388
 389    def is_closed(self) -> bool:
 390        """
 391        Return `True` if the mesh is watertight.
 392        Note that if the mesh contains coincident points the result may be flase.
 393        Use in this case `mesh.clean()` to merge coincident points.
 394        """
 395        fe = vtki.new("FeatureEdges")
 396        fe.BoundaryEdgesOn()
 397        fe.FeatureEdgesOff()
 398        fe.NonManifoldEdgesOn()
 399        fe.SetInputData(self.dataset)
 400        fe.Update()
 401        ne = fe.GetOutput().GetNumberOfCells()
 402        return not bool(ne)
 403
 404    def is_manifold(self) -> bool:
 405        """Return `True` if the mesh is manifold."""
 406        fe = vtki.new("FeatureEdges")
 407        fe.BoundaryEdgesOff()
 408        fe.FeatureEdgesOff()
 409        fe.NonManifoldEdgesOn()
 410        fe.SetInputData(self.dataset)
 411        fe.Update()
 412        ne = fe.GetOutput().GetNumberOfCells()
 413        return not bool(ne)
 414
 415    def non_manifold_faces(self, remove=True, tol="auto") -> Self:
 416        """
 417        Detect and (try to) remove non-manifold faces of a triangular mesh:
 418
 419            - set `remove` to `False` to mark cells without removing them.
 420            - set `tol=0` for zero-tolerance, the result will be manifold but with holes.
 421            - set `tol>0` to cut off non-manifold faces, and try to recover the good ones.
 422            - set `tol="auto"` to make an automatic choice of the tolerance.
 423        """
 424        # mark original point and cell ids
 425        self.add_ids()
 426        toremove = self.boundaries(
 427            boundary_edges=False,
 428            non_manifold_edges=True,
 429            cell_edge=True,
 430            return_cell_ids=True,
 431        )
 432        if len(toremove) == 0: # type: ignore
 433            return self
 434
 435        points = self.vertices
 436        faces = self.cells
 437        centers = self.cell_centers
 438
 439        copy = self.clone()
 440        copy.delete_cells(toremove).clean()
 441        copy.compute_normals(cells=False)
 442        normals = copy.vertex_normals
 443        deltas, deltas_i = [], []
 444
 445        for i in vedo.utils.progressbar(toremove, delay=3, title="recover faces"):
 446            pids = copy.closest_point(centers[i], n=3, return_point_id=True)
 447            norms = normals[pids]
 448            n = np.mean(norms, axis=0)
 449            dn = np.linalg.norm(n)
 450            if not dn:
 451                continue
 452            n = n / dn
 453
 454            p0, p1, p2 = points[faces[i]][:3]
 455            v = np.cross(p1 - p0, p2 - p0)
 456            lv = np.linalg.norm(v)
 457            if not lv:
 458                continue
 459            v = v / lv
 460
 461            cosa = 1 - np.dot(n, v)
 462            deltas.append(cosa)
 463            deltas_i.append(i)
 464
 465        recover = []
 466        if len(deltas) > 0:
 467            mean_delta = np.mean(deltas)
 468            err_delta = np.std(deltas)
 469            txt = ""
 470            if tol == "auto":  # automatic choice
 471                tol = mean_delta / 5
 472                txt = f"\n Automatic tol. : {tol: .4f}"
 473            for i, cosa in zip(deltas_i, deltas):
 474                if cosa < tol:
 475                    recover.append(i)
 476
 477            vedo.logger.info(
 478                f"\n --------- Non manifold faces ---------"
 479                f"\n Average tol.   : {mean_delta: .4f} +- {err_delta: .4f}{txt}"
 480                f"\n Removed faces  : {len(toremove)}" # type: ignore
 481                f"\n Recovered faces: {len(recover)}"
 482            )
 483
 484        toremove = list(set(toremove) - set(recover)) # type: ignore
 485
 486        if not remove:
 487            mark = np.zeros(self.ncells, dtype=np.uint8)
 488            mark[recover] = 1
 489            mark[toremove] = 2
 490            self.celldata["NonManifoldCell"] = mark
 491        else:
 492            self.delete_cells(toremove) # type: ignore
 493
 494        self.pipeline = OperationNode(
 495            "non_manifold_faces",
 496            parents=[self],
 497            comment=f"#cells {self.dataset.GetNumberOfCells()}",
 498        )
 499        return self
 500
 501
 502    def euler_characteristic(self) -> int:
 503        """
 504        Compute the Euler characteristic of the mesh.
 505        The Euler characteristic is a topological invariant for surfaces.
 506        """
 507        return self.npoints - len(self.edges) + self.ncells
 508
 509    def genus(self) -> int:
 510        """
 511        Compute the genus of the mesh.
 512        The genus is a topological invariant for surfaces.
 513        """
 514        nb = len(self.boundaries().split()) - 1
 515        return (2 - self.euler_characteristic() - nb ) / 2
 516    
 517    def to_reeb_graph(self, field_id=0):
 518        """
 519        Convert the mesh into a Reeb graph.
 520        The Reeb graph is a topological structure that captures the evolution
 521        of the level sets of a scalar field.
 522
 523        Arguments:
 524            field_id : (int)
 525                the id of the scalar field to use.
 526        
 527        Example:
 528            ```python
 529            from vedo import *
 530            mesh = Mesh("https://discourse.paraview.org/uploads/short-url/qVuZ1fiRjwhE1qYtgGE2HGXybgo.stl")
 531            mesh.rotate_x(10).rotate_y(15).alpha(0.5)
 532            mesh.pointdata["scalars"] = mesh.vertices[:, 2]
 533
 534            printc("is_closed  :", mesh.is_closed())
 535            printc("is_manifold:", mesh.is_manifold())
 536            printc("euler_char :", mesh.euler_characteristic())
 537            printc("genus      :", mesh.genus())
 538
 539            reeb = mesh.to_reeb_graph()
 540            ids = reeb[0].pointdata["Vertex Ids"]
 541            pts = Points(mesh.vertices[ids], r=10)
 542
 543            show([[mesh, pts], reeb], N=2, sharecam=False)
 544            ```
 545        """
 546        rg = vtki.new("PolyDataToReebGraphFilter")
 547        rg.SetInputData(self.dataset)
 548        rg.SetFieldId(field_id)
 549        rg.Update()
 550        gr = vedo.pyplot.DirectedGraph()
 551        gr.mdg = rg.GetOutput()
 552        gr.build()
 553        return gr
 554
 555
 556    def shrink(self, fraction=0.85) -> Self:
 557        """
 558        Shrink the triangle polydata in the representation of the input mesh.
 559
 560        Examples:
 561            - [shrink.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/shrink.py)
 562
 563            ![](https://vedo.embl.es/images/basic/shrink.png)
 564        """
 565        # Overriding base class method core.shrink()
 566        shrink = vtki.new("ShrinkPolyData")
 567        shrink.SetInputData(self.dataset)
 568        shrink.SetShrinkFactor(fraction)
 569        shrink.Update()
 570        self._update(shrink.GetOutput())
 571        self.pipeline = OperationNode("shrink", parents=[self])
 572        return self
 573
 574    def cap(self, return_cap=False) -> Self:
 575        """
 576        Generate a "cap" on a clipped mesh, or caps sharp edges.
 577
 578        Examples:
 579            - [cut_and_cap.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/cut_and_cap.py)
 580
 581            ![](https://vedo.embl.es/images/advanced/cutAndCap.png)
 582
 583        See also: `join()`, `join_segments()`, `slice()`.
 584        """
 585        fe = vtki.new("FeatureEdges")
 586        fe.SetInputData(self.dataset)
 587        fe.BoundaryEdgesOn()
 588        fe.FeatureEdgesOff()
 589        fe.NonManifoldEdgesOff()
 590        fe.ManifoldEdgesOff()
 591        fe.Update()
 592
 593        stripper = vtki.new("Stripper")
 594        stripper.SetInputData(fe.GetOutput())
 595        stripper.JoinContiguousSegmentsOn()
 596        stripper.Update()
 597
 598        boundary_poly = vtki.vtkPolyData()
 599        boundary_poly.SetPoints(stripper.GetOutput().GetPoints())
 600        boundary_poly.SetPolys(stripper.GetOutput().GetLines())
 601
 602        rev = vtki.new("ReverseSense")
 603        rev.ReverseCellsOn()
 604        rev.SetInputData(boundary_poly)
 605        rev.Update()
 606
 607        tf = vtki.new("TriangleFilter")
 608        tf.SetInputData(rev.GetOutput())
 609        tf.Update()
 610
 611        if return_cap:
 612            m = Mesh(tf.GetOutput())
 613            m.pipeline = OperationNode(
 614                "cap", parents=[self], comment=f"#pts {m.dataset.GetNumberOfPoints()}"
 615            )
 616            m.name = "MeshCap"
 617            return m
 618
 619        polyapp = vtki.new("AppendPolyData")
 620        polyapp.AddInputData(self.dataset)
 621        polyapp.AddInputData(tf.GetOutput())
 622        polyapp.Update()
 623
 624        self._update(polyapp.GetOutput())
 625        self.clean()
 626
 627        self.pipeline = OperationNode(
 628            "capped", parents=[self], comment=f"#pts {self.dataset.GetNumberOfPoints()}"
 629        )
 630        return self
 631
 632    def join(self, polys=True, reset=False) -> Self:
 633        """
 634        Generate triangle strips and/or polylines from
 635        input polygons, triangle strips, and lines.
 636
 637        Input polygons are assembled into triangle strips only if they are triangles;
 638        other types of polygons are passed through to the output and not stripped.
 639        Use mesh.triangulate() to triangulate non-triangular polygons prior to running
 640        this filter if you need to strip all the data.
 641
 642        Also note that if triangle strips or polylines are present in the input
 643        they are passed through and not joined nor extended.
 644        If you wish to strip these use mesh.triangulate() to fragment the input
 645        into triangles and lines prior to applying join().
 646
 647        Arguments:
 648            polys : (bool)
 649                polygonal segments will be joined if they are contiguous
 650            reset : (bool)
 651                reset points ordering
 652
 653        Warning:
 654            If triangle strips or polylines exist in the input data
 655            they will be passed through to the output data.
 656            This filter will only construct triangle strips if triangle polygons
 657            are available; and will only construct polylines if lines are available.
 658
 659        Example:
 660            ```python
 661            from vedo import *
 662            c1 = Cylinder(pos=(0,0,0), r=2, height=3, axis=(1,.0,0), alpha=.1).triangulate()
 663            c2 = Cylinder(pos=(0,0,2), r=1, height=2, axis=(0,.3,1), alpha=.1).triangulate()
 664            intersect = c1.intersect_with(c2).join(reset=True)
 665            spline = Spline(intersect).c('blue').lw(5)
 666            show(c1, c2, spline, intersect.labels('id'), axes=1).close()
 667            ```
 668            ![](https://vedo.embl.es/images/feats/line_join.png)
 669        """
 670        sf = vtki.new("Stripper")
 671        sf.SetPassThroughCellIds(True)
 672        sf.SetPassThroughPointIds(True)
 673        sf.SetJoinContiguousSegments(polys)
 674        sf.SetInputData(self.dataset)
 675        sf.Update()
 676        if reset:
 677            poly = sf.GetOutput()
 678            cpd = vtki.new("CleanPolyData")
 679            cpd.PointMergingOn()
 680            cpd.ConvertLinesToPointsOn()
 681            cpd.ConvertPolysToLinesOn()
 682            cpd.ConvertStripsToPolysOn()
 683            cpd.SetInputData(poly)
 684            cpd.Update()
 685            poly = cpd.GetOutput()
 686            vpts = poly.GetCell(0).GetPoints().GetData()
 687            poly.GetPoints().SetData(vpts)
 688        else:
 689            poly = sf.GetOutput()
 690
 691        self._update(poly)
 692
 693        self.pipeline = OperationNode(
 694            "join", parents=[self], comment=f"#pts {self.dataset.GetNumberOfPoints()}"
 695        )
 696        return self
 697
 698    def join_segments(self, closed=True, tol=1e-03) -> list:
 699        """
 700        Join line segments into contiguous lines.
 701        Useful to call with `triangulate()` method.
 702
 703        Returns:
 704            list of `shapes.Lines`
 705
 706        Example:
 707            ```python
 708            from vedo import *
 709            msh = Torus().alpha(0.1).wireframe()
 710            intersection = msh.intersect_with_plane(normal=[1,1,1]).c('purple5')
 711            slices = [s.triangulate() for s in intersection.join_segments()]
 712            show(msh, intersection, merge(slices), axes=1, viewup='z')
 713            ```
 714            ![](https://vedo.embl.es/images/feats/join_segments.jpg)
 715        """
 716        vlines = []
 717        for ipiece, outline in enumerate(self.split(must_share_edge=False)): # type: ignore
 718
 719            outline.clean()
 720            pts = outline.vertices
 721            if len(pts) < 3:
 722                continue
 723            avesize = outline.average_size()
 724            lines = outline.lines
 725            # print("---lines", lines, "in piece", ipiece)
 726            tol = avesize / pts.shape[0] * tol
 727
 728            k = 0
 729            joinedpts = [pts[k]]
 730            for _ in range(len(pts)):
 731                pk = pts[k]
 732                for j, line in enumerate(lines):
 733
 734                    id0, id1 = line[0], line[-1]
 735                    p0, p1 = pts[id0], pts[id1]
 736
 737                    if np.linalg.norm(p0 - pk) < tol:
 738                        n = len(line)
 739                        for m in range(1, n):
 740                            joinedpts.append(pts[line[m]])
 741                        # joinedpts.append(p1)
 742                        k = id1
 743                        lines.pop(j)
 744                        break
 745
 746                    elif np.linalg.norm(p1 - pk) < tol:
 747                        n = len(line)
 748                        for m in reversed(range(0, n - 1)):
 749                            joinedpts.append(pts[line[m]])
 750                        # joinedpts.append(p0)
 751                        k = id0
 752                        lines.pop(j)
 753                        break
 754
 755            if len(joinedpts) > 1:
 756                newline = vedo.shapes.Line(joinedpts, closed=closed)
 757                newline.clean()
 758                newline.actor.SetProperty(self.properties)
 759                newline.properties = self.properties
 760                newline.pipeline = OperationNode(
 761                    "join_segments",
 762                    parents=[self],
 763                    comment=f"#pts {newline.dataset.GetNumberOfPoints()}",
 764                )
 765                vlines.append(newline)
 766
 767        return vlines
 768
 769    def join_with_strips(self, b1, closed=True) -> Self:
 770        """
 771        Join booundary lines by creating a triangle strip between them.
 772
 773        Example:
 774        ```python
 775        from vedo import *
 776        m1 = Cylinder(cap=False).boundaries()
 777        m2 = Cylinder(cap=False).boundaries().pos(0.2,0,1)
 778        strips = m1.join_with_strips(m2)
 779        show(m1, m2, strips, axes=1).close()
 780        ```
 781        """
 782        b0 = self.clone().join()
 783        b1 = b1.clone().join()
 784
 785        vertices0 = b0.vertices.tolist()
 786        vertices1 = b1.vertices.tolist()
 787
 788        lines0 = b0.lines
 789        lines1 = b1.lines
 790        m =  len(lines0)
 791        assert m == len(lines1), (
 792            "lines must have the same number of points\n"
 793            f"line has {m} points in b0 and {len(lines1)} in b1"
 794        )
 795
 796        strips = []
 797        points: List[Any] = []
 798
 799        for j in range(m):
 800
 801            ids0j = list(lines0[j])
 802            ids1j = list(lines1[j])
 803
 804            n = len(ids0j)
 805            assert n == len(ids1j), (
 806                "lines must have the same number of points\n"
 807                f"line {j} has {n} points in b0 and {len(ids1j)} in b1"
 808            )
 809
 810            if closed:
 811                ids0j.append(ids0j[0])
 812                ids1j.append(ids1j[0])
 813                vertices0.append(vertices0[ids0j[0]])
 814                vertices1.append(vertices1[ids1j[0]])
 815                n = n + 1
 816
 817            strip = []  # create a triangle strip
 818            npt = len(points)
 819            for ipt in range(n):
 820                points.append(vertices0[ids0j[ipt]])
 821                points.append(vertices1[ids1j[ipt]])
 822
 823            strip = list(range(npt, npt + 2*n))
 824            strips.append(strip)
 825
 826        return Mesh([points, [], [], strips], c="k6")
 827
 828    def split_polylines(self) -> Self:
 829        """Split polylines into separate segments."""
 830        tf = vtki.new("TriangleFilter")
 831        tf.SetPassLines(True)
 832        tf.SetPassVerts(False)
 833        tf.SetInputData(self.dataset)
 834        tf.Update()
 835        self._update(tf.GetOutput(), reset_locators=False)
 836        self.lw(0).lighting("default").pickable()
 837        self.pipeline = OperationNode(
 838            "split_polylines", parents=[self], 
 839            comment=f"#lines {self.dataset.GetNumberOfLines()}"
 840        )
 841        return self
 842
 843    def slice(self, origin=(0, 0, 0), normal=(1, 0, 0)) -> Self:
 844        """
 845        Slice a mesh with a plane and fill the contour.
 846
 847        Example:
 848            ```python
 849            from vedo import *
 850            msh = Mesh(dataurl+"bunny.obj").alpha(0.1).wireframe()
 851            mslice = msh.slice(normal=[0,1,0.3], origin=[0,0.16,0])
 852            mslice.c('purple5')
 853            show(msh, mslice, axes=1)
 854            ```
 855            ![](https://vedo.embl.es/images/feats/mesh_slice.jpg)
 856
 857        See also: `join()`, `join_segments()`, `cap()`, `cut_with_plane()`.
 858        """
 859        intersection = self.intersect_with_plane(origin=origin, normal=normal)
 860        slices = [s.triangulate() for s in intersection.join_segments()]
 861        mslices = vedo.pointcloud.merge(slices)
 862        if mslices:
 863            mslices.name = "MeshSlice"
 864            mslices.pipeline = OperationNode("slice", parents=[self], comment=f"normal = {normal}")
 865        return mslices
 866
 867    def triangulate(self, verts=True, lines=True) -> Self:
 868        """
 869        Converts mesh polygons into triangles.
 870
 871        If the input mesh is only made of 2D lines (no faces) the output will be a triangulation
 872        that fills the internal area. The contours may be concave, and may even contain holes,
 873        i.e. a contour may contain an internal contour winding in the opposite
 874        direction to indicate that it is a hole.
 875
 876        Arguments:
 877            verts : (bool)
 878                if True, break input vertex cells into individual vertex cells (one point per cell).
 879                If False, the input vertex cells will be ignored.
 880            lines : (bool)
 881                if True, break input polylines into line segments.
 882                If False, input lines will be ignored and the output will have no lines.
 883        """
 884        if self.dataset.GetNumberOfPolys() or self.dataset.GetNumberOfStrips():
 885            # print("Using vtkTriangleFilter")
 886            tf = vtki.new("TriangleFilter")
 887            tf.SetPassLines(lines)
 888            tf.SetPassVerts(verts)
 889
 890        elif self.dataset.GetNumberOfLines():
 891            # print("Using vtkContourTriangulator")
 892            tf = vtki.new("ContourTriangulator")
 893            tf.TriangulationErrorDisplayOn()
 894
 895        else:
 896            vedo.logger.debug("input in triangulate() seems to be void! Skip.")
 897            return self
 898
 899        tf.SetInputData(self.dataset)
 900        tf.Update()
 901        self._update(tf.GetOutput(), reset_locators=False)
 902        self.lw(0).lighting("default").pickable()
 903
 904        self.pipeline = OperationNode(
 905            "triangulate", parents=[self], comment=f"#cells {self.dataset.GetNumberOfCells()}"
 906        )
 907        return self
 908
 909    def compute_cell_vertex_count(self) -> Self:
 910        """
 911        Add to this mesh a cell data array containing the nr of vertices that a polygonal face has.
 912        """
 913        csf = vtki.new("CellSizeFilter")
 914        csf.SetInputData(self.dataset)
 915        csf.SetComputeArea(False)
 916        csf.SetComputeVolume(False)
 917        csf.SetComputeLength(False)
 918        csf.SetComputeVertexCount(True)
 919        csf.SetVertexCountArrayName("VertexCount")
 920        csf.Update()
 921        self.dataset.GetCellData().AddArray(
 922            csf.GetOutput().GetCellData().GetArray("VertexCount")
 923        )
 924        return self
 925
 926    def compute_quality(self, metric=6) -> Self:
 927        """
 928        Calculate metrics of quality for the elements of a triangular mesh.
 929        This method adds to the mesh a cell array named "Quality".
 930        See class 
 931        [vtkMeshQuality](https://vtk.org/doc/nightly/html/classvtkMeshQuality.html).
 932
 933        Arguments:
 934            metric : (int)
 935                type of available estimators are:
 936                - EDGE RATIO, 0
 937                - ASPECT RATIO, 1
 938                - RADIUS RATIO, 2
 939                - ASPECT FROBENIUS, 3
 940                - MED ASPECT FROBENIUS, 4
 941                - MAX ASPECT FROBENIUS, 5
 942                - MIN_ANGLE, 6
 943                - COLLAPSE RATIO, 7
 944                - MAX ANGLE, 8
 945                - CONDITION, 9
 946                - SCALED JACOBIAN, 10
 947                - SHEAR, 11
 948                - RELATIVE SIZE SQUARED, 12
 949                - SHAPE, 13
 950                - SHAPE AND SIZE, 14
 951                - DISTORTION, 15
 952                - MAX EDGE RATIO, 16
 953                - SKEW, 17
 954                - TAPER, 18
 955                - VOLUME, 19
 956                - STRETCH, 20
 957                - DIAGONAL, 21
 958                - DIMENSION, 22
 959                - ODDY, 23
 960                - SHEAR AND SIZE, 24
 961                - JACOBIAN, 25
 962                - WARPAGE, 26
 963                - ASPECT GAMMA, 27
 964                - AREA, 28
 965                - ASPECT BETA, 29
 966
 967        Examples:
 968            - [meshquality.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/meshquality.py)
 969
 970            ![](https://vedo.embl.es/images/advanced/meshquality.png)
 971        """
 972        qf = vtki.new("MeshQuality")
 973        qf.SetInputData(self.dataset)
 974        qf.SetTriangleQualityMeasure(metric)
 975        qf.SaveCellQualityOn()
 976        qf.Update()
 977        self._update(qf.GetOutput(), reset_locators=False)
 978        self.mapper.SetScalarModeToUseCellData()
 979        self.pipeline = OperationNode("compute_quality", parents=[self])
 980        return self
 981
 982    def count_vertices(self) -> np.ndarray:
 983        """Count the number of vertices each cell has and return it as a numpy array"""
 984        vc = vtki.new("CountVertices")
 985        vc.SetInputData(self.dataset)
 986        vc.SetOutputArrayName("VertexCount")
 987        vc.Update()
 988        varr = vc.GetOutput().GetCellData().GetArray("VertexCount")
 989        return vtk2numpy(varr)
 990
 991    def check_validity(self, tol=0) -> np.ndarray:
 992        """
 993        Return a numpy array of possible problematic faces following this convention:
 994        - Valid               =  0
 995        - WrongNumberOfPoints =  1
 996        - IntersectingEdges   =  2
 997        - IntersectingFaces   =  4
 998        - NoncontiguousEdges  =  8
 999        - Nonconvex           = 10
1000        - OrientedIncorrectly = 20
1001
1002        Arguments:
1003            tol : (float)
1004                value is used as an epsilon for floating point
1005                equality checks throughout the cell checking process.
1006        """
1007        vald = vtki.new("CellValidator")
1008        if tol:
1009            vald.SetTolerance(tol)
1010        vald.SetInputData(self.dataset)
1011        vald.Update()
1012        varr = vald.GetOutput().GetCellData().GetArray("ValidityState")
1013        return vtk2numpy(varr)
1014
1015    def compute_curvature(self, method=0) -> Self:
1016        """
1017        Add scalars to `Mesh` that contains the curvature calculated in three different ways.
1018
1019        Variable `method` can be:
1020        - 0 = gaussian
1021        - 1 = mean curvature
1022        - 2 = max curvature
1023        - 3 = min curvature
1024
1025        Example:
1026            ```python
1027            from vedo import Torus
1028            Torus().compute_curvature().add_scalarbar().show().close()
1029            ```
1030            ![](https://vedo.embl.es/images/advanced/torus_curv.png)
1031        """
1032        curve = vtki.new("Curvatures")
1033        curve.SetInputData(self.dataset)
1034        curve.SetCurvatureType(method)
1035        curve.Update()
1036        self._update(curve.GetOutput(), reset_locators=False)
1037        self.mapper.ScalarVisibilityOn()
1038        return self
1039
1040    def compute_elevation(self, low=(0, 0, 0), high=(0, 0, 1), vrange=(0, 1)) -> Self:
1041        """
1042        Add to `Mesh` a scalar array that contains distance along a specified direction.
1043
1044        Arguments:
1045            low : (list)
1046                one end of the line (small scalar values)
1047            high : (list)
1048                other end of the line (large scalar values)
1049            vrange : (list)
1050                set the range of the scalar
1051
1052        Example:
1053            ```python
1054            from vedo import Sphere
1055            s = Sphere().compute_elevation(low=(0,0,0), high=(1,1,1))
1056            s.add_scalarbar().show(axes=1).close()
1057            ```
1058            ![](https://vedo.embl.es/images/basic/compute_elevation.png)
1059        """
1060        ef = vtki.new("ElevationFilter")
1061        ef.SetInputData(self.dataset)
1062        ef.SetLowPoint(low)
1063        ef.SetHighPoint(high)
1064        ef.SetScalarRange(vrange)
1065        ef.Update()
1066        self._update(ef.GetOutput(), reset_locators=False)
1067        self.mapper.ScalarVisibilityOn()
1068        return self
1069
1070
1071    def laplacian_diffusion(self, array_name, dt, num_steps) -> Self:
1072        """
1073        Apply a diffusion process to a scalar array defined on the points of a mesh.
1074
1075        Arguments:
1076            array_name : (str)
1077                name of the array to diffuse.
1078            dt : (float)
1079                time step.
1080            num_steps : (int)
1081                number of iterations.
1082        """
1083        try:
1084            import scipy.sparse
1085            import scipy.sparse.linalg
1086        except ImportError:
1087            vedo.logger.error("scipy not found. Cannot run laplacian_diffusion()")
1088            return self
1089        
1090        def build_laplacian():
1091            rows = []
1092            cols = []
1093            data = []
1094            n_points = points.shape[0]
1095            avg_area = np.mean(areas) * 10000
1096            # print("avg_area", avg_area)
1097
1098            for triangle in cells:
1099                for i in range(3):
1100                    for j in range(i + 1, 3):
1101                        u = triangle[i]
1102                        v = triangle[j]
1103                        rows.append(u)
1104                        cols.append(v)
1105                        rows.append(v)
1106                        cols.append(u)
1107                        data.append(-1/avg_area)
1108                        data.append(-1/avg_area)
1109
1110            L = scipy.sparse.coo_matrix(
1111                (data, (rows, cols)), shape=(n_points, n_points)
1112            ).tocsc()
1113
1114            degree = -np.array(L.sum(axis=1)).flatten() # adjust the diagonal
1115            # print("degree", degree)
1116            L.setdiag(degree)
1117            return L
1118
1119        def _diffuse(u0, L, dt, num_steps):
1120            # mean_area = np.mean(areas) * 10000
1121            # print("mean_area", mean_area)
1122            mean_area = 1
1123            I = scipy.sparse.eye(L.shape[0], format="csc")
1124            A = I - (dt/mean_area) * L 
1125            u = u0
1126            for _ in range(int(num_steps)):
1127                u = A.dot(u)
1128            return u
1129
1130        self.compute_cell_size()
1131        areas = self.celldata["Area"]
1132        points = self.vertices
1133        cells = self.cells
1134        u0 = self.pointdata[array_name]
1135
1136        # Simulate diffusion
1137        L = build_laplacian()
1138        u = _diffuse(u0, L, dt, num_steps)
1139        self.pointdata[array_name] = u
1140        return self
1141
1142
1143    def subdivide(self, n=1, method=0, mel=None) -> Self:
1144        """
1145        Increase the number of vertices of a surface mesh.
1146
1147        Arguments:
1148            n : (int)
1149                number of subdivisions.
1150            method : (int)
1151                Loop(0), Linear(1), Adaptive(2), Butterfly(3), Centroid(4)
1152            mel : (float)
1153                Maximum Edge Length (applicable to Adaptive method only).
1154        """
1155        triangles = vtki.new("TriangleFilter")
1156        triangles.SetInputData(self.dataset)
1157        triangles.Update()
1158        tri_mesh = triangles.GetOutput()
1159        if method == 0:
1160            sdf = vtki.new("LoopSubdivisionFilter")
1161        elif method == 1:
1162            sdf = vtki.new("LinearSubdivisionFilter")
1163        elif method == 2:
1164            sdf = vtki.new("AdaptiveSubdivisionFilter")
1165            if mel is None:
1166                mel = self.diagonal_size() / np.sqrt(self.dataset.GetNumberOfPoints()) / n
1167            sdf.SetMaximumEdgeLength(mel)
1168        elif method == 3:
1169            sdf = vtki.new("ButterflySubdivisionFilter")
1170        elif method == 4:
1171            sdf = vtki.new("DensifyPolyData")
1172        else:
1173            vedo.logger.error(f"in subdivide() unknown method {method}")
1174            raise RuntimeError()
1175
1176        if method != 2:
1177            sdf.SetNumberOfSubdivisions(n)
1178
1179        sdf.SetInputData(tri_mesh)
1180        sdf.Update()
1181
1182        self._update(sdf.GetOutput())
1183
1184        self.pipeline = OperationNode(
1185            "subdivide",
1186            parents=[self],
1187            comment=f"#pts {self.dataset.GetNumberOfPoints()}",
1188        )
1189        return self
1190
1191
1192    def decimate(self, fraction=0.5, n=None, preserve_volume=True, regularization=0.0) -> Self:
1193        """
1194        Downsample the number of vertices in a mesh to `fraction`.
1195
1196        This filter preserves the `pointdata` of the input dataset. In previous versions
1197        of vedo, this decimation algorithm was referred to as quadric decimation.
1198
1199        Arguments:
1200            fraction : (float)
1201                the desired target of reduction.
1202            n : (int)
1203                the desired number of final points
1204                (`fraction` is recalculated based on it).
1205            preserve_volume : (bool)
1206                Decide whether to activate volume preservation which greatly
1207                reduces errors in triangle normal direction.
1208            regularization : (float)
1209                regularize the point finding algorithm so as to have better quality
1210                mesh elements at the cost of a slightly lower precision on the
1211                geometry potentially (mostly at sharp edges).
1212                Can be useful for decimating meshes that have been triangulated on noisy data.
1213
1214        Note:
1215            Setting `fraction=0.1` leaves 10% of the original number of vertices.
1216            Internally the VTK class
1217            [vtkQuadricDecimation](https://vtk.org/doc/nightly/html/classvtkQuadricDecimation.html)
1218            is used for this operation.
1219        
1220        See also: `decimate_binned()` and `decimate_pro()`.
1221        """
1222        poly = self.dataset
1223        if n:  # N = desired number of points
1224            npt = poly.GetNumberOfPoints()
1225            fraction = n / npt
1226            if fraction >= 1:
1227                return self
1228
1229        decimate = vtki.new("QuadricDecimation")
1230        decimate.SetVolumePreservation(preserve_volume)
1231        # decimate.AttributeErrorMetricOn()
1232        if regularization:
1233            decimate.SetRegularize(True)
1234            decimate.SetRegularization(regularization)
1235
1236        try:
1237            decimate.MapPointDataOn()
1238        except AttributeError:
1239            pass
1240
1241        decimate.SetTargetReduction(1 - fraction)
1242        decimate.SetInputData(poly)
1243        decimate.Update()
1244
1245        self._update(decimate.GetOutput())
1246        self.metadata["decimate_actual_fraction"] = 1 - decimate.GetActualReduction()
1247
1248        self.pipeline = OperationNode(
1249            "decimate",
1250            parents=[self],
1251            comment=f"#pts {self.dataset.GetNumberOfPoints()}",
1252        )
1253        return self
1254    
1255    def decimate_pro(
1256            self,
1257            fraction=0.5,
1258            n=None,
1259            preserve_topology=True,
1260            preserve_boundaries=True,
1261            splitting=False,
1262            splitting_angle=75,
1263            feature_angle=0,
1264            inflection_point_ratio=10,
1265            vertex_degree=0,
1266        ) -> Self:
1267        """
1268        Downsample the number of vertices in a mesh to `fraction`.
1269
1270        This filter preserves the `pointdata` of the input dataset.
1271
1272        Arguments:
1273            fraction : (float)
1274                The desired target of reduction.
1275                Setting `fraction=0.1` leaves 10% of the original number of vertices.
1276            n : (int)
1277                the desired number of final points (`fraction` is recalculated based on it).
1278            preserve_topology : (bool)
1279                If on, mesh splitting and hole elimination will not occur.
1280                This may limit the maximum reduction that may be achieved.
1281            preserve_boundaries : (bool)
1282                Turn on/off the deletion of vertices on the boundary of a mesh.
1283                Control whether mesh boundaries are preserved during decimation.
1284            feature_angle : (float)
1285                Specify the angle that defines a feature.
1286                This angle is used to define what an edge is
1287                (i.e., if the surface normal between two adjacent triangles
1288                is >= FeatureAngle, an edge exists).
1289            splitting : (bool)
1290                Turn on/off the splitting of the mesh at corners,
1291                along edges, at non-manifold points, or anywhere else a split is required.
1292                Turning splitting off will better preserve the original topology of the mesh,
1293                but you may not obtain the requested reduction.
1294            splitting_angle : (float)
1295                Specify the angle that defines a sharp edge.
1296                This angle is used to control the splitting of the mesh.
1297                A split line exists when the surface normals between two edge connected triangles
1298                are >= `splitting_angle`.
1299            inflection_point_ratio : (float)
1300                An inflection point occurs when the ratio of reduction error between two iterations
1301                is greater than or equal to the `inflection_point_ratio` value.
1302            vertex_degree : (int)
1303                If the number of triangles connected to a vertex exceeds it then the vertex will be split.
1304
1305        Note:
1306            Setting `fraction=0.1` leaves 10% of the original number of vertices
1307        
1308        See also:
1309            `decimate()` and `decimate_binned()`.
1310        """
1311        poly = self.dataset
1312        if n:  # N = desired number of points
1313            npt = poly.GetNumberOfPoints()
1314            fraction = n / npt
1315            if fraction >= 1:
1316                return self
1317
1318        decimate = vtki.new("DecimatePro")
1319        decimate.SetPreserveTopology(preserve_topology)
1320        decimate.SetBoundaryVertexDeletion(preserve_boundaries)
1321        if feature_angle:
1322            decimate.SetFeatureAngle(feature_angle)
1323        decimate.SetSplitting(splitting)
1324        decimate.SetSplitAngle(splitting_angle)
1325        decimate.SetInflectionPointRatio(inflection_point_ratio)
1326        if vertex_degree:
1327            decimate.SetDegree(vertex_degree)
1328
1329        decimate.SetTargetReduction(1 - fraction)
1330        decimate.SetInputData(poly)
1331        decimate.Update()
1332        self._update(decimate.GetOutput())
1333
1334        self.pipeline = OperationNode(
1335            "decimate_pro",
1336            parents=[self],
1337            comment=f"#pts {self.dataset.GetNumberOfPoints()}",
1338        )
1339        return self
1340    
1341    def decimate_binned(self, divisions=(), use_clustering=False) -> Self:
1342        """
1343        Downsample the number of vertices in a mesh.
1344        
1345        This filter preserves the `celldata` of the input dataset,
1346        if `use_clustering=True` also the `pointdata` will be preserved in the result.
1347
1348        Arguments:
1349            divisions : (list)
1350                number of divisions along x, y and z axes.
1351            auto_adjust : (bool)
1352                if True, the number of divisions is automatically adjusted to
1353                create more uniform cells.
1354            use_clustering : (bool)
1355                use [vtkQuadricClustering](https://vtk.org/doc/nightly/html/classvtkQuadricClustering.html)
1356                instead of 
1357                [vtkBinnedDecimation](https://vtk.org/doc/nightly/html/classvtkBinnedDecimation.html).
1358        
1359        See also: `decimate()` and `decimate_pro()`.
1360        """
1361        if use_clustering:
1362            decimate = vtki.new("QuadricClustering")
1363            decimate.CopyCellDataOn()
1364        else:
1365            decimate = vtki.new("BinnedDecimation")
1366            decimate.ProducePointDataOn()
1367            decimate.ProduceCellDataOn()
1368
1369        decimate.SetInputData(self.dataset)
1370
1371        if len(divisions) == 0:
1372            decimate.SetAutoAdjustNumberOfDivisions(1)
1373        else:
1374            decimate.SetAutoAdjustNumberOfDivisions(0)
1375            decimate.SetNumberOfDivisions(divisions)
1376        decimate.Update()
1377
1378        self._update(decimate.GetOutput())
1379        self.metadata["decimate_binned_divisions"] = decimate.GetNumberOfDivisions()
1380        self.pipeline = OperationNode(
1381            "decimate_binned",
1382            parents=[self],
1383            comment=f"#pts {self.dataset.GetNumberOfPoints()}",
1384        )
1385        return self
1386
1387    def generate_random_points(self, n: int, min_radius=0.0) -> "Points":
1388        """
1389        Generate `n` uniformly distributed random points
1390        inside the polygonal mesh.
1391
1392        A new point data array is added to the output points
1393        called "OriginalCellID" which contains the index of
1394        the cell ID in which the point was generated.
1395
1396        Arguments:
1397            n : (int)
1398                number of points to generate.
1399            min_radius: (float)
1400                impose a minimum distance between points.
1401                If `min_radius` is set to 0, the points are
1402                generated uniformly at random inside the mesh.
1403                If `min_radius` is set to a positive value,
1404                the points are generated uniformly at random
1405                inside the mesh, but points closer than `min_radius`
1406                to any other point are discarded.
1407
1408        Returns a `vedo.Points` object.
1409
1410        Note:
1411            Consider using `points.probe(msh)` or
1412            `points.interpolate_data_from(msh)`
1413            to interpolate existing mesh data onto the new points.
1414
1415        Example:
1416        ```python
1417        from vedo import *
1418        msh = Mesh(dataurl + "panther.stl").lw(2)
1419        pts = msh.generate_random_points(20000, min_radius=0.5)
1420        print("Original cell ids:", pts.pointdata["OriginalCellID"])
1421        show(pts, msh, axes=1).close()
1422        ```
1423        """
1424        cmesh = self.clone().clean().triangulate().compute_cell_size()
1425        triangles = cmesh.cells
1426        vertices = cmesh.vertices
1427        cumul = np.cumsum(cmesh.celldata["Area"])
1428
1429        out_pts = []
1430        orig_cell = []
1431        for _ in range(n):
1432            # choose a triangle based on area
1433            random_area = np.random.random() * cumul[-1]
1434            it = np.searchsorted(cumul, random_area)
1435            A, B, C = vertices[triangles[it]]
1436            # calculate the random point in the triangle
1437            r1, r2 = np.random.random(2)
1438            if r1 + r2 > 1:
1439                r1 = 1 - r1
1440                r2 = 1 - r2
1441            out_pts.append((1 - r1 - r2) * A + r1 * B + r2 * C)
1442            orig_cell.append(it)
1443        nporig_cell = np.array(orig_cell, dtype=np.uint32)
1444
1445        vpts = Points(out_pts)
1446        vpts.pointdata["OriginalCellID"] = nporig_cell
1447
1448        if min_radius > 0:
1449            vpts.subsample(min_radius, absolute=True)
1450
1451        vpts.point_size(5).color("k1")
1452        vpts.name = "RandomPoints"
1453        vpts.pipeline = OperationNode(
1454            "generate_random_points", c="#edabab", parents=[self])
1455        return vpts
1456
1457    def delete_cells(self, ids: List[int]) -> Self:
1458        """
1459        Remove cells from the mesh object by their ID.
1460        Points (vertices) are not removed (you may use `clean()` to remove those).
1461        """
1462        self.dataset.BuildLinks()
1463        for cid in ids:
1464            self.dataset.DeleteCell(cid)
1465        self.dataset.RemoveDeletedCells()
1466        self.dataset.Modified()
1467        self.mapper.Modified()
1468        self.pipeline = OperationNode(
1469            "delete_cells",
1470            parents=[self],
1471            comment=f"#cells {self.dataset.GetNumberOfCells()}",
1472        )
1473        return self
1474
1475    def delete_cells_by_point_index(self, indices: List[int]) -> Self:
1476        """
1477        Delete a list of vertices identified by any of their vertex index.
1478
1479        See also `delete_cells()`.
1480
1481        Examples:
1482            - [delete_mesh_pts.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/delete_mesh_pts.py)
1483
1484                ![](https://vedo.embl.es/images/basic/deleteMeshPoints.png)
1485        """
1486        cell_ids = vtki.vtkIdList()
1487        self.dataset.BuildLinks()
1488        n = 0
1489        for i in np.unique(indices):
1490            self.dataset.GetPointCells(i, cell_ids)
1491            for j in range(cell_ids.GetNumberOfIds()):
1492                self.dataset.DeleteCell(cell_ids.GetId(j))  # flag cell
1493                n += 1
1494
1495        self.dataset.RemoveDeletedCells()
1496        self.dataset.Modified()
1497        self.pipeline = OperationNode("delete_cells_by_point_index", parents=[self])
1498        return self
1499
1500    def collapse_edges(self, distance: float, iterations=1) -> Self:
1501        """
1502        Collapse mesh edges so that are all above `distance`.
1503        
1504        Example:
1505            ```python
1506            from vedo import *
1507            np.random.seed(2)
1508            grid1 = Grid().add_gaussian_noise(0.8).triangulate().lw(1)
1509            grid1.celldata['scalar'] = grid1.cell_centers[:,1]
1510            grid2 = grid1.clone().collapse_edges(0.1)
1511            show(grid1, grid2, N=2, axes=1)
1512            ```
1513        """
1514        for _ in range(iterations):
1515            medges = self.edges
1516            pts = self.vertices
1517            newpts = np.array(pts)
1518            moved = []
1519            for e in medges:
1520                if len(e) == 2:
1521                    id0, id1 = e
1522                    p0, p1 = pts[id0], pts[id1]
1523                    if (np.linalg.norm(p1-p0) < distance 
1524                        and id0 not in moved
1525                        and id1 not in moved
1526                    ):
1527                        p = (p0 + p1) / 2
1528                        newpts[id0] = p
1529                        newpts[id1] = p
1530                        moved += [id0, id1]
1531            self.vertices = newpts
1532            cpd = vtki.new("CleanPolyData")
1533            cpd.ConvertLinesToPointsOff()
1534            cpd.ConvertPolysToLinesOff()
1535            cpd.ConvertStripsToPolysOff()
1536            cpd.SetInputData(self.dataset)
1537            cpd.Update()
1538            self._update(cpd.GetOutput())
1539
1540        self.pipeline = OperationNode(
1541            "collapse_edges",
1542            parents=[self],
1543            comment=f"#pts {self.dataset.GetNumberOfPoints()}",
1544        )
1545        return self
1546
1547    def adjacency_list(self) -> List[set]:
1548        """
1549        Computes the adjacency list for mesh edge-graph.
1550
1551        Returns: 
1552            a list with i-th entry being the set if indices of vertices connected by an edge to i-th vertex
1553        """
1554        inc = [set()] * self.nvertices
1555        for cell in self.cells:
1556            nc = len(cell)
1557            if nc > 1:
1558                for i in range(nc-1):
1559                    ci = cell[i]
1560                    inc[ci] = inc[ci].union({cell[i-1], cell[i+1]})
1561        return inc
1562
1563    def graph_ball(self, index, n: int) -> set:
1564        """
1565        Computes the ball of radius `n` in the mesh' edge-graph metric centred in vertex `index`.
1566
1567        Arguments:
1568            index : (int)
1569                index of the vertex
1570            n : (int)
1571                radius in the graph metric
1572
1573        Returns:
1574            the set of indices of the vertices which are at most `n` edges from vertex `index`.
1575        """
1576        if n == 0:
1577            return {index}
1578        else:
1579            al = self.adjacency_list()
1580            ball = {index}
1581            i = 0
1582            while i < n and len(ball) < self.nvertices:
1583                for v in ball:
1584                    ball = ball.union(al[v])
1585                i += 1
1586            return ball
1587
1588    def smooth(self, niter=15, pass_band=0.1, edge_angle=15, feature_angle=60, boundary=False) -> Self:
1589        """
1590        Adjust mesh point positions using the so-called "Windowed Sinc" method.
1591
1592        Arguments:
1593            niter : (int)
1594                number of iterations.
1595            pass_band : (float)
1596                set the pass_band value for the windowed sinc filter.
1597            edge_angle : (float)
1598                edge angle to control smoothing along edges (either interior or boundary).
1599            feature_angle : (float)
1600                specifies the feature angle for sharp edge identification.
1601            boundary : (bool)
1602                specify if boundary should also be smoothed or kept unmodified
1603
1604        Examples:
1605            - [mesh_smoother1.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/mesh_smoother1.py)
1606
1607            ![](https://vedo.embl.es/images/advanced/mesh_smoother2.png)
1608        """
1609        cl = vtki.new("CleanPolyData")
1610        cl.SetInputData(self.dataset)
1611        cl.Update()
1612        smf = vtki.new("WindowedSincPolyDataFilter")
1613        smf.SetInputData(cl.GetOutput())
1614        smf.SetNumberOfIterations(niter)
1615        smf.SetEdgeAngle(edge_angle)
1616        smf.SetFeatureAngle(feature_angle)
1617        smf.SetPassBand(pass_band)
1618        smf.NormalizeCoordinatesOn()
1619        smf.NonManifoldSmoothingOn()
1620        smf.FeatureEdgeSmoothingOn()
1621        smf.SetBoundarySmoothing(boundary)
1622        smf.Update()
1623
1624        self._update(smf.GetOutput())
1625
1626        self.pipeline = OperationNode(
1627            "smooth", parents=[self], comment=f"#pts {self.dataset.GetNumberOfPoints()}"
1628        )
1629        return self
1630
1631    def fill_holes(self, size=None) -> Self:
1632        """
1633        Identifies and fills holes in the input mesh.
1634        Holes are identified by locating boundary edges, linking them together
1635        into loops, and then triangulating the resulting loops.
1636
1637        Arguments:
1638            size : (float)
1639                Approximate limit to the size of the hole that can be filled.
1640
1641        Examples:
1642            - [fillholes.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/fillholes.py)
1643        """
1644        fh = vtki.new("FillHolesFilter")
1645        if not size:
1646            mb = self.diagonal_size()
1647            size = mb / 10
1648        fh.SetHoleSize(size)
1649        fh.SetInputData(self.dataset)
1650        fh.Update()
1651
1652        self._update(fh.GetOutput())
1653
1654        self.pipeline = OperationNode(
1655            "fill_holes",
1656            parents=[self],
1657            comment=f"#pts {self.dataset.GetNumberOfPoints()}",
1658        )
1659        return self
1660
1661    def contains(self, point: tuple, tol=1e-05) -> bool:
1662        """
1663        Return True if point is inside a polydata closed surface.
1664        
1665        Note:
1666            if you have many points to check use `inside_points()` instead.
1667        
1668        Example:
1669            ```python
1670            from vedo import *
1671            s = Sphere().c('green5').alpha(0.5)
1672            pt  = [0.1, 0.2, 0.3]
1673            print("Sphere contains", pt, s.contains(pt))
1674            show(s, Point(pt), axes=1).close()
1675            ```      
1676        """
1677        points = vtki.vtkPoints()
1678        points.InsertNextPoint(point)
1679        poly = vtki.vtkPolyData()
1680        poly.SetPoints(points)
1681        sep = vtki.new("SelectEnclosedPoints")
1682        sep.SetTolerance(tol)
1683        sep.CheckSurfaceOff()
1684        sep.SetInputData(poly)
1685        sep.SetSurfaceData(self.dataset)
1686        sep.Update()
1687        return bool(sep.IsInside(0))
1688
1689    def inside_points(self, pts: Union["Points", list], invert=False, tol=1e-05, return_ids=False) -> Union["Points", np.ndarray]:
1690        """
1691        Return the point cloud that is inside mesh surface as a new Points object.
1692
1693        If return_ids is True a list of IDs is returned and in addition input points
1694        are marked by a pointdata array named "IsInside".
1695
1696        Example:
1697            `print(pts.pointdata["IsInside"])`
1698
1699        Examples:
1700            - [pca_ellipsoid.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/pca_ellipsoid.py)
1701
1702            ![](https://vedo.embl.es/images/basic/pca.png)
1703        """
1704        if isinstance(pts, Points):
1705            poly = pts.dataset
1706            ptsa = pts.vertices
1707        else:
1708            ptsa = np.asarray(pts)
1709            vpoints = vtki.vtkPoints()
1710            vpoints.SetData(numpy2vtk(ptsa, dtype=np.float32))
1711            poly = vtki.vtkPolyData()
1712            poly.SetPoints(vpoints)
1713
1714        sep = vtki.new("SelectEnclosedPoints")
1715        # sep = vtki.new("ExtractEnclosedPoints()
1716        sep.SetTolerance(tol)
1717        sep.SetInputData(poly)
1718        sep.SetSurfaceData(self.dataset)
1719        sep.SetInsideOut(invert)
1720        sep.Update()
1721
1722        varr = sep.GetOutput().GetPointData().GetArray("SelectedPoints")
1723        mask = vtk2numpy(varr).astype(bool)
1724        ids = np.array(range(len(ptsa)), dtype=int)[mask]
1725
1726        if isinstance(pts, Points):
1727            varr.SetName("IsInside")
1728            pts.dataset.GetPointData().AddArray(varr)
1729
1730        if return_ids:
1731            return ids
1732
1733        pcl = Points(ptsa[ids])
1734        pcl.name = "InsidePoints"
1735
1736        pcl.pipeline = OperationNode(
1737            "inside_points",
1738            parents=[self, ptsa],
1739            comment=f"#pts {pcl.dataset.GetNumberOfPoints()}",
1740        )
1741        return pcl
1742
1743    def boundaries(
1744        self,
1745        boundary_edges=True,
1746        manifold_edges=False,
1747        non_manifold_edges=False,
1748        feature_angle=None,
1749        return_point_ids=False,
1750        return_cell_ids=False,
1751        cell_edge=False,
1752    ) -> Union[Self, np.ndarray]:
1753        """
1754        Return the boundary lines of an input mesh.
1755        Check also `vedo.core.CommonAlgorithms.mark_boundaries()` method.
1756
1757        Arguments:
1758            boundary_edges : (bool)
1759                Turn on/off the extraction of boundary edges.
1760            manifold_edges : (bool)
1761                Turn on/off the extraction of manifold edges.
1762            non_manifold_edges : (bool)
1763                Turn on/off the extraction of non-manifold edges.
1764            feature_angle : (bool)
1765                Specify the min angle btw 2 faces for extracting edges.
1766            return_point_ids : (bool)
1767                return a numpy array of point indices
1768            return_cell_ids : (bool)
1769                return a numpy array of cell indices
1770            cell_edge : (bool)
1771                set to `True` if a cell need to share an edge with
1772                the boundary line, or `False` if a single vertex is enough
1773
1774        Examples:
1775            - [boundaries.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/boundaries.py)
1776
1777            ![](https://vedo.embl.es/images/basic/boundaries.png)
1778        """
1779        fe = vtki.new("FeatureEdges")
1780        fe.SetBoundaryEdges(boundary_edges)
1781        fe.SetNonManifoldEdges(non_manifold_edges)
1782        fe.SetManifoldEdges(manifold_edges)
1783        try:
1784            fe.SetPassLines(True) # vtk9.2
1785        except AttributeError:
1786            pass
1787        fe.ColoringOff()
1788        fe.SetFeatureEdges(False)
1789        if feature_angle is not None:
1790            fe.SetFeatureEdges(True)
1791            fe.SetFeatureAngle(feature_angle)
1792
1793        if return_point_ids or return_cell_ids:
1794            idf = vtki.new("IdFilter")
1795            idf.SetInputData(self.dataset)
1796            idf.SetPointIdsArrayName("BoundaryIds")
1797            idf.SetPointIds(True)
1798            idf.Update()
1799
1800            fe.SetInputData(idf.GetOutput())
1801            fe.Update()
1802
1803            vid = fe.GetOutput().GetPointData().GetArray("BoundaryIds")
1804            npid = vtk2numpy(vid).astype(int)
1805
1806            if return_point_ids:
1807                return npid
1808
1809            if return_cell_ids:
1810                n = 1 if cell_edge else 0
1811                inface = []
1812                for i, face in enumerate(self.cells):
1813                    # isin = np.any([vtx in npid for vtx in face])
1814                    isin = 0
1815                    for vtx in face:
1816                        isin += int(vtx in npid)
1817                        if isin > n:
1818                            break
1819                    if isin > n:
1820                        inface.append(i)
1821                return np.array(inface).astype(int)
1822
1823            return self
1824
1825        else:
1826
1827            fe.SetInputData(self.dataset)
1828            fe.Update()
1829            msh = Mesh(fe.GetOutput(), c="p").lw(5).lighting("off")
1830            msh.name = "MeshBoundaries"
1831
1832            msh.pipeline = OperationNode(
1833                "boundaries",
1834                parents=[self],
1835                shape="octagon",
1836                comment=f"#pts {msh.dataset.GetNumberOfPoints()}",
1837            )
1838            return msh
1839
1840    def imprint(self, loopline, tol=0.01) -> Self:
1841        """
1842        Imprint the contact surface of one object onto another surface.
1843
1844        Arguments:
1845            loopline : (vedo.Line)
1846                a Line object to be imprinted onto the mesh.
1847            tol : (float)
1848                projection tolerance which controls how close the imprint
1849                surface must be to the target.
1850
1851        Example:
1852            ```python
1853            from vedo import *
1854            grid = Grid()#.triangulate()
1855            circle = Circle(r=0.3, res=24).pos(0.11,0.12)
1856            line = Line(circle, closed=True, lw=4, c='r4')
1857            grid.imprint(line)
1858            show(grid, line, axes=1).close()
1859            ```
1860            ![](https://vedo.embl.es/images/feats/imprint.png)
1861        """
1862        loop = vtki.new("ContourLoopExtraction")
1863        loop.SetInputData(loopline.dataset)
1864        loop.Update()
1865
1866        clean_loop = vtki.new("CleanPolyData")
1867        clean_loop.SetInputData(loop.GetOutput())
1868        clean_loop.Update()
1869
1870        imp = vtki.new("ImprintFilter")
1871        imp.SetTargetData(self.dataset)
1872        imp.SetImprintData(clean_loop.GetOutput())
1873        imp.SetTolerance(tol)
1874        imp.BoundaryEdgeInsertionOn()
1875        imp.TriangulateOutputOn()
1876        imp.Update()
1877
1878        self._update(imp.GetOutput())
1879
1880        self.pipeline = OperationNode(
1881            "imprint",
1882            parents=[self],
1883            comment=f"#pts {self.dataset.GetNumberOfPoints()}",
1884        )
1885        return self
1886
1887    def connected_vertices(self, index: int) -> List[int]:
1888        """Find all vertices connected to an input vertex specified by its index.
1889
1890        Examples:
1891            - [connected_vtx.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/connected_vtx.py)
1892
1893            ![](https://vedo.embl.es/images/basic/connVtx.png)
1894        """
1895        poly = self.dataset
1896
1897        cell_idlist = vtki.vtkIdList()
1898        poly.GetPointCells(index, cell_idlist)
1899
1900        idxs = []
1901        for i in range(cell_idlist.GetNumberOfIds()):
1902            point_idlist = vtki.vtkIdList()
1903            poly.GetCellPoints(cell_idlist.GetId(i), point_idlist)
1904            for j in range(point_idlist.GetNumberOfIds()):
1905                idj = point_idlist.GetId(j)
1906                if idj == index:
1907                    continue
1908                if idj in idxs:
1909                    continue
1910                idxs.append(idj)
1911
1912        return idxs
1913
1914    def extract_cells(self, ids: List[int]) -> Self:
1915        """
1916        Extract a subset of cells from a mesh and return it as a new mesh.
1917        """
1918        selectCells = vtki.new("SelectionNode")
1919        selectCells.SetFieldType(vtki.get_class("SelectionNode").CELL)
1920        selectCells.SetContentType(vtki.get_class("SelectionNode").INDICES)
1921        idarr = vtki.vtkIdTypeArray()
1922        idarr.SetNumberOfComponents(1)
1923        idarr.SetNumberOfValues(len(ids))
1924        for i, v in enumerate(ids):
1925            idarr.SetValue(i, v)
1926        selectCells.SetSelectionList(idarr)
1927
1928        selection = vtki.new("Selection")
1929        selection.AddNode(selectCells)
1930
1931        extractSelection = vtki.new("ExtractSelection")
1932        extractSelection.SetInputData(0, self.dataset)
1933        extractSelection.SetInputData(1, selection)
1934        extractSelection.Update()
1935
1936        gf = vtki.new("GeometryFilter")
1937        gf.SetInputData(extractSelection.GetOutput())
1938        gf.Update()
1939        msh = Mesh(gf.GetOutput())
1940        msh.copy_properties_from(self)
1941        return msh
1942
1943    def connected_cells(self, index: int, return_ids=False) -> Union[Self, List[int]]:
1944        """Find all cellls connected to an input vertex specified by its index."""
1945
1946        # Find all cells connected to point index
1947        dpoly = self.dataset
1948        idlist = vtki.vtkIdList()
1949        dpoly.GetPointCells(index, idlist)
1950
1951        ids = vtki.vtkIdTypeArray()
1952        ids.SetNumberOfComponents(1)
1953        rids = []
1954        for k in range(idlist.GetNumberOfIds()):
1955            cid = idlist.GetId(k)
1956            ids.InsertNextValue(cid)
1957            rids.append(int(cid))
1958        if return_ids:
1959            return rids
1960
1961        selection_node = vtki.new("SelectionNode")
1962        selection_node.SetFieldType(vtki.get_class("SelectionNode").CELL)
1963        selection_node.SetContentType(vtki.get_class("SelectionNode").INDICES)
1964        selection_node.SetSelectionList(ids)
1965        selection = vtki.new("Selection")
1966        selection.AddNode(selection_node)
1967        extractSelection = vtki.new("ExtractSelection")
1968        extractSelection.SetInputData(0, dpoly)
1969        extractSelection.SetInputData(1, selection)
1970        extractSelection.Update()
1971        gf = vtki.new("GeometryFilter")
1972        gf.SetInputData(extractSelection.GetOutput())
1973        gf.Update()
1974        return Mesh(gf.GetOutput()).lw(1)
1975
1976    def silhouette(self, direction=None, border_edges=True, feature_angle=False) -> Self:
1977        """
1978        Return a new line `Mesh` which corresponds to the outer `silhouette`
1979        of the input as seen along a specified `direction`, this can also be
1980        a `vtkCamera` object.
1981
1982        Arguments:
1983            direction : (list)
1984                viewpoint direction vector.
1985                If `None` this is guessed by looking at the minimum
1986                of the sides of the bounding box.
1987            border_edges : (bool)
1988                enable or disable generation of border edges
1989            feature_angle : (float)
1990                minimal angle for sharp edges detection.
1991                If set to `False` the functionality is disabled.
1992
1993        Examples:
1994            - [silhouette1.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/silhouette1.py)
1995
1996            ![](https://vedo.embl.es/images/basic/silhouette1.png)
1997        """
1998        sil = vtki.new("PolyDataSilhouette")
1999        sil.SetInputData(self.dataset)
2000        sil.SetBorderEdges(border_edges)
2001        if feature_angle is False:
2002            sil.SetEnableFeatureAngle(0)
2003        else:
2004            sil.SetEnableFeatureAngle(1)
2005            sil.SetFeatureAngle(feature_angle)
2006
2007        if direction is None and vedo.plotter_instance and vedo.plotter_instance.camera:
2008            sil.SetCamera(vedo.plotter_instance.camera)
2009            m = Mesh()
2010            m.mapper.SetInputConnection(sil.GetOutputPort())
2011
2012        elif isinstance(direction, vtki.vtkCamera):
2013            sil.SetCamera(direction)
2014            m = Mesh()
2015            m.mapper.SetInputConnection(sil.GetOutputPort())
2016
2017        elif direction == "2d":
2018            sil.SetVector(3.4, 4.5, 5.6)  # random
2019            sil.SetDirectionToSpecifiedVector()
2020            sil.Update()
2021            m = Mesh(sil.GetOutput())
2022
2023        elif is_sequence(direction):
2024            sil.SetVector(direction)
2025            sil.SetDirectionToSpecifiedVector()
2026            sil.Update()
2027            m = Mesh(sil.GetOutput())
2028        else:
2029            vedo.logger.error(f"in silhouette() unknown direction type {type(direction)}")
2030            vedo.logger.error("first render the scene with show() or specify camera/direction")
2031            return self
2032
2033        m.lw(2).c((0, 0, 0)).lighting("off")
2034        m.mapper.SetResolveCoincidentTopologyToPolygonOffset()
2035        m.pipeline = OperationNode("silhouette", parents=[self])
2036        m.name = "Silhouette"
2037        return m
2038
2039    def isobands(self, n=10, vmin=None, vmax=None) -> Self:
2040        """
2041        Return a new `Mesh` representing the isobands of the active scalars.
2042        This is a new mesh where the scalar is now associated to cell faces and
2043        used to colorize the mesh.
2044
2045        Arguments:
2046            n : (int)
2047                number of isobands in the range
2048            vmin : (float)
2049                minimum of the range
2050            vmax : (float)
2051                maximum of the range
2052
2053        Examples:
2054            - [isolines.py](https://github.com/marcomusy/vedo/tree/master/examples/pyplot/isolines.py)
2055        """
2056        r0, r1 = self.dataset.GetScalarRange()
2057        if vmin is None:
2058            vmin = r0
2059        if vmax is None:
2060            vmax = r1
2061
2062        # --------------------------------
2063        bands = []
2064        dx = (vmax - vmin) / float(n)
2065        b = [vmin, vmin + dx / 2.0, vmin + dx]
2066        i = 0
2067        while i < n:
2068            bands.append(b)
2069            b = [b[0] + dx, b[1] + dx, b[2] + dx]
2070            i += 1
2071
2072        # annotate, use the midpoint of the band as the label
2073        lut = self.mapper.GetLookupTable()
2074        labels = []
2075        for b in bands:
2076            labels.append("{:4.2f}".format(b[1]))
2077        values = vtki.vtkVariantArray()
2078        for la in labels:
2079            values.InsertNextValue(vtki.vtkVariant(la))
2080        for i in range(values.GetNumberOfTuples()):
2081            lut.SetAnnotation(i, values.GetValue(i).ToString())
2082
2083        bcf = vtki.new("BandedPolyDataContourFilter")
2084        bcf.SetInputData(self.dataset)
2085        # Use either the minimum or maximum value for each band.
2086        for i, band in enumerate(bands):
2087            bcf.SetValue(i, band[2])
2088        # We will use an indexed lookup table.
2089        bcf.SetScalarModeToIndex()
2090        bcf.GenerateContourEdgesOff()
2091        bcf.Update()
2092        bcf.GetOutput().GetCellData().GetScalars().SetName("IsoBands")
2093
2094        m1 = Mesh(bcf.GetOutput()).compute_normals(cells=True)
2095        m1.mapper.SetLookupTable(lut)
2096        m1.mapper.SetScalarRange(lut.GetRange())
2097        m1.pipeline = OperationNode("isobands", parents=[self])
2098        m1.name = "IsoBands"
2099        return m1
2100
2101    def isolines(self, n=10, vmin=None, vmax=None) -> Self:
2102        """
2103        Return a new `Mesh` representing the isolines of the active scalars.
2104
2105        Arguments:
2106            n : (int)
2107                number of isolines in the range
2108            vmin : (float)
2109                minimum of the range
2110            vmax : (float)
2111                maximum of the range
2112
2113        Examples:
2114            - [isolines.py](https://github.com/marcomusy/vedo/tree/master/examples/pyplot/isolines.py)
2115
2116            ![](https://vedo.embl.es/images/pyplot/isolines.png)
2117        """
2118        bcf = vtki.new("ContourFilter")
2119        bcf.SetInputData(self.dataset)
2120        r0, r1 = self.dataset.GetScalarRange()
2121        if vmin is None:
2122            vmin = r0
2123        if vmax is None:
2124            vmax = r1
2125        bcf.GenerateValues(n, vmin, vmax)
2126        bcf.Update()
2127        sf = vtki.new("Stripper")
2128        sf.SetJoinContiguousSegments(True)
2129        sf.SetInputData(bcf.GetOutput())
2130        sf.Update()
2131        cl = vtki.new("CleanPolyData")
2132        cl.SetInputData(sf.GetOutput())
2133        cl.Update()
2134        msh = Mesh(cl.GetOutput(), c="k").lighting("off")
2135        msh.mapper.SetResolveCoincidentTopologyToPolygonOffset()
2136        msh.pipeline = OperationNode("isolines", parents=[self])
2137        msh.name = "IsoLines"
2138        return msh
2139
2140    def extrude(self, zshift=1.0, direction=(), rotation=0.0, dr=0.0, cap=True, res=1) -> Self:
2141        """
2142        Sweep a polygonal data creating a "skirt" from free edges and lines, and lines from vertices.
2143        The input dataset is swept around the z-axis to create new polygonal primitives.
2144        For example, sweeping a line results in a cylindrical shell, and sweeping a circle creates a torus.
2145
2146        You can control whether the sweep of a 2D object (i.e., polygon or triangle strip)
2147        is capped with the generating geometry.
2148        Also, you can control the angle of rotation, and whether translation along the z-axis
2149        is performed along with the rotation. (Translation is useful for creating "springs").
2150        You also can adjust the radius of the generating geometry using the "dR" keyword.
2151
2152        The skirt is generated by locating certain topological features.
2153        Free edges (edges of polygons or triangle strips only used by one polygon or triangle strips)
2154        generate surfaces. This is true also of lines or polylines. Vertices generate lines.
2155
2156        This filter can be used to model axisymmetric objects like cylinders, bottles, and wine glasses;
2157        or translational/rotational symmetric objects like springs or corkscrews.
2158
2159        Arguments:
2160            zshift : (float)
2161                shift along z axis.
2162            direction : (list)
2163                extrusion direction in the xy plane. 
2164                note that zshift is forced to be the 3rd component of direction,
2165                which is therefore ignored.
2166            rotation : (float)
2167                set the angle of rotation.
2168            dr : (float)
2169                set the radius variation in absolute units.
2170            cap : (bool)
2171                enable or disable capping.
2172            res : (int)
2173                set the resolution of the generating geometry.
2174
2175        Warning:
2176            Some polygonal objects have no free edges (e.g., sphere). When swept, this will result
2177            in two separate surfaces if capping is on, or no surface if capping is off.
2178
2179        Examples:
2180            - [extrude.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/extrude.py)
2181
2182            ![](https://vedo.embl.es/images/basic/extrude.png)
2183        """
2184        rf = vtki.new("RotationalExtrusionFilter")
2185        # rf = vtki.new("LinearExtrusionFilter")
2186        rf.SetInputData(self.dataset)  # must not be transformed
2187        rf.SetResolution(res)
2188        rf.SetCapping(cap)
2189        rf.SetAngle(rotation)
2190        rf.SetTranslation(zshift)
2191        rf.SetDeltaRadius(dr)
2192        rf.Update()
2193
2194        # convert triangle strips to polygonal data
2195        tris = vtki.new("TriangleFilter")
2196        tris.SetInputData(rf.GetOutput())
2197        tris.Update()
2198
2199        m = Mesh(tris.GetOutput())
2200
2201        if len(direction) > 1:
2202            p = self.pos()
2203            LT = vedo.LinearTransform()
2204            LT.translate(-p)
2205            LT.concatenate([
2206                [1, 0, direction[0]],
2207                [0, 1, direction[1]],
2208                [0, 0, 1]
2209            ])
2210            LT.translate(p)
2211            m.apply_transform(LT)
2212
2213        m.copy_properties_from(self).flat().lighting("default")
2214        m.pipeline = OperationNode(
2215            "extrude", parents=[self], 
2216            comment=f"#pts {m.dataset.GetNumberOfPoints()}"
2217        )
2218        m.name = "ExtrudedMesh"
2219        return m
2220
2221    def extrude_and_trim_with(
2222            self,
2223            surface: "Mesh",
2224            direction=(),
2225            strategy="all",
2226            cap=True,
2227            cap_strategy="max",
2228    ) -> Self:
2229        """
2230        Extrude a Mesh and trim it with an input surface mesh.
2231
2232        Arguments:
2233            surface : (Mesh)
2234                the surface mesh to trim with.
2235            direction : (list)
2236                extrusion direction in the xy plane.
2237            strategy : (str)
2238                either "boundary_edges" or "all_edges".
2239            cap : (bool)
2240                enable or disable capping.
2241            cap_strategy : (str)
2242                either "intersection", "minimum_distance", "maximum_distance", "average_distance".
2243
2244        The input Mesh is swept along a specified direction forming a "skirt"
2245        from the boundary edges 2D primitives (i.e., edges used by only one polygon);
2246        and/or from vertices and lines.
2247        The extent of the sweeping is limited by a second input: defined where
2248        the sweep intersects a user-specified surface.
2249
2250        Capping of the extrusion can be enabled.
2251        In this case the input, generating primitive is copied inplace as well
2252        as to the end of the extrusion skirt.
2253        (See warnings below on what happens if the intersecting sweep does not
2254        intersect, or partially intersects the trim surface.)
2255
2256        Note that this method operates in two fundamentally different modes
2257        based on the extrusion strategy. 
2258        If the strategy is "boundary_edges", then only the boundary edges of the input's
2259        2D primitives are extruded (verts and lines are extruded to generate lines and quads).
2260        However, if the extrusions strategy is "all_edges", then every edge of the 2D primitives
2261        is used to sweep out a quadrilateral polygon (again verts and lines are swept to produce lines and quads).
2262
2263        Warning:
2264            The extrusion direction is assumed to define an infinite line.
2265            The intersection with the trim surface is along a ray from the - to + direction,
2266            however only the first intersection is taken.
2267            Some polygonal objects have no free edges (e.g., sphere). When swept, this will result in two separate
2268            surfaces if capping is on and "boundary_edges" enabled,
2269            or no surface if capping is off and "boundary_edges" is enabled.
2270            If all the extrusion lines emanating from an extruding primitive do not intersect the trim surface,
2271            then no output for that primitive will be generated. In extreme cases, it is possible that no output
2272            whatsoever will be generated.
2273        
2274        Example:
2275            ```python
2276            from vedo import *
2277            sphere = Sphere([-1,0,4]).rotate_x(25).wireframe().color('red5')
2278            circle = Circle([0,0,0], r=2, res=100).color('b6')
2279            extruded_circle = circle.extrude_and_trim_with(
2280                sphere, 
2281                direction=[0,-0.2,1],
2282                strategy="bound",
2283                cap=True,
2284                cap_strategy="intersection",
2285            )
2286            circle.lw(3).color("tomato").shift(dz=-0.1)
2287            show(circle, sphere, extruded_circle, axes=1).close()
2288            ```
2289        """
2290        trimmer = vtki.new("TrimmedExtrusionFilter")
2291        trimmer.SetInputData(self.dataset)
2292        trimmer.SetCapping(cap)
2293        trimmer.SetExtrusionDirection(direction)
2294        trimmer.SetTrimSurfaceData(surface.dataset)
2295        if "bound" in strategy:
2296            trimmer.SetExtrusionStrategyToBoundaryEdges()
2297        elif "all" in strategy:
2298            trimmer.SetExtrusionStrategyToAllEdges()
2299        else:
2300            vedo.logger.warning(f"extrude_and_trim(): unknown strategy {strategy}")
2301        # print (trimmer.GetExtrusionStrategy())
2302        
2303        if "intersect" in cap_strategy:
2304            trimmer.SetCappingStrategyToIntersection()
2305        elif "min" in cap_strategy:
2306            trimmer.SetCappingStrategyToMinimumDistance()
2307        elif "max" in cap_strategy:
2308            trimmer.SetCappingStrategyToMaximumDistance()
2309        elif "ave" in cap_strategy:
2310            trimmer.SetCappingStrategyToAverageDistance()
2311        else:
2312            vedo.logger.warning(f"extrude_and_trim(): unknown cap_strategy {cap_strategy}")
2313        # print (trimmer.GetCappingStrategy())
2314
2315        trimmer.Update()
2316
2317        m = Mesh(trimmer.GetOutput())
2318        m.copy_properties_from(self).flat().lighting("default")
2319        m.pipeline = OperationNode(
2320            "extrude_and_trim", parents=[self, surface],
2321            comment=f"#pts {m.dataset.GetNumberOfPoints()}"
2322        )
2323        m.name = "ExtrudedAndTrimmedMesh"
2324        return m
2325
2326    def split(
2327        self, maxdepth=1000, flag=False, must_share_edge=False, sort_by_area=True
2328    ) -> List[Self]:
2329        """
2330        Split a mesh by connectivity and order the pieces by increasing area.
2331
2332        Arguments:
2333            maxdepth : (int)
2334                only consider this maximum number of mesh parts.
2335            flag : (bool)
2336                if set to True return the same single object,
2337                but add a "RegionId" array to flag the mesh subparts
2338            must_share_edge : (bool)
2339                if True, mesh regions that only share single points will be split.
2340            sort_by_area : (bool)
2341                if True, sort the mesh parts by decreasing area.
2342
2343        Examples:
2344            - [splitmesh.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/splitmesh.py)
2345
2346            ![](https://vedo.embl.es/images/advanced/splitmesh.png)
2347        """
2348        pd = self.dataset
2349        if must_share_edge:
2350            if pd.GetNumberOfPolys() == 0:
2351                vedo.logger.warning("in split(): no polygons found. Skip.")
2352                return [self]
2353            cf = vtki.new("PolyDataEdgeConnectivityFilter")
2354            cf.BarrierEdgesOff()
2355        else:
2356            cf = vtki.new("PolyDataConnectivityFilter")
2357
2358        cf.SetInputData(pd)
2359        cf.SetExtractionModeToAllRegions()
2360        cf.SetColorRegions(True)
2361        cf.Update()
2362        out = cf.GetOutput()
2363
2364        if not out.GetNumberOfPoints():
2365            return [self]
2366
2367        if flag:
2368            self.pipeline = OperationNode("split mesh", parents=[self])
2369            self._update(out)
2370            return [self]
2371
2372        msh = Mesh(out)
2373        if must_share_edge:
2374            arr = msh.celldata["RegionId"]
2375            on = "cells"
2376        else:
2377            arr = msh.pointdata["RegionId"]
2378            on = "points"
2379
2380        alist = []
2381        for t in range(max(arr) + 1):
2382            if t == maxdepth:
2383                break
2384            suba = msh.clone().threshold("RegionId", t, t, on=on)
2385            if sort_by_area:
2386                area = suba.area()
2387            else:
2388                area = 0  # dummy
2389            suba.name = "MeshRegion" + str(t)
2390            alist.append([suba, area])
2391
2392        if sort_by_area:
2393            alist.sort(key=lambda x: x[1])
2394            alist.reverse()
2395
2396        blist = []
2397        for i, l in enumerate(alist):
2398            l[0].color(i + 1).phong()
2399            l[0].mapper.ScalarVisibilityOff()
2400            blist.append(l[0])
2401            if i < 10:
2402                l[0].pipeline = OperationNode(
2403                    f"split mesh {i}",
2404                    parents=[self],
2405                    comment=f"#pts {l[0].dataset.GetNumberOfPoints()}",
2406                )
2407        return blist
2408
2409    def extract_largest_region(self) -> Self:
2410        """
2411        Extract the largest connected part of a mesh and discard all the smaller pieces.
2412
2413        Examples:
2414            - [largestregion.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/largestregion.py)
2415        """
2416        conn = vtki.new("PolyDataConnectivityFilter")
2417        conn.SetExtractionModeToLargestRegion()
2418        conn.ScalarConnectivityOff()
2419        conn.SetInputData(self.dataset)
2420        conn.Update()
2421
2422        m = Mesh(conn.GetOutput())
2423        m.copy_properties_from(self)
2424        m.pipeline = OperationNode(
2425            "extract_largest_region",
2426            parents=[self],
2427            comment=f"#pts {m.dataset.GetNumberOfPoints()}",
2428        )
2429        m.name = "MeshLargestRegion"
2430        return m
2431
2432    def boolean(self, operation: str, mesh2, method=0, tol=None) -> Self:
2433        """Volumetric union, intersection and subtraction of surfaces.
2434
2435        Use `operation` for the allowed operations `['plus', 'intersect', 'minus']`.
2436
2437        Two possible algorithms are available.
2438        Setting `method` to 0 (the default) uses the boolean operation algorithm
2439        written by Cory Quammen, Chris Weigle, and Russ Taylor (https://doi.org/10.54294/216g01);
2440        setting `method` to 1 will use the "loop" boolean algorithm
2441        written by Adam Updegrove (https://doi.org/10.1016/j.advengsoft.2016.01.015).
2442
2443        Use `tol` to specify the absolute tolerance used to determine
2444        when the distance between two points is considered to be zero (defaults to 1e-6).
2445
2446        Example:
2447            - [boolean.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/boolean.py)
2448
2449            ![](https://vedo.embl.es/images/basic/boolean.png)
2450        """
2451        if method == 0:
2452            bf = vtki.new("BooleanOperationPolyDataFilter")
2453        elif method == 1:
2454            bf = vtki.new("LoopBooleanPolyDataFilter")
2455        else:
2456            raise ValueError(f"Unknown method={method}")
2457
2458        poly1 = self.compute_normals().dataset
2459        poly2 = mesh2.compute_normals().dataset
2460
2461        if operation.lower() in ("plus", "+"):
2462            bf.SetOperationToUnion()
2463        elif operation.lower() == "intersect":
2464            bf.SetOperationToIntersection()
2465        elif operation.lower() in ("minus", "-"):
2466            bf.SetOperationToDifference()
2467
2468        if tol:
2469            bf.SetTolerance(tol)
2470
2471        bf.SetInputData(0, poly1)
2472        bf.SetInputData(1, poly2)
2473        bf.Update()
2474
2475        msh = Mesh(bf.GetOutput(), c=None)
2476        msh.flat()
2477
2478        msh.pipeline = OperationNode(
2479            "boolean " + operation,
2480            parents=[self, mesh2],
2481            shape="cylinder",
2482            comment=f"#pts {msh.dataset.GetNumberOfPoints()}",
2483        )
2484        msh.name = self.name + operation + mesh2.name
2485        return msh
2486
2487    def intersect_with(self, mesh2, tol=1e-06) -> Self:
2488        """
2489        Intersect this Mesh with the input surface to return a set of lines.
2490
2491        Examples:
2492            - [surf_intersect.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/surf_intersect.py)
2493
2494                ![](https://vedo.embl.es/images/basic/surfIntersect.png)
2495        """
2496        bf = vtki.new("IntersectionPolyDataFilter")
2497        bf.SetGlobalWarningDisplay(0)
2498        bf.SetTolerance(tol)
2499        bf.SetInputData(0, self.dataset)
2500        bf.SetInputData(1, mesh2.dataset)
2501        bf.Update()
2502        msh = Mesh(bf.GetOutput(), c="k", alpha=1).lighting("off")
2503        msh.properties.SetLineWidth(3)
2504        msh.pipeline = OperationNode(
2505            "intersect_with", parents=[self, mesh2], comment=f"#pts {msh.npoints}"
2506        )
2507        msh.name = "SurfaceIntersection"
2508        return msh
2509
2510    def intersect_with_line(self, p0, p1=None, return_ids=False, tol=0) -> Union[np.ndarray, Tuple[np.ndarray, np.ndarray]]:
2511        """
2512        Return the list of points intersecting the mesh
2513        along the segment defined by two points `p0` and `p1`.
2514
2515        Use `return_ids` to return the cell ids along with point coords
2516
2517        Example:
2518            ```python
2519            from vedo import *
2520            s = Spring()
2521            pts = s.intersect_with_line([0,0,0], [1,0.1,0])
2522            ln = Line([0,0,0], [1,0.1,0], c='blue')
2523            ps = Points(pts, r=10, c='r')
2524            show(s, ln, ps, bg='white').close()
2525            ```
2526            ![](https://user-images.githubusercontent.com/32848391/55967065-eee08300-5c79-11e9-8933-265e1bab9f7e.png)
2527        """
2528        if isinstance(p0, Points):
2529            p0, p1 = p0.vertices
2530
2531        if not self.line_locator:
2532            self.line_locator = vtki.new("OBBTree")
2533            self.line_locator.SetDataSet(self.dataset)
2534            if not tol:
2535                tol = mag(np.asarray(p1) - np.asarray(p0)) / 10000
2536            self.line_locator.SetTolerance(tol)
2537            self.line_locator.BuildLocator()
2538
2539        vpts = vtki.vtkPoints()
2540        idlist = vtki.vtkIdList()
2541        self.line_locator.IntersectWithLine(p0, p1, vpts, idlist)
2542        pts = []
2543        for i in range(vpts.GetNumberOfPoints()):
2544            intersection: MutableSequence[float] = [0, 0, 0]
2545            vpts.GetPoint(i, intersection)
2546            pts.append(intersection)
2547        pts2 = np.array(pts)
2548
2549        if return_ids:
2550            pts_ids = []
2551            for i in range(idlist.GetNumberOfIds()):
2552                cid = idlist.GetId(i)
2553                pts_ids.append(cid)
2554            return (pts2, np.array(pts_ids).astype(np.uint32))
2555
2556        return pts2
2557
2558    def intersect_with_plane(self, origin=(0, 0, 0), normal=(1, 0, 0)) -> Self:
2559        """
2560        Intersect this Mesh with a plane to return a set of lines.
2561
2562        Example:
2563            ```python
2564            from vedo import *
2565            sph = Sphere()
2566            mi = sph.clone().intersect_with_plane().join()
2567            print(mi.lines)
2568            show(sph, mi, axes=1).close()
2569            ```
2570            ![](https://vedo.embl.es/images/feats/intersect_plane.png)
2571        """
2572        plane = vtki.new("Plane")
2573        plane.SetOrigin(origin)
2574        plane.SetNormal(normal)
2575
2576        cutter = vtki.new("PolyDataPlaneCutter")
2577        cutter.SetInputData(self.dataset)
2578        cutter.SetPlane(plane)
2579        cutter.InterpolateAttributesOn()
2580        cutter.ComputeNormalsOff()
2581        cutter.Update()
2582
2583        msh = Mesh(cutter.GetOutput())
2584        msh.c('k').lw(3).lighting("off")
2585        msh.pipeline = OperationNode(
2586            "intersect_with_plan",
2587            parents=[self],
2588            comment=f"#pts {msh.dataset.GetNumberOfPoints()}",
2589        )
2590        msh.name = "PlaneIntersection"
2591        return msh
2592    
2593    def cut_closed_surface(self, origins, normals, invert=False, return_assembly=False) -> Union[Self, "vedo.Assembly"]:
2594        """
2595        Cut/clip a closed surface mesh with a collection of planes.
2596        This will produce a new closed surface by creating new polygonal
2597        faces where the input surface hits the planes.
2598
2599        The orientation of the polygons that form the surface is important.
2600        Polygons have a front face and a back face, and it's the back face that defines
2601        the interior or "solid" region of the closed surface.
2602        When a plane cuts through a "solid" region, a new cut face is generated,
2603        but not when a clipping plane cuts through a hole or "empty" region.
2604        This distinction is crucial when dealing with complex surfaces.
2605        Note that if a simple surface has its back faces pointing outwards,
2606        then that surface defines a hole in a potentially infinite solid.
2607
2608        Non-manifold surfaces should not be used with this method. 
2609
2610        Arguments:
2611            origins : (list)
2612                list of plane origins
2613            normals : (list)
2614                list of plane normals
2615            invert : (bool)
2616                invert the clipping.
2617            return_assembly : (bool)
2618                return the cap and the clipped surfaces as a `vedo.Assembly`.
2619        
2620        Example:
2621            ```python
2622            from vedo import *
2623            s = Sphere(res=50).linewidth(1)
2624            origins = [[-0.7, 0, 0], [0, -0.6, 0]]
2625            normals = [[-1, 0, 0], [0, -1, 0]]
2626            s.cut_closed_surface(origins, normals)
2627            show(s, axes=1).close()
2628            ```
2629        """        
2630        planes = vtki.new("PlaneCollection")
2631        for p, s in zip(origins, normals):
2632            plane = vtki.vtkPlane()
2633            plane.SetOrigin(vedo.utils.make3d(p))
2634            plane.SetNormal(vedo.utils.make3d(s))
2635            planes.AddItem(plane)
2636        clipper = vtki.new("ClipClosedSurface")
2637        clipper.SetInputData(self.dataset)
2638        clipper.SetClippingPlanes(planes)
2639        clipper.PassPointDataOn()
2640        clipper.GenerateFacesOn()
2641        clipper.SetScalarModeToLabels()
2642        clipper.TriangulationErrorDisplayOn()
2643        clipper.SetInsideOut(not invert)
2644
2645        if return_assembly:
2646            clipper.GenerateClipFaceOutputOn()
2647            clipper.Update()
2648            parts = []
2649            for i in range(clipper.GetNumberOfOutputPorts()):
2650                msh = Mesh(clipper.GetOutput(i))
2651                msh.copy_properties_from(self)
2652                msh.name = "CutClosedSurface"
2653                msh.pipeline = OperationNode(
2654                    "cut_closed_surface",
2655                    parents=[self],
2656                    comment=f"#pts {msh.dataset.GetNumberOfPoints()}",
2657                )
2658                parts.append(msh)
2659            asse = vedo.Assembly(parts)
2660            asse.name = "CutClosedSurface"
2661            return asse
2662
2663        else:
2664            clipper.GenerateClipFaceOutputOff()
2665            clipper.Update()
2666            self._update(clipper.GetOutput())
2667            self.flat()
2668            self.name = "CutClosedSurface"
2669            self.pipeline = OperationNode(
2670                "cut_closed_surface",
2671                parents=[self],
2672                comment=f"#pts {self.dataset.GetNumberOfPoints()}",
2673            )
2674            return self
2675
2676    def collide_with(self, mesh2, tol=0, return_bool=False) -> Union[Self, bool]:
2677        """
2678        Collide this Mesh with the input surface.
2679        Information is stored in `ContactCells1` and `ContactCells2`.
2680        """
2681        ipdf = vtki.new("CollisionDetectionFilter")
2682        # ipdf.SetGlobalWarningDisplay(0)
2683
2684        transform0 = vtki.vtkTransform()
2685        transform1 = vtki.vtkTransform()
2686
2687        # ipdf.SetBoxTolerance(tol)
2688        ipdf.SetCellTolerance(tol)
2689        ipdf.SetInputData(0, self.dataset)
2690        ipdf.SetInputData(1, mesh2.dataset)
2691        ipdf.SetTransform(0, transform0)
2692        ipdf.SetTransform(1, transform1)
2693        if return_bool:
2694            ipdf.SetCollisionModeToFirstContact()
2695        else:
2696            ipdf.SetCollisionModeToAllContacts()
2697        ipdf.Update()
2698
2699        if return_bool:
2700            return bool(ipdf.GetNumberOfContacts())
2701
2702        msh = Mesh(ipdf.GetContactsOutput(), "k", 1).lighting("off")
2703        msh.metadata["ContactCells1"] = vtk2numpy(
2704            ipdf.GetOutput(0).GetFieldData().GetArray("ContactCells")
2705        )
2706        msh.metadata["ContactCells2"] = vtk2numpy(
2707            ipdf.GetOutput(1).GetFieldData().GetArray("ContactCells")
2708        )
2709        msh.properties.SetLineWidth(3)
2710
2711        msh.pipeline = OperationNode(
2712            "collide_with",
2713            parents=[self, mesh2],
2714            comment=f"#pts {msh.dataset.GetNumberOfPoints()}",
2715        )
2716        msh.name = "SurfaceCollision"
2717        return msh
2718
2719    def geodesic(self, start, end) -> Self:
2720        """
2721        Dijkstra algorithm to compute the geodesic line.
2722        Takes as input a polygonal mesh and performs a single source shortest path calculation.
2723
2724        The output mesh contains the array "VertexIDs" that contains the ordered list of vertices
2725        traversed to get from the start vertex to the end vertex.
2726        
2727        Arguments:
2728            start : (int, list)
2729                start vertex index or close point `[x,y,z]`
2730            end :  (int, list)
2731                end vertex index or close point `[x,y,z]`
2732
2733        Examples:
2734            - [geodesic_curve.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/geodesic_curve.py)
2735
2736                ![](https://vedo.embl.es/images/advanced/geodesic.png)
2737        """
2738        if is_sequence(start):
2739            cc = self.vertices
2740            pa = Points(cc)
2741            start = pa.closest_point(start, return_point_id=True)
2742            end = pa.closest_point(end, return_point_id=True)
2743
2744        dijkstra = vtki.new("DijkstraGraphGeodesicPath")
2745        dijkstra.SetInputData(self.dataset)
2746        dijkstra.SetStartVertex(end)  # inverted in vtk
2747        dijkstra.SetEndVertex(start)
2748        dijkstra.Update()
2749
2750        weights = vtki.vtkDoubleArray()
2751        dijkstra.GetCumulativeWeights(weights)
2752
2753        idlist = dijkstra.GetIdList()
2754        ids = [idlist.GetId(i) for i in range(idlist.GetNumberOfIds())]
2755
2756        length = weights.GetMaxId() + 1
2757        arr = np.zeros(length)
2758        for i in range(length):
2759            arr[i] = weights.GetTuple(i)[0]
2760
2761        poly = dijkstra.GetOutput()
2762
2763        vdata = numpy2vtk(arr)
2764        vdata.SetName("CumulativeWeights")
2765        poly.GetPointData().AddArray(vdata)
2766
2767        vdata2 = numpy2vtk(ids, dtype=np.uint)
2768        vdata2.SetName("VertexIDs")
2769        poly.GetPointData().AddArray(vdata2)
2770        poly.GetPointData().Modified()
2771
2772        dmesh = Mesh(poly).copy_properties_from(self)
2773        dmesh.lw(3).alpha(1).lighting("off")
2774        dmesh.name = "GeodesicLine"
2775
2776        dmesh.pipeline = OperationNode(
2777            "GeodesicLine",
2778            parents=[self],
2779            comment=f"#steps {poly.GetNumberOfPoints()}",
2780        )
2781        return dmesh
2782
2783    #####################################################################
2784    ### Stuff returning a Volume object
2785    #####################################################################
2786    def binarize(
2787        self,
2788        values=(255, 0),
2789        spacing=None,
2790        dims=None,
2791        origin=None,
2792    ) -> "vedo.Volume":
2793        """
2794        Convert a `Mesh` into a `Volume` where
2795        the interior voxels value is set to `values[0]` (255 by default), while
2796        the exterior voxels value is set to `values[1]` (0 by default).
2797
2798        Arguments:
2799            values : (list)
2800                background and foreground values.
2801            spacing : (list)
2802                voxel spacing in x, y and z.
2803            dims : (list)
2804                dimensions (nr. of voxels) of the output volume.
2805            origin : (list)
2806                position in space of the (0,0,0) voxel.
2807
2808        Examples:
2809            - [mesh2volume.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/mesh2volume.py)
2810
2811                ![](https://vedo.embl.es/images/volumetric/mesh2volume.png)
2812        """
2813        assert len(values) == 2, "values must be a list of 2 values"
2814        fg_value, bg_value = values
2815
2816        bounds = self.bounds()
2817        if spacing is None:  # compute spacing
2818            spacing = [0, 0, 0]
2819            diagonal = np.sqrt(
2820                  (bounds[1] - bounds[0]) ** 2
2821                + (bounds[3] - bounds[2]) ** 2
2822                + (bounds[5] - bounds[4]) ** 2
2823            )
2824            spacing[0] = spacing[1] = spacing[2] = diagonal / 250.0
2825
2826        if dims is None:  # compute dimensions
2827            dim = [0, 0, 0]
2828            for i in [0, 1, 2]:
2829                dim[i] = int(np.ceil((bounds[i*2+1] - bounds[i*2]) / spacing[i]))
2830        else:
2831            dim = dims
2832        
2833        white_img = vtki.vtkImageData()
2834        white_img.SetDimensions(dim)
2835        white_img.SetSpacing(spacing)
2836        white_img.SetExtent(0, dim[0]-1, 0, dim[1]-1, 0, dim[2]-1)
2837
2838        if origin is None:
2839            origin = [0, 0, 0]
2840            origin[0] = bounds[0] + spacing[0]
2841            origin[1] = bounds[2] + spacing[1]
2842            origin[2] = bounds[4] + spacing[2]
2843        white_img.SetOrigin(origin)
2844
2845        # if direction_matrix is not None:
2846        #     white_img.SetDirectionMatrix(direction_matrix)
2847
2848        white_img.AllocateScalars(vtki.VTK_UNSIGNED_CHAR, 1)
2849
2850        # fill the image with foreground voxels:
2851        white_img.GetPointData().GetScalars().Fill(fg_value)
2852
2853        # polygonal data --> image stencil:
2854        pol2stenc = vtki.new("PolyDataToImageStencil")
2855        pol2stenc.SetInputData(self.dataset)
2856        pol2stenc.SetOutputOrigin(white_img.GetOrigin())
2857        pol2stenc.SetOutputSpacing(white_img.GetSpacing())
2858        pol2stenc.SetOutputWholeExtent(white_img.GetExtent())
2859        pol2stenc.Update()
2860
2861        # cut the corresponding white image and set the background:
2862        imgstenc = vtki.new("ImageStencil")
2863        imgstenc.SetInputData(white_img)
2864        imgstenc.SetStencilConnection(pol2stenc.GetOutputPort())
2865        # imgstenc.SetReverseStencil(True)
2866        imgstenc.SetBackgroundValue(bg_value)
2867        imgstenc.Update()
2868
2869        vol = vedo.Volume(imgstenc.GetOutput())
2870        vol.name = "BinarizedVolume"
2871        vol.pipeline = OperationNode(
2872            "binarize",
2873            parents=[self],
2874            comment=f"dims={tuple(vol.dimensions())}",
2875            c="#e9c46a:#0096c7",
2876        )
2877        return vol
2878
2879    def signed_distance(self, bounds=None, dims=(20, 20, 20), invert=False, maxradius=None) -> "vedo.Volume":
2880        """
2881        Compute the `Volume` object whose voxels contains 
2882        the signed distance from the mesh.
2883
2884        Arguments:
2885            bounds : (list)
2886                bounds of the output volume
2887            dims : (list)
2888                dimensions (nr. of voxels) of the output volume
2889            invert : (bool)
2890                flip the sign
2891
2892        Examples:
2893            - [volume_from_mesh.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/volume_from_mesh.py)
2894        """
2895        if maxradius is not None:
2896            vedo.logger.warning(
2897                "in signedDistance(maxradius=...) is ignored. (Only valid for pointclouds)."
2898            )
2899        if bounds is None:
2900            bounds = self.bounds()
2901        sx = (bounds[1] - bounds[0]) / dims[0]
2902        sy = (bounds[3] - bounds[2]) / dims[1]
2903        sz = (bounds[5] - bounds[4]) / dims[2]
2904
2905        img = vtki.vtkImageData()
2906        img.SetDimensions(dims)
2907        img.SetSpacing(sx, sy, sz)
2908        img.SetOrigin(bounds[0], bounds[2], bounds[4])
2909        img.AllocateScalars(vtki.VTK_FLOAT, 1)
2910
2911        imp = vtki.new("ImplicitPolyDataDistance")
2912        imp.SetInput(self.dataset)
2913        b2 = bounds[2]
2914        b4 = bounds[4]
2915        d0, d1, d2 = dims
2916
2917        for i in range(d0):
2918            x = i * sx + bounds[0]
2919            for j in range(d1):
2920                y = j * sy + b2
2921                for k in range(d2):
2922                    v = imp.EvaluateFunction((x, y, k * sz + b4))
2923                    if invert:
2924                        v = -v
2925                    img.SetScalarComponentFromFloat(i, j, k, 0, v)
2926
2927        vol = vedo.Volume(img)
2928        vol.name = "SignedVolume"
2929
2930        vol.pipeline = OperationNode(
2931            "signed_distance",
2932            parents=[self],
2933            comment=f"dims={tuple(vol.dimensions())}",
2934            c="#e9c46a:#0096c7",
2935        )
2936        return vol
2937
2938    def tetralize(
2939        self,
2940        side=0.02,
2941        nmax=300_000,
2942        gap=None,
2943        subsample=False,
2944        uniform=True,
2945        seed=0,
2946        debug=False,
2947    ) -> "vedo.TetMesh":
2948        """
2949        Tetralize a closed polygonal mesh. Return a `TetMesh`.
2950
2951        Arguments:
2952            side : (float)
2953                desired side of the single tetras as fraction of the bounding box diagonal.
2954                Typical values are in the range (0.01 - 0.03)
2955            nmax : (int)
2956                maximum random numbers to be sampled in the bounding box
2957            gap : (float)
2958                keep this minimum distance from the surface,
2959                if None an automatic choice is made.
2960            subsample : (bool)
2961                subsample input surface, the geometry might be affected
2962                (the number of original faces reduceed), but higher tet quality might be obtained.
2963            uniform : (bool)
2964                generate tets more uniformly packed in the interior of the mesh
2965            seed : (int)
2966                random number generator seed
2967            debug : (bool)
2968                show an intermediate plot with sampled points
2969
2970        Examples:
2971            - [tetralize_surface.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/tetralize_surface.py)
2972
2973                ![](https://vedo.embl.es/images/volumetric/tetralize_surface.jpg)
2974        """
2975        surf = self.clone().clean().compute_normals()
2976        d = surf.diagonal_size()
2977        if gap is None:
2978            gap = side * d * np.sqrt(2 / 3)
2979        n = int(min((1 / side) ** 3, nmax))
2980
2981        # fill the space w/ points
2982        x0, x1, y0, y1, z0, z1 = surf.bounds()
2983
2984        if uniform:
2985            pts = vedo.utils.pack_spheres([x0, x1, y0, y1, z0, z1], side * d * 1.42)
2986            pts += np.random.randn(len(pts), 3) * side * d * 1.42 / 100  # some small jitter
2987        else:
2988            disp = np.array([x0 + x1, y0 + y1, z0 + z1]) / 2
2989            np.random.seed(seed)
2990            pts = (np.random.rand(n, 3) - 0.5) * np.array([x1 - x0, y1 - y0, z1 - z0]) + disp
2991
2992        normals = surf.celldata["Normals"]
2993        cc = surf.cell_centers
2994        subpts = cc - normals * gap * 1.05
2995        pts = pts.tolist() + subpts.tolist()
2996
2997        if debug:
2998            print(".. tetralize(): subsampling and cleaning")
2999
3000        fillpts = surf.inside_points(pts)
3001        fillpts.subsample(side)
3002
3003        if gap:
3004            fillpts.distance_to(surf)
3005            fillpts.threshold("Distance", above=gap)
3006
3007        if subsample:
3008            surf.subsample(side)
3009
3010        merged_fs = vedo.merge(fillpts, surf)
3011        tmesh = merged_fs.generate_delaunay3d()
3012        tcenters = tmesh.cell_centers
3013
3014        ids = surf.inside_points(tcenters, return_ids=True)
3015        ins = np.zeros(tmesh.ncells)
3016        ins[ids] = 1
3017
3018        if debug:
3019            # vedo.pyplot.histogram(fillpts.pointdata["Distance"], xtitle=f"gap={gap}").show().close()
3020            edges = self.edges
3021            points = self.vertices
3022            elen = mag(points[edges][:, 0, :] - points[edges][:, 1, :])
3023            histo = vedo.pyplot.histogram(elen, xtitle="edge length", xlim=(0, 3 * side * d))
3024            print(".. edges min, max", elen.min(), elen.max())
3025            fillpts.cmap("bone")
3026            vedo.show(
3027                [
3028                    [
3029                        f"This is a debug plot.\n\nGenerated points: {n}\ngap: {gap}",
3030                        surf.wireframe().alpha(0.2),
3031                        vedo.addons.Axes(surf),
3032                        fillpts,
3033                        Points(subpts).c("r4").ps(3),
3034                    ],
3035                    [f"Edges mean length: {np.mean(elen)}\n\nPress q to continue", histo],
3036                ],
3037                N=2,
3038                sharecam=False,
3039                new=True,
3040            ).close()
3041            print(".. thresholding")
3042
3043        tmesh.celldata["inside"] = ins.astype(np.uint8)
3044        tmesh.threshold("inside", above=0.9)
3045        tmesh.celldata.remove("inside")
3046
3047        if debug:
3048            print(f".. tetralize() completed, ntets = {tmesh.ncells}")
3049
3050        tmesh.pipeline = OperationNode(
3051            "tetralize",
3052            parents=[self],
3053            comment=f"#tets = {tmesh.ncells}",
3054            c="#e9c46a:#9e2a2b",
3055        )
3056        return tmesh
  29class Mesh(MeshVisual, Points):
  30    """
  31    Build an instance of object `Mesh` derived from `vedo.PointCloud`.
  32    """
  33
  34    def __init__(self, inputobj=None, c="gold", alpha=1):
  35        """
  36        Initialize a ``Mesh`` object.
  37
  38        Arguments:
  39            inputobj : (str, vtkPolyData, vtkActor, vedo.Mesh)
  40                If inputobj is `None` an empty mesh is created.
  41                If inputobj is a `str` then it is interpreted as the name of a file to load as mesh.
  42                If inputobj is an `vtkPolyData` or `vtkActor` or `vedo.Mesh`
  43                then a shallow copy of it is created.
  44                If inputobj is a `vedo.Mesh` then a shallow copy of it is created.
  45
  46        Examples:
  47            - [buildmesh.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/buildmesh.py)
  48            (and many others!)
  49
  50            ![](https://vedo.embl.es/images/basic/buildmesh.png)
  51        """
  52        # print("INIT MESH", super())
  53        super().__init__()
  54
  55        self.name = "Mesh"
  56
  57        if inputobj is None:
  58            # self.dataset = vtki.vtkPolyData()
  59            pass
  60
  61        elif isinstance(inputobj, str):
  62            self.dataset = vedo.file_io.load(inputobj).dataset
  63            self.filename = inputobj
  64
  65        elif isinstance(inputobj, vtki.vtkPolyData):
  66            # self.dataset.DeepCopy(inputobj) # NO
  67            self.dataset = inputobj
  68            if self.dataset.GetNumberOfCells() == 0:
  69                carr = vtki.vtkCellArray()
  70                for i in range(inputobj.GetNumberOfPoints()):
  71                    carr.InsertNextCell(1)
  72                    carr.InsertCellPoint(i)
  73                self.dataset.SetVerts(carr)
  74
  75        elif isinstance(inputobj, Mesh):
  76            self.dataset = inputobj.dataset
  77
  78        elif is_sequence(inputobj):
  79            ninp = len(inputobj)
  80            if   ninp == 4:  # assume input is [vertices, faces, lines, strips]
  81                self.dataset = buildPolyData(inputobj[0], inputobj[1], inputobj[2], inputobj[3])
  82            elif ninp == 3:  # assume input is [vertices, faces, lines]
  83                self.dataset = buildPolyData(inputobj[0], inputobj[1], inputobj[2])
  84            elif ninp == 2:  # assume input is [vertices, faces]
  85                self.dataset = buildPolyData(inputobj[0], inputobj[1])
  86            elif ninp == 1:  # assume input is [vertices]
  87                self.dataset = buildPolyData(inputobj[0])
  88            else:
  89                vedo.logger.error("input must be a list of max 4 elements.")
  90                raise ValueError()
  91
  92        elif isinstance(inputobj, vtki.vtkActor):
  93            self.dataset.DeepCopy(inputobj.GetMapper().GetInput())
  94            v = inputobj.GetMapper().GetScalarVisibility()
  95            self.mapper.SetScalarVisibility(v)
  96            pr = vtki.vtkProperty()
  97            pr.DeepCopy(inputobj.GetProperty())
  98            self.actor.SetProperty(pr)
  99            self.properties = pr
 100
 101        elif isinstance(inputobj, (vtki.vtkStructuredGrid, vtki.vtkRectilinearGrid)):
 102            gf = vtki.new("GeometryFilter")
 103            gf.SetInputData(inputobj)
 104            gf.Update()
 105            self.dataset = gf.GetOutput()
 106
 107        elif "meshlab" in str(type(inputobj)):
 108            self.dataset = vedo.utils.meshlab2vedo(inputobj).dataset
 109
 110        elif "meshlib" in str(type(inputobj)):
 111            import meshlib.mrmeshnumpy as mrmeshnumpy
 112            self.dataset = buildPolyData(
 113                mrmeshnumpy.getNumpyVerts(inputobj),
 114                mrmeshnumpy.getNumpyFaces(inputobj.topology),
 115            )
 116
 117        elif "trimesh" in str(type(inputobj)):
 118            self.dataset = vedo.utils.trimesh2vedo(inputobj).dataset
 119
 120        elif "meshio" in str(type(inputobj)):
 121            # self.dataset = vedo.utils.meshio2vedo(inputobj) ##TODO
 122            if len(inputobj.cells) > 0:
 123                mcells = []
 124                for cellblock in inputobj.cells:
 125                    if cellblock.type in ("triangle", "quad"):
 126                        mcells += cellblock.data.tolist()
 127                self.dataset = buildPolyData(inputobj.points, mcells)
 128            else:
 129                self.dataset = buildPolyData(inputobj.points, None)
 130            # add arrays:
 131            try:
 132                if len(inputobj.point_data) > 0:
 133                    for k in inputobj.point_data.keys():
 134                        vdata = numpy2vtk(inputobj.point_data[k])
 135                        vdata.SetName(str(k))
 136                        self.dataset.GetPointData().AddArray(vdata)
 137            except AssertionError:
 138                print("Could not add meshio point data, skip.")
 139
 140        else:
 141            try:
 142                gf = vtki.new("GeometryFilter")
 143                gf.SetInputData(inputobj)
 144                gf.Update()
 145                self.dataset = gf.GetOutput()
 146            except:
 147                vedo.logger.error(f"cannot build mesh from type {type(inputobj)}")
 148                raise RuntimeError()
 149
 150        self.mapper.SetInputData(self.dataset)
 151        self.actor.SetMapper(self.mapper)
 152
 153        self.properties.SetInterpolationToPhong()
 154        self.properties.SetColor(get_color(c))
 155
 156        if alpha is not None:
 157            self.properties.SetOpacity(alpha)
 158
 159        self.mapper.SetInterpolateScalarsBeforeMapping(
 160            vedo.settings.interpolate_scalars_before_mapping
 161        )
 162
 163        if vedo.settings.use_polygon_offset:
 164            self.mapper.SetResolveCoincidentTopologyToPolygonOffset()
 165            pof = vedo.settings.polygon_offset_factor
 166            pou = vedo.settings.polygon_offset_units
 167            self.mapper.SetResolveCoincidentTopologyPolygonOffsetParameters(pof, pou)
 168
 169        n = self.dataset.GetNumberOfPoints()
 170        self.pipeline = OperationNode(self, comment=f"#pts {n}")
 171
 172    def _repr_html_(self):
 173        """
 174        HTML representation of the Mesh object for Jupyter Notebooks.
 175
 176        Returns:
 177            HTML text with the image and some properties.
 178        """
 179        import io
 180        import base64
 181        from PIL import Image
 182
 183        library_name = "vedo.mesh.Mesh"
 184        help_url = "https://vedo.embl.es/docs/vedo/mesh.html#Mesh"
 185
 186        arr = self.thumbnail()
 187        im = Image.fromarray(arr)
 188        buffered = io.BytesIO()
 189        im.save(buffered, format="PNG", quality=100)
 190        encoded = base64.b64encode(buffered.getvalue()).decode("utf-8")
 191        url = "data:image/png;base64," + encoded
 192        image = f"<img src='{url}'></img>"
 193
 194        bounds = "<br/>".join(
 195            [
 196                precision(min_x, 4) + " ... " + precision(max_x, 4)
 197                for min_x, max_x in zip(self.bounds()[::2], self.bounds()[1::2])
 198            ]
 199        )
 200        average_size = "{size:.3f}".format(size=self.average_size())
 201
 202        help_text = ""
 203        if self.name:
 204            help_text += f"<b> {self.name}: &nbsp&nbsp</b>"
 205        help_text += '<b><a href="' + help_url + '" target="_blank">' + library_name + "</a></b>"
 206        if self.filename:
 207            dots = ""
 208            if len(self.filename) > 30:
 209                dots = "..."
 210            help_text += f"<br/><code><i>({dots}{self.filename[-30:]})</i></code>"
 211
 212        pdata = ""
 213        if self.dataset.GetPointData().GetScalars():
 214            if self.dataset.GetPointData().GetScalars().GetName():
 215                name = self.dataset.GetPointData().GetScalars().GetName()
 216                pdata = "<tr><td><b> point data array </b></td><td>" + name + "</td></tr>"
 217
 218        cdata = ""
 219        if self.dataset.GetCellData().GetScalars():
 220            if self.dataset.GetCellData().GetScalars().GetName():
 221                name = self.dataset.GetCellData().GetScalars().GetName()
 222                cdata = "<tr><td><b> cell data array </b></td><td>" + name + "</td></tr>"
 223
 224        allt = [
 225            "<table>",
 226            "<tr>",
 227            "<td>",
 228            image,
 229            "</td>",
 230            "<td style='text-align: center; vertical-align: center;'><br/>",
 231            help_text,
 232            "<table>",
 233            "<tr><td><b> bounds </b> <br/> (x/y/z) </td><td>" + str(bounds) + "</td></tr>",
 234            "<tr><td><b> center of mass </b></td><td>"
 235            + precision(self.center_of_mass(), 3)
 236            + "</td></tr>",
 237            "<tr><td><b> average size </b></td><td>" + str(average_size) + "</td></tr>",
 238            "<tr><td><b> nr. points&nbsp/&nbspfaces </b></td><td>"
 239            + str(self.npoints)
 240            + "&nbsp/&nbsp"
 241            + str(self.ncells)
 242            + "</td></tr>",
 243            pdata,
 244            cdata,
 245            "</table>",
 246            "</table>",
 247        ]
 248        return "\n".join(allt)
 249
 250    def faces(self, ids=()):
 251        """DEPRECATED. Use property `mesh.cells` instead."""
 252        vedo.printc("WARNING: use property mesh.cells instead of mesh.faces()",c='y')
 253        return self.cells
 254    
 255    @property
 256    def edges(self):
 257        """Return an array containing the edges connectivity."""
 258        extractEdges = vtki.new("ExtractEdges")
 259        extractEdges.SetInputData(self.dataset)
 260        # eed.UseAllPointsOn()
 261        extractEdges.Update()
 262        lpoly = extractEdges.GetOutput()
 263
 264        arr1d = vtk2numpy(lpoly.GetLines().GetData())
 265        # [nids1, id0 ... idn, niids2, id0 ... idm,  etc].
 266
 267        i = 0
 268        conn = []
 269        n = len(arr1d)
 270        for _ in range(n):
 271            cell = [arr1d[i + k + 1] for k in range(arr1d[i])]
 272            conn.append(cell)
 273            i += arr1d[i] + 1
 274            if i >= n:
 275                break
 276        return conn  # cannot always make a numpy array of it!
 277
 278    @property
 279    def cell_normals(self):
 280        """
 281        Retrieve face normals as a numpy array.
 282        Check out also `compute_normals(cells=True)` and `compute_normals_with_pca()`.
 283        """
 284        vtknormals = self.dataset.GetCellData().GetNormals()
 285        numpy_normals = vtk2numpy(vtknormals)
 286        if len(numpy_normals) == 0 and len(self.cells) != 0:
 287            raise ValueError("VTK failed to return any normal vectors. You may need to call `Mesh.compute_normals()` before accessing `Mesh.cell_normals`.")
 288        return numpy_normals
 289
 290    def compute_normals(self, points=True, cells=True, feature_angle=None, consistency=True) -> Self:
 291        """
 292        Compute cell and vertex normals for the mesh.
 293
 294        Arguments:
 295            points : (bool)
 296                do the computation for the vertices too
 297            cells : (bool)
 298                do the computation for the cells too
 299            feature_angle : (float)
 300                specify the angle that defines a sharp edge.
 301                If the difference in angle across neighboring polygons is greater than this value,
 302                the shared edge is considered "sharp" and it is split.
 303            consistency : (bool)
 304                turn on/off the enforcement of consistent polygon ordering.
 305
 306        .. warning::
 307            If `feature_angle` is set then the Mesh can be modified, and it
 308            can have a different nr. of vertices from the original.
 309
 310            Note that the appearance of the mesh may change if the normals are computed,
 311            as shading is automatically enabled when such information is present.
 312            Use `mesh.flat()` to avoid smoothing effects.
 313        """
 314        pdnorm = vtki.new("PolyDataNormals")
 315        pdnorm.SetInputData(self.dataset)
 316        pdnorm.SetComputePointNormals(points)
 317        pdnorm.SetComputeCellNormals(cells)
 318        pdnorm.SetConsistency(consistency)
 319        pdnorm.FlipNormalsOff()
 320        if feature_angle:
 321            pdnorm.SetSplitting(True)
 322            pdnorm.SetFeatureAngle(feature_angle)
 323        else:
 324            pdnorm.SetSplitting(False)
 325        pdnorm.Update()
 326        out = pdnorm.GetOutput()
 327        self._update(out, reset_locators=False)
 328        return self
 329
 330    def reverse(self, cells=True, normals=False) -> Self:
 331        """
 332        Reverse the order of polygonal cells
 333        and/or reverse the direction of point and cell normals.
 334
 335        Two flags are used to control these operations:
 336            - `cells=True` reverses the order of the indices in the cell connectivity list.
 337                If cell is a list of IDs only those cells will be reversed.
 338            - `normals=True` reverses the normals by multiplying the normal vector by -1
 339                (both point and cell normals, if present).
 340        """
 341        poly = self.dataset
 342
 343        if is_sequence(cells):
 344            for cell in cells:
 345                poly.ReverseCell(cell)
 346            poly.GetCellData().Modified()
 347            return self  ##############
 348
 349        rev = vtki.new("ReverseSense")
 350        if cells:
 351            rev.ReverseCellsOn()
 352        else:
 353            rev.ReverseCellsOff()
 354        if normals:
 355            rev.ReverseNormalsOn()
 356        else:
 357            rev.ReverseNormalsOff()
 358        rev.SetInputData(poly)
 359        rev.Update()
 360        self._update(rev.GetOutput(), reset_locators=False)
 361        self.pipeline = OperationNode("reverse", parents=[self])
 362        return self
 363
 364    def volume(self) -> float:
 365        """
 366        Compute the volume occupied by mesh.
 367        The mesh must be triangular for this to work.
 368        To triangulate a mesh use `mesh.triangulate()`.
 369        """
 370        mass = vtki.new("MassProperties")
 371        mass.SetGlobalWarningDisplay(0)
 372        mass.SetInputData(self.dataset)
 373        mass.Update()
 374        mass.SetGlobalWarningDisplay(1)
 375        return mass.GetVolume()
 376
 377    def area(self) -> float:
 378        """
 379        Compute the surface area of the mesh.
 380        The mesh must be triangular for this to work.
 381        To triangulate a mesh use `mesh.triangulate()`.
 382        """
 383        mass = vtki.new("MassProperties")
 384        mass.SetGlobalWarningDisplay(0)
 385        mass.SetInputData(self.dataset)
 386        mass.Update()
 387        mass.SetGlobalWarningDisplay(1)
 388        return mass.GetSurfaceArea()
 389
 390    def is_closed(self) -> bool:
 391        """
 392        Return `True` if the mesh is watertight.
 393        Note that if the mesh contains coincident points the result may be flase.
 394        Use in this case `mesh.clean()` to merge coincident points.
 395        """
 396        fe = vtki.new("FeatureEdges")
 397        fe.BoundaryEdgesOn()
 398        fe.FeatureEdgesOff()
 399        fe.NonManifoldEdgesOn()
 400        fe.SetInputData(self.dataset)
 401        fe.Update()
 402        ne = fe.GetOutput().GetNumberOfCells()
 403        return not bool(ne)
 404
 405    def is_manifold(self) -> bool:
 406        """Return `True` if the mesh is manifold."""
 407        fe = vtki.new("FeatureEdges")
 408        fe.BoundaryEdgesOff()
 409        fe.FeatureEdgesOff()
 410        fe.NonManifoldEdgesOn()
 411        fe.SetInputData(self.dataset)
 412        fe.Update()
 413        ne = fe.GetOutput().GetNumberOfCells()
 414        return not bool(ne)
 415
 416    def non_manifold_faces(self, remove=True, tol="auto") -> Self:
 417        """
 418        Detect and (try to) remove non-manifold faces of a triangular mesh:
 419
 420            - set `remove` to `False` to mark cells without removing them.
 421            - set `tol=0` for zero-tolerance, the result will be manifold but with holes.
 422            - set `tol>0` to cut off non-manifold faces, and try to recover the good ones.
 423            - set `tol="auto"` to make an automatic choice of the tolerance.
 424        """
 425        # mark original point and cell ids
 426        self.add_ids()
 427        toremove = self.boundaries(
 428            boundary_edges=False,
 429            non_manifold_edges=True,
 430            cell_edge=True,
 431            return_cell_ids=True,
 432        )
 433        if len(toremove) == 0: # type: ignore
 434            return self
 435
 436        points = self.vertices
 437        faces = self.cells
 438        centers = self.cell_centers
 439
 440        copy = self.clone()
 441        copy.delete_cells(toremove).clean()
 442        copy.compute_normals(cells=False)
 443        normals = copy.vertex_normals
 444        deltas, deltas_i = [], []
 445
 446        for i in vedo.utils.progressbar(toremove, delay=3, title="recover faces"):
 447            pids = copy.closest_point(centers[i], n=3, return_point_id=True)
 448            norms = normals[pids]
 449            n = np.mean(norms, axis=0)
 450            dn = np.linalg.norm(n)
 451            if not dn:
 452                continue
 453            n = n / dn
 454
 455            p0, p1, p2 = points[faces[i]][:3]
 456            v = np.cross(p1 - p0, p2 - p0)
 457            lv = np.linalg.norm(v)
 458            if not lv:
 459                continue
 460            v = v / lv
 461
 462            cosa = 1 - np.dot(n, v)
 463            deltas.append(cosa)
 464            deltas_i.append(i)
 465
 466        recover = []
 467        if len(deltas) > 0:
 468            mean_delta = np.mean(deltas)
 469            err_delta = np.std(deltas)
 470            txt = ""
 471            if tol == "auto":  # automatic choice
 472                tol = mean_delta / 5
 473                txt = f"\n Automatic tol. : {tol: .4f}"
 474            for i, cosa in zip(deltas_i, deltas):
 475                if cosa < tol:
 476                    recover.append(i)
 477
 478            vedo.logger.info(
 479                f"\n --------- Non manifold faces ---------"
 480                f"\n Average tol.   : {mean_delta: .4f} +- {err_delta: .4f}{txt}"
 481                f"\n Removed faces  : {len(toremove)}" # type: ignore
 482                f"\n Recovered faces: {len(recover)}"
 483            )
 484
 485        toremove = list(set(toremove) - set(recover)) # type: ignore
 486
 487        if not remove:
 488            mark = np.zeros(self.ncells, dtype=np.uint8)
 489            mark[recover] = 1
 490            mark[toremove] = 2
 491            self.celldata["NonManifoldCell"] = mark
 492        else:
 493            self.delete_cells(toremove) # type: ignore
 494
 495        self.pipeline = OperationNode(
 496            "non_manifold_faces",
 497            parents=[self],
 498            comment=f"#cells {self.dataset.GetNumberOfCells()}",
 499        )
 500        return self
 501
 502
 503    def euler_characteristic(self) -> int:
 504        """
 505        Compute the Euler characteristic of the mesh.
 506        The Euler characteristic is a topological invariant for surfaces.
 507        """
 508        return self.npoints - len(self.edges) + self.ncells
 509
 510    def genus(self) -> int:
 511        """
 512        Compute the genus of the mesh.
 513        The genus is a topological invariant for surfaces.
 514        """
 515        nb = len(self.boundaries().split()) - 1
 516        return (2 - self.euler_characteristic() - nb ) / 2
 517    
 518    def to_reeb_graph(self, field_id=0):
 519        """
 520        Convert the mesh into a Reeb graph.
 521        The Reeb graph is a topological structure that captures the evolution
 522        of the level sets of a scalar field.
 523
 524        Arguments:
 525            field_id : (int)
 526                the id of the scalar field to use.
 527        
 528        Example:
 529            ```python
 530            from vedo import *
 531            mesh = Mesh("https://discourse.paraview.org/uploads/short-url/qVuZ1fiRjwhE1qYtgGE2HGXybgo.stl")
 532            mesh.rotate_x(10).rotate_y(15).alpha(0.5)
 533            mesh.pointdata["scalars"] = mesh.vertices[:, 2]
 534
 535            printc("is_closed  :", mesh.is_closed())
 536            printc("is_manifold:", mesh.is_manifold())
 537            printc("euler_char :", mesh.euler_characteristic())
 538            printc("genus      :", mesh.genus())
 539
 540            reeb = mesh.to_reeb_graph()
 541            ids = reeb[0].pointdata["Vertex Ids"]
 542            pts = Points(mesh.vertices[ids], r=10)
 543
 544            show([[mesh, pts], reeb], N=2, sharecam=False)
 545            ```
 546        """
 547        rg = vtki.new("PolyDataToReebGraphFilter")
 548        rg.SetInputData(self.dataset)
 549        rg.SetFieldId(field_id)
 550        rg.Update()
 551        gr = vedo.pyplot.DirectedGraph()
 552        gr.mdg = rg.GetOutput()
 553        gr.build()
 554        return gr
 555
 556
 557    def shrink(self, fraction=0.85) -> Self:
 558        """
 559        Shrink the triangle polydata in the representation of the input mesh.
 560
 561        Examples:
 562            - [shrink.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/shrink.py)
 563
 564            ![](https://vedo.embl.es/images/basic/shrink.png)
 565        """
 566        # Overriding base class method core.shrink()
 567        shrink = vtki.new("ShrinkPolyData")
 568        shrink.SetInputData(self.dataset)
 569        shrink.SetShrinkFactor(fraction)
 570        shrink.Update()
 571        self._update(shrink.GetOutput())
 572        self.pipeline = OperationNode("shrink", parents=[self])
 573        return self
 574
 575    def cap(self, return_cap=False) -> Self:
 576        """
 577        Generate a "cap" on a clipped mesh, or caps sharp edges.
 578
 579        Examples:
 580            - [cut_and_cap.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/cut_and_cap.py)
 581
 582            ![](https://vedo.embl.es/images/advanced/cutAndCap.png)
 583
 584        See also: `join()`, `join_segments()`, `slice()`.
 585        """
 586        fe = vtki.new("FeatureEdges")
 587        fe.SetInputData(self.dataset)
 588        fe.BoundaryEdgesOn()
 589        fe.FeatureEdgesOff()
 590        fe.NonManifoldEdgesOff()
 591        fe.ManifoldEdgesOff()
 592        fe.Update()
 593
 594        stripper = vtki.new("Stripper")
 595        stripper.SetInputData(fe.GetOutput())
 596        stripper.JoinContiguousSegmentsOn()
 597        stripper.Update()
 598
 599        boundary_poly = vtki.vtkPolyData()
 600        boundary_poly.SetPoints(stripper.GetOutput().GetPoints())
 601        boundary_poly.SetPolys(stripper.GetOutput().GetLines())
 602
 603        rev = vtki.new("ReverseSense")
 604        rev.ReverseCellsOn()
 605        rev.SetInputData(boundary_poly)
 606        rev.Update()
 607
 608        tf = vtki.new("TriangleFilter")
 609        tf.SetInputData(rev.GetOutput())
 610        tf.Update()
 611
 612        if return_cap:
 613            m = Mesh(tf.GetOutput())
 614            m.pipeline = OperationNode(
 615                "cap", parents=[self], comment=f"#pts {m.dataset.GetNumberOfPoints()}"
 616            )
 617            m.name = "MeshCap"
 618            return m
 619
 620        polyapp = vtki.new("AppendPolyData")
 621        polyapp.AddInputData(self.dataset)
 622        polyapp.AddInputData(tf.GetOutput())
 623        polyapp.Update()
 624
 625        self._update(polyapp.GetOutput())
 626        self.clean()
 627
 628        self.pipeline = OperationNode(
 629            "capped", parents=[self], comment=f"#pts {self.dataset.GetNumberOfPoints()}"
 630        )
 631        return self
 632
 633    def join(self, polys=True, reset=False) -> Self:
 634        """
 635        Generate triangle strips and/or polylines from
 636        input polygons, triangle strips, and lines.
 637
 638        Input polygons are assembled into triangle strips only if they are triangles;
 639        other types of polygons are passed through to the output and not stripped.
 640        Use mesh.triangulate() to triangulate non-triangular polygons prior to running
 641        this filter if you need to strip all the data.
 642
 643        Also note that if triangle strips or polylines are present in the input
 644        they are passed through and not joined nor extended.
 645        If you wish to strip these use mesh.triangulate() to fragment the input
 646        into triangles and lines prior to applying join().
 647
 648        Arguments:
 649            polys : (bool)
 650                polygonal segments will be joined if they are contiguous
 651            reset : (bool)
 652                reset points ordering
 653
 654        Warning:
 655            If triangle strips or polylines exist in the input data
 656            they will be passed through to the output data.
 657            This filter will only construct triangle strips if triangle polygons
 658            are available; and will only construct polylines if lines are available.
 659
 660        Example:
 661            ```python
 662            from vedo import *
 663            c1 = Cylinder(pos=(0,0,0), r=2, height=3, axis=(1,.0,0), alpha=.1).triangulate()
 664            c2 = Cylinder(pos=(0,0,2), r=1, height=2, axis=(0,.3,1), alpha=.1).triangulate()
 665            intersect = c1.intersect_with(c2).join(reset=True)
 666            spline = Spline(intersect).c('blue').lw(5)
 667            show(c1, c2, spline, intersect.labels('id'), axes=1).close()
 668            ```
 669            ![](https://vedo.embl.es/images/feats/line_join.png)
 670        """
 671        sf = vtki.new("Stripper")
 672        sf.SetPassThroughCellIds(True)
 673        sf.SetPassThroughPointIds(True)
 674        sf.SetJoinContiguousSegments(polys)
 675        sf.SetInputData(self.dataset)
 676        sf.Update()
 677        if reset:
 678            poly = sf.GetOutput()
 679            cpd = vtki.new("CleanPolyData")
 680            cpd.PointMergingOn()
 681            cpd.ConvertLinesToPointsOn()
 682            cpd.ConvertPolysToLinesOn()
 683            cpd.ConvertStripsToPolysOn()
 684            cpd.SetInputData(poly)
 685            cpd.Update()
 686            poly = cpd.GetOutput()
 687            vpts = poly.GetCell(0).GetPoints().GetData()
 688            poly.GetPoints().SetData(vpts)
 689        else:
 690            poly = sf.GetOutput()
 691
 692        self._update(poly)
 693
 694        self.pipeline = OperationNode(
 695            "join", parents=[self], comment=f"#pts {self.dataset.GetNumberOfPoints()}"
 696        )
 697        return self
 698
 699    def join_segments(self, closed=True, tol=1e-03) -> list:
 700        """
 701        Join line segments into contiguous lines.
 702        Useful to call with `triangulate()` method.
 703
 704        Returns:
 705            list of `shapes.Lines`
 706
 707        Example:
 708            ```python
 709            from vedo import *
 710            msh = Torus().alpha(0.1).wireframe()
 711            intersection = msh.intersect_with_plane(normal=[1,1,1]).c('purple5')
 712            slices = [s.triangulate() for s in intersection.join_segments()]
 713            show(msh, intersection, merge(slices), axes=1, viewup='z')
 714            ```
 715            ![](https://vedo.embl.es/images/feats/join_segments.jpg)
 716        """
 717        vlines = []
 718        for ipiece, outline in enumerate(self.split(must_share_edge=False)): # type: ignore
 719
 720            outline.clean()
 721            pts = outline.vertices
 722            if len(pts) < 3:
 723                continue
 724            avesize = outline.average_size()
 725            lines = outline.lines
 726            # print("---lines", lines, "in piece", ipiece)
 727            tol = avesize / pts.shape[0] * tol
 728
 729            k = 0
 730            joinedpts = [pts[k]]
 731            for _ in range(len(pts)):
 732                pk = pts[k]
 733                for j, line in enumerate(lines):
 734
 735                    id0, id1 = line[0], line[-1]
 736                    p0, p1 = pts[id0], pts[id1]
 737
 738                    if np.linalg.norm(p0 - pk) < tol:
 739                        n = len(line)
 740                        for m in range(1, n):
 741                            joinedpts.append(pts[line[m]])
 742                        # joinedpts.append(p1)
 743                        k = id1
 744                        lines.pop(j)
 745                        break
 746
 747                    elif np.linalg.norm(p1 - pk) < tol:
 748                        n = len(line)
 749                        for m in reversed(range(0, n - 1)):
 750                            joinedpts.append(pts[line[m]])
 751                        # joinedpts.append(p0)
 752                        k = id0
 753                        lines.pop(j)
 754                        break
 755
 756            if len(joinedpts) > 1:
 757                newline = vedo.shapes.Line(joinedpts, closed=closed)
 758                newline.clean()
 759                newline.actor.SetProperty(self.properties)
 760                newline.properties = self.properties
 761                newline.pipeline = OperationNode(
 762                    "join_segments",
 763                    parents=[self],
 764                    comment=f"#pts {newline.dataset.GetNumberOfPoints()}",
 765                )
 766                vlines.append(newline)
 767
 768        return vlines
 769
 770    def join_with_strips(self, b1, closed=True) -> Self:
 771        """
 772        Join booundary lines by creating a triangle strip between them.
 773
 774        Example:
 775        ```python
 776        from vedo import *
 777        m1 = Cylinder(cap=False).boundaries()
 778        m2 = Cylinder(cap=False).boundaries().pos(0.2,0,1)
 779        strips = m1.join_with_strips(m2)
 780        show(m1, m2, strips, axes=1).close()
 781        ```
 782        """
 783        b0 = self.clone().join()
 784        b1 = b1.clone().join()
 785
 786        vertices0 = b0.vertices.tolist()
 787        vertices1 = b1.vertices.tolist()
 788
 789        lines0 = b0.lines
 790        lines1 = b1.lines
 791        m =  len(lines0)
 792        assert m == len(lines1), (
 793            "lines must have the same number of points\n"
 794            f"line has {m} points in b0 and {len(lines1)} in b1"
 795        )
 796
 797        strips = []
 798        points: List[Any] = []
 799
 800        for j in range(m):
 801
 802            ids0j = list(lines0[j])
 803            ids1j = list(lines1[j])
 804
 805            n = len(ids0j)
 806            assert n == len(ids1j), (
 807                "lines must have the same number of points\n"
 808                f"line {j} has {n} points in b0 and {len(ids1j)} in b1"
 809            )
 810
 811            if closed:
 812                ids0j.append(ids0j[0])
 813                ids1j.append(ids1j[0])
 814                vertices0.append(vertices0[ids0j[0]])
 815                vertices1.append(vertices1[ids1j[0]])
 816                n = n + 1
 817
 818            strip = []  # create a triangle strip
 819            npt = len(points)
 820            for ipt in range(n):
 821                points.append(vertices0[ids0j[ipt]])
 822                points.append(vertices1[ids1j[ipt]])
 823
 824            strip = list(range(npt, npt + 2*n))
 825            strips.append(strip)
 826
 827        return Mesh([points, [], [], strips], c="k6")
 828
 829    def split_polylines(self) -> Self:
 830        """Split polylines into separate segments."""
 831        tf = vtki.new("TriangleFilter")
 832        tf.SetPassLines(True)
 833        tf.SetPassVerts(False)
 834        tf.SetInputData(self.dataset)
 835        tf.Update()
 836        self._update(tf.GetOutput(), reset_locators=False)
 837        self.lw(0).lighting("default").pickable()
 838        self.pipeline = OperationNode(
 839            "split_polylines", parents=[self], 
 840            comment=f"#lines {self.dataset.GetNumberOfLines()}"
 841        )
 842        return self
 843
 844    def slice(self, origin=(0, 0, 0), normal=(1, 0, 0)) -> Self:
 845        """
 846        Slice a mesh with a plane and fill the contour.
 847
 848        Example:
 849            ```python
 850            from vedo import *
 851            msh = Mesh(dataurl+"bunny.obj").alpha(0.1).wireframe()
 852            mslice = msh.slice(normal=[0,1,0.3], origin=[0,0.16,0])
 853            mslice.c('purple5')
 854            show(msh, mslice, axes=1)
 855            ```
 856            ![](https://vedo.embl.es/images/feats/mesh_slice.jpg)
 857
 858        See also: `join()`, `join_segments()`, `cap()`, `cut_with_plane()`.
 859        """
 860        intersection = self.intersect_with_plane(origin=origin, normal=normal)
 861        slices = [s.triangulate() for s in intersection.join_segments()]
 862        mslices = vedo.pointcloud.merge(slices)
 863        if mslices:
 864            mslices.name = "MeshSlice"
 865            mslices.pipeline = OperationNode("slice", parents=[self], comment=f"normal = {normal}")
 866        return mslices
 867
 868    def triangulate(self, verts=True, lines=True) -> Self:
 869        """
 870        Converts mesh polygons into triangles.
 871
 872        If the input mesh is only made of 2D lines (no faces) the output will be a triangulation
 873        that fills the internal area. The contours may be concave, and may even contain holes,
 874        i.e. a contour may contain an internal contour winding in the opposite
 875        direction to indicate that it is a hole.
 876
 877        Arguments:
 878            verts : (bool)
 879                if True, break input vertex cells into individual vertex cells (one point per cell).
 880                If False, the input vertex cells will be ignored.
 881            lines : (bool)
 882                if True, break input polylines into line segments.
 883                If False, input lines will be ignored and the output will have no lines.
 884        """
 885        if self.dataset.GetNumberOfPolys() or self.dataset.GetNumberOfStrips():
 886            # print("Using vtkTriangleFilter")
 887            tf = vtki.new("TriangleFilter")
 888            tf.SetPassLines(lines)
 889            tf.SetPassVerts(verts)
 890
 891        elif self.dataset.GetNumberOfLines():
 892            # print("Using vtkContourTriangulator")
 893            tf = vtki.new("ContourTriangulator")
 894            tf.TriangulationErrorDisplayOn()
 895
 896        else:
 897            vedo.logger.debug("input in triangulate() seems to be void! Skip.")
 898            return self
 899
 900        tf.SetInputData(self.dataset)
 901        tf.Update()
 902        self._update(tf.GetOutput(), reset_locators=False)
 903        self.lw(0).lighting("default").pickable()
 904
 905        self.pipeline = OperationNode(
 906            "triangulate", parents=[self], comment=f"#cells {self.dataset.GetNumberOfCells()}"
 907        )
 908        return self
 909
 910    def compute_cell_vertex_count(self) -> Self:
 911        """
 912        Add to this mesh a cell data array containing the nr of vertices that a polygonal face has.
 913        """
 914        csf = vtki.new("CellSizeFilter")
 915        csf.SetInputData(self.dataset)
 916        csf.SetComputeArea(False)
 917        csf.SetComputeVolume(False)
 918        csf.SetComputeLength(False)
 919        csf.SetComputeVertexCount(True)
 920        csf.SetVertexCountArrayName("VertexCount")
 921        csf.Update()
 922        self.dataset.GetCellData().AddArray(
 923            csf.GetOutput().GetCellData().GetArray("VertexCount")
 924        )
 925        return self
 926
 927    def compute_quality(self, metric=6) -> Self:
 928        """
 929        Calculate metrics of quality for the elements of a triangular mesh.
 930        This method adds to the mesh a cell array named "Quality".
 931        See class 
 932        [vtkMeshQuality](https://vtk.org/doc/nightly/html/classvtkMeshQuality.html).
 933
 934        Arguments:
 935            metric : (int)
 936                type of available estimators are:
 937                - EDGE RATIO, 0
 938                - ASPECT RATIO, 1
 939                - RADIUS RATIO, 2
 940                - ASPECT FROBENIUS, 3
 941                - MED ASPECT FROBENIUS, 4
 942                - MAX ASPECT FROBENIUS, 5
 943                - MIN_ANGLE, 6
 944                - COLLAPSE RATIO, 7
 945                - MAX ANGLE, 8
 946                - CONDITION, 9
 947                - SCALED JACOBIAN, 10
 948                - SHEAR, 11
 949                - RELATIVE SIZE SQUARED, 12
 950                - SHAPE, 13
 951                - SHAPE AND SIZE, 14
 952                - DISTORTION, 15
 953                - MAX EDGE RATIO, 16
 954                - SKEW, 17
 955                - TAPER, 18
 956                - VOLUME, 19
 957                - STRETCH, 20
 958                - DIAGONAL, 21
 959                - DIMENSION, 22
 960                - ODDY, 23
 961                - SHEAR AND SIZE, 24
 962                - JACOBIAN, 25
 963                - WARPAGE, 26
 964                - ASPECT GAMMA, 27
 965                - AREA, 28
 966                - ASPECT BETA, 29
 967
 968        Examples:
 969            - [meshquality.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/meshquality.py)
 970
 971            ![](https://vedo.embl.es/images/advanced/meshquality.png)
 972        """
 973        qf = vtki.new("MeshQuality")
 974        qf.SetInputData(self.dataset)
 975        qf.SetTriangleQualityMeasure(metric)
 976        qf.SaveCellQualityOn()
 977        qf.Update()
 978        self._update(qf.GetOutput(), reset_locators=False)
 979        self.mapper.SetScalarModeToUseCellData()
 980        self.pipeline = OperationNode("compute_quality", parents=[self])
 981        return self
 982
 983    def count_vertices(self) -> np.ndarray:
 984        """Count the number of vertices each cell has and return it as a numpy array"""
 985        vc = vtki.new("CountVertices")
 986        vc.SetInputData(self.dataset)
 987        vc.SetOutputArrayName("VertexCount")
 988        vc.Update()
 989        varr = vc.GetOutput().GetCellData().GetArray("VertexCount")
 990        return vtk2numpy(varr)
 991
 992    def check_validity(self, tol=0) -> np.ndarray:
 993        """
 994        Return a numpy array of possible problematic faces following this convention:
 995        - Valid               =  0
 996        - WrongNumberOfPoints =  1
 997        - IntersectingEdges   =  2
 998        - IntersectingFaces   =  4
 999        - NoncontiguousEdges  =  8
1000        - Nonconvex           = 10
1001        - OrientedIncorrectly = 20
1002
1003        Arguments:
1004            tol : (float)
1005                value is used as an epsilon for floating point
1006                equality checks throughout the cell checking process.
1007        """
1008        vald = vtki.new("CellValidator")
1009        if tol:
1010            vald.SetTolerance(tol)
1011        vald.SetInputData(self.dataset)
1012        vald.Update()
1013        varr = vald.GetOutput().GetCellData().GetArray("ValidityState")
1014        return vtk2numpy(varr)
1015
1016    def compute_curvature(self, method=0) -> Self:
1017        """
1018        Add scalars to `Mesh` that contains the curvature calculated in three different ways.
1019
1020        Variable `method` can be:
1021        - 0 = gaussian
1022        - 1 = mean curvature
1023        - 2 = max curvature
1024        - 3 = min curvature
1025
1026        Example:
1027            ```python
1028            from vedo import Torus
1029            Torus().compute_curvature().add_scalarbar().show().close()
1030            ```
1031            ![](https://vedo.embl.es/images/advanced/torus_curv.png)
1032        """
1033        curve = vtki.new("Curvatures")
1034        curve.SetInputData(self.dataset)
1035        curve.SetCurvatureType(method)
1036        curve.Update()
1037        self._update(curve.GetOutput(), reset_locators=False)
1038        self.mapper.ScalarVisibilityOn()
1039        return self
1040
1041    def compute_elevation(self, low=(0, 0, 0), high=(0, 0, 1), vrange=(0, 1)) -> Self:
1042        """
1043        Add to `Mesh` a scalar array that contains distance along a specified direction.
1044
1045        Arguments:
1046            low : (list)
1047                one end of the line (small scalar values)
1048            high : (list)
1049                other end of the line (large scalar values)
1050            vrange : (list)
1051                set the range of the scalar
1052
1053        Example:
1054            ```python
1055            from vedo import Sphere
1056            s = Sphere().compute_elevation(low=(0,0,0), high=(1,1,1))
1057            s.add_scalarbar().show(axes=1).close()
1058            ```
1059            ![](https://vedo.embl.es/images/basic/compute_elevation.png)
1060        """
1061        ef = vtki.new("ElevationFilter")
1062        ef.SetInputData(self.dataset)
1063        ef.SetLowPoint(low)
1064        ef.SetHighPoint(high)
1065        ef.SetScalarRange(vrange)
1066        ef.Update()
1067        self._update(ef.GetOutput(), reset_locators=False)
1068        self.mapper.ScalarVisibilityOn()
1069        return self
1070
1071
1072    def laplacian_diffusion(self, array_name, dt, num_steps) -> Self:
1073        """
1074        Apply a diffusion process to a scalar array defined on the points of a mesh.
1075
1076        Arguments:
1077            array_name : (str)
1078                name of the array to diffuse.
1079            dt : (float)
1080                time step.
1081            num_steps : (int)
1082                number of iterations.
1083        """
1084        try:
1085            import scipy.sparse
1086            import scipy.sparse.linalg
1087        except ImportError:
1088            vedo.logger.error("scipy not found. Cannot run laplacian_diffusion()")
1089            return self
1090        
1091        def build_laplacian():
1092            rows = []
1093            cols = []
1094            data = []
1095            n_points = points.shape[0]
1096            avg_area = np.mean(areas) * 10000
1097            # print("avg_area", avg_area)
1098
1099            for triangle in cells:
1100                for i in range(3):
1101                    for j in range(i + 1, 3):
1102                        u = triangle[i]
1103                        v = triangle[j]
1104                        rows.append(u)
1105                        cols.append(v)
1106                        rows.append(v)
1107                        cols.append(u)
1108                        data.append(-1/avg_area)
1109                        data.append(-1/avg_area)
1110
1111            L = scipy.sparse.coo_matrix(
1112                (data, (rows, cols)), shape=(n_points, n_points)
1113            ).tocsc()
1114
1115            degree = -np.array(L.sum(axis=1)).flatten() # adjust the diagonal
1116            # print("degree", degree)
1117            L.setdiag(degree)
1118            return L
1119
1120        def _diffuse(u0, L, dt, num_steps):
1121            # mean_area = np.mean(areas) * 10000
1122            # print("mean_area", mean_area)
1123            mean_area = 1
1124            I = scipy.sparse.eye(L.shape[0], format="csc")
1125            A = I - (dt/mean_area) * L 
1126            u = u0
1127            for _ in range(int(num_steps)):
1128                u = A.dot(u)
1129            return u
1130
1131        self.compute_cell_size()
1132        areas = self.celldata["Area"]
1133        points = self.vertices
1134        cells = self.cells
1135        u0 = self.pointdata[array_name]
1136
1137        # Simulate diffusion
1138        L = build_laplacian()
1139        u = _diffuse(u0, L, dt, num_steps)
1140        self.pointdata[array_name] = u
1141        return self
1142
1143
1144    def subdivide(self, n=1, method=0, mel=None) -> Self:
1145        """
1146        Increase the number of vertices of a surface mesh.
1147
1148        Arguments:
1149            n : (int)
1150                number of subdivisions.
1151            method : (int)
1152                Loop(0), Linear(1), Adaptive(2), Butterfly(3), Centroid(4)
1153            mel : (float)
1154                Maximum Edge Length (applicable to Adaptive method only).
1155        """
1156        triangles = vtki.new("TriangleFilter")
1157        triangles.SetInputData(self.dataset)
1158        triangles.Update()
1159        tri_mesh = triangles.GetOutput()
1160        if method == 0:
1161            sdf = vtki.new("LoopSubdivisionFilter")
1162        elif method == 1:
1163            sdf = vtki.new("LinearSubdivisionFilter")
1164        elif method == 2:
1165            sdf = vtki.new("AdaptiveSubdivisionFilter")
1166            if mel is None:
1167                mel = self.diagonal_size() / np.sqrt(self.dataset.GetNumberOfPoints()) / n
1168            sdf.SetMaximumEdgeLength(mel)
1169        elif method == 3:
1170            sdf = vtki.new("ButterflySubdivisionFilter")
1171        elif method == 4:
1172            sdf = vtki.new("DensifyPolyData")
1173        else:
1174            vedo.logger.error(f"in subdivide() unknown method {method}")
1175            raise RuntimeError()
1176
1177        if method != 2:
1178            sdf.SetNumberOfSubdivisions(n)
1179
1180        sdf.SetInputData(tri_mesh)
1181        sdf.Update()
1182
1183        self._update(sdf.GetOutput())
1184
1185        self.pipeline = OperationNode(
1186            "subdivide",
1187            parents=[self],
1188            comment=f"#pts {self.dataset.GetNumberOfPoints()}",
1189        )
1190        return self
1191
1192
1193    def decimate(self, fraction=0.5, n=None, preserve_volume=True, regularization=0.0) -> Self:
1194        """
1195        Downsample the number of vertices in a mesh to `fraction`.
1196
1197        This filter preserves the `pointdata` of the input dataset. In previous versions
1198        of vedo, this decimation algorithm was referred to as quadric decimation.
1199
1200        Arguments:
1201            fraction : (float)
1202                the desired target of reduction.
1203            n : (int)
1204                the desired number of final points
1205                (`fraction` is recalculated based on it).
1206            preserve_volume : (bool)
1207                Decide whether to activate volume preservation which greatly
1208                reduces errors in triangle normal direction.
1209            regularization : (float)
1210                regularize the point finding algorithm so as to have better quality
1211                mesh elements at the cost of a slightly lower precision on the
1212                geometry potentially (mostly at sharp edges).
1213                Can be useful for decimating meshes that have been triangulated on noisy data.
1214
1215        Note:
1216            Setting `fraction=0.1` leaves 10% of the original number of vertices.
1217            Internally the VTK class
1218            [vtkQuadricDecimation](https://vtk.org/doc/nightly/html/classvtkQuadricDecimation.html)
1219            is used for this operation.
1220        
1221        See also: `decimate_binned()` and `decimate_pro()`.
1222        """
1223        poly = self.dataset
1224        if n:  # N = desired number of points
1225            npt = poly.GetNumberOfPoints()
1226            fraction = n / npt
1227            if fraction >= 1:
1228                return self
1229
1230        decimate = vtki.new("QuadricDecimation")
1231        decimate.SetVolumePreservation(preserve_volume)
1232        # decimate.AttributeErrorMetricOn()
1233        if regularization:
1234            decimate.SetRegularize(True)
1235            decimate.SetRegularization(regularization)
1236
1237        try:
1238            decimate.MapPointDataOn()
1239        except AttributeError:
1240            pass
1241
1242        decimate.SetTargetReduction(1 - fraction)
1243        decimate.SetInputData(poly)
1244        decimate.Update()
1245
1246        self._update(decimate.GetOutput())
1247        self.metadata["decimate_actual_fraction"] = 1 - decimate.GetActualReduction()
1248
1249        self.pipeline = OperationNode(
1250            "decimate",
1251            parents=[self],
1252            comment=f"#pts {self.dataset.GetNumberOfPoints()}",
1253        )
1254        return self
1255    
1256    def decimate_pro(
1257            self,
1258            fraction=0.5,
1259            n=None,
1260            preserve_topology=True,
1261            preserve_boundaries=True,
1262            splitting=False,
1263            splitting_angle=75,
1264            feature_angle=0,
1265            inflection_point_ratio=10,
1266            vertex_degree=0,
1267        ) -> Self:
1268        """
1269        Downsample the number of vertices in a mesh to `fraction`.
1270
1271        This filter preserves the `pointdata` of the input dataset.
1272
1273        Arguments:
1274            fraction : (float)
1275                The desired target of reduction.
1276                Setting `fraction=0.1` leaves 10% of the original number of vertices.
1277            n : (int)
1278                the desired number of final points (`fraction` is recalculated based on it).
1279            preserve_topology : (bool)
1280                If on, mesh splitting and hole elimination will not occur.
1281                This may limit the maximum reduction that may be achieved.
1282            preserve_boundaries : (bool)
1283                Turn on/off the deletion of vertices on the boundary of a mesh.
1284                Control whether mesh boundaries are preserved during decimation.
1285            feature_angle : (float)
1286                Specify the angle that defines a feature.
1287                This angle is used to define what an edge is
1288                (i.e., if the surface normal between two adjacent triangles
1289                is >= FeatureAngle, an edge exists).
1290            splitting : (bool)
1291                Turn on/off the splitting of the mesh at corners,
1292                along edges, at non-manifold points, or anywhere else a split is required.
1293                Turning splitting off will better preserve the original topology of the mesh,
1294                but you may not obtain the requested reduction.
1295            splitting_angle : (float)
1296                Specify the angle that defines a sharp edge.
1297                This angle is used to control the splitting of the mesh.
1298                A split line exists when the surface normals between two edge connected triangles
1299                are >= `splitting_angle`.
1300            inflection_point_ratio : (float)
1301                An inflection point occurs when the ratio of reduction error between two iterations
1302                is greater than or equal to the `inflection_point_ratio` value.
1303            vertex_degree : (int)
1304                If the number of triangles connected to a vertex exceeds it then the vertex will be split.
1305
1306        Note:
1307            Setting `fraction=0.1` leaves 10% of the original number of vertices
1308        
1309        See also:
1310            `decimate()` and `decimate_binned()`.
1311        """
1312        poly = self.dataset
1313        if n:  # N = desired number of points
1314            npt = poly.GetNumberOfPoints()
1315            fraction = n / npt
1316            if fraction >= 1:
1317                return self
1318
1319        decimate = vtki.new("DecimatePro")
1320        decimate.SetPreserveTopology(preserve_topology)
1321        decimate.SetBoundaryVertexDeletion(preserve_boundaries)
1322        if feature_angle:
1323            decimate.SetFeatureAngle(feature_angle)
1324        decimate.SetSplitting(splitting)
1325        decimate.SetSplitAngle(splitting_angle)
1326        decimate.SetInflectionPointRatio(inflection_point_ratio)
1327        if vertex_degree:
1328            decimate.SetDegree(vertex_degree)
1329
1330        decimate.SetTargetReduction(1 - fraction)
1331        decimate.SetInputData(poly)
1332        decimate.Update()
1333        self._update(decimate.GetOutput())
1334
1335        self.pipeline = OperationNode(
1336            "decimate_pro",
1337            parents=[self],
1338            comment=f"#pts {self.dataset.GetNumberOfPoints()}",
1339        )
1340        return self
1341    
1342    def decimate_binned(self, divisions=(), use_clustering=False) -> Self:
1343        """
1344        Downsample the number of vertices in a mesh.
1345        
1346        This filter preserves the `celldata` of the input dataset,
1347        if `use_clustering=True` also the `pointdata` will be preserved in the result.
1348
1349        Arguments:
1350            divisions : (list)
1351                number of divisions along x, y and z axes.
1352            auto_adjust : (bool)
1353                if True, the number of divisions is automatically adjusted to
1354                create more uniform cells.
1355            use_clustering : (bool)
1356                use [vtkQuadricClustering](https://vtk.org/doc/nightly/html/classvtkQuadricClustering.html)
1357                instead of 
1358                [vtkBinnedDecimation](https://vtk.org/doc/nightly/html/classvtkBinnedDecimation.html).
1359        
1360        See also: `decimate()` and `decimate_pro()`.
1361        """
1362        if use_clustering:
1363            decimate = vtki.new("QuadricClustering")
1364            decimate.CopyCellDataOn()
1365        else:
1366            decimate = vtki.new("BinnedDecimation")
1367            decimate.ProducePointDataOn()
1368            decimate.ProduceCellDataOn()
1369
1370        decimate.SetInputData(self.dataset)
1371
1372        if len(divisions) == 0:
1373            decimate.SetAutoAdjustNumberOfDivisions(1)
1374        else:
1375            decimate.SetAutoAdjustNumberOfDivisions(0)
1376            decimate.SetNumberOfDivisions(divisions)
1377        decimate.Update()
1378
1379        self._update(decimate.GetOutput())
1380        self.metadata["decimate_binned_divisions"] = decimate.GetNumberOfDivisions()
1381        self.pipeline = OperationNode(
1382            "decimate_binned",
1383            parents=[self],
1384            comment=f"#pts {self.dataset.GetNumberOfPoints()}",
1385        )
1386        return self
1387
1388    def generate_random_points(self, n: int, min_radius=0.0) -> "Points":
1389        """
1390        Generate `n` uniformly distributed random points
1391        inside the polygonal mesh.
1392
1393        A new point data array is added to the output points
1394        called "OriginalCellID" which contains the index of
1395        the cell ID in which the point was generated.
1396
1397        Arguments:
1398            n : (int)
1399                number of points to generate.
1400            min_radius: (float)
1401                impose a minimum distance between points.
1402                If `min_radius` is set to 0, the points are
1403                generated uniformly at random inside the mesh.
1404                If `min_radius` is set to a positive value,
1405                the points are generated uniformly at random
1406                inside the mesh, but points closer than `min_radius`
1407                to any other point are discarded.
1408
1409        Returns a `vedo.Points` object.
1410
1411        Note:
1412            Consider using `points.probe(msh)` or
1413            `points.interpolate_data_from(msh)`
1414            to interpolate existing mesh data onto the new points.
1415
1416        Example:
1417        ```python
1418        from vedo import *
1419        msh = Mesh(dataurl + "panther.stl").lw(2)
1420        pts = msh.generate_random_points(20000, min_radius=0.5)
1421        print("Original cell ids:", pts.pointdata["OriginalCellID"])
1422        show(pts, msh, axes=1).close()
1423        ```
1424        """
1425        cmesh = self.clone().clean().triangulate().compute_cell_size()
1426        triangles = cmesh.cells
1427        vertices = cmesh.vertices
1428        cumul = np.cumsum(cmesh.celldata["Area"])
1429
1430        out_pts = []
1431        orig_cell = []
1432        for _ in range(n):
1433            # choose a triangle based on area
1434            random_area = np.random.random() * cumul[-1]
1435            it = np.searchsorted(cumul, random_area)
1436            A, B, C = vertices[triangles[it]]
1437            # calculate the random point in the triangle
1438            r1, r2 = np.random.random(2)
1439            if r1 + r2 > 1:
1440                r1 = 1 - r1
1441                r2 = 1 - r2
1442            out_pts.append((1 - r1 - r2) * A + r1 * B + r2 * C)
1443            orig_cell.append(it)
1444        nporig_cell = np.array(orig_cell, dtype=np.uint32)
1445
1446        vpts = Points(out_pts)
1447        vpts.pointdata["OriginalCellID"] = nporig_cell
1448
1449        if min_radius > 0:
1450            vpts.subsample(min_radius, absolute=True)
1451
1452        vpts.point_size(5).color("k1")
1453        vpts.name = "RandomPoints"
1454        vpts.pipeline = OperationNode(
1455            "generate_random_points", c="#edabab", parents=[self])
1456        return vpts
1457
1458    def delete_cells(self, ids: List[int]) -> Self:
1459        """
1460        Remove cells from the mesh object by their ID.
1461        Points (vertices) are not removed (you may use `clean()` to remove those).
1462        """
1463        self.dataset.BuildLinks()
1464        for cid in ids:
1465            self.dataset.DeleteCell(cid)
1466        self.dataset.RemoveDeletedCells()
1467        self.dataset.Modified()
1468        self.mapper.Modified()
1469        self.pipeline = OperationNode(
1470            "delete_cells",
1471            parents=[self],
1472            comment=f"#cells {self.dataset.GetNumberOfCells()}",
1473        )
1474        return self
1475
1476    def delete_cells_by_point_index(self, indices: List[int]) -> Self:
1477        """
1478        Delete a list of vertices identified by any of their vertex index.
1479
1480        See also `delete_cells()`.
1481
1482        Examples:
1483            - [delete_mesh_pts.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/delete_mesh_pts.py)
1484
1485                ![](https://vedo.embl.es/images/basic/deleteMeshPoints.png)
1486        """
1487        cell_ids = vtki.vtkIdList()
1488        self.dataset.BuildLinks()
1489        n = 0
1490        for i in np.unique(indices):
1491            self.dataset.GetPointCells(i, cell_ids)
1492            for j in range(cell_ids.GetNumberOfIds()):
1493                self.dataset.DeleteCell(cell_ids.GetId(j))  # flag cell
1494                n += 1
1495
1496        self.dataset.RemoveDeletedCells()
1497        self.dataset.Modified()
1498        self.pipeline = OperationNode("delete_cells_by_point_index", parents=[self])
1499        return self
1500
1501    def collapse_edges(self, distance: float, iterations=1) -> Self:
1502        """
1503        Collapse mesh edges so that are all above `distance`.
1504        
1505        Example:
1506            ```python
1507            from vedo import *
1508            np.random.seed(2)
1509            grid1 = Grid().add_gaussian_noise(0.8).triangulate().lw(1)
1510            grid1.celldata['scalar'] = grid1.cell_centers[:,1]
1511            grid2 = grid1.clone().collapse_edges(0.1)
1512            show(grid1, grid2, N=2, axes=1)
1513            ```
1514        """
1515        for _ in range(iterations):
1516            medges = self.edges
1517            pts = self.vertices
1518            newpts = np.array(pts)
1519            moved = []
1520            for e in medges:
1521                if len(e) == 2:
1522                    id0, id1 = e
1523                    p0, p1 = pts[id0], pts[id1]
1524                    if (np.linalg.norm(p1-p0) < distance 
1525                        and id0 not in moved
1526                        and id1 not in moved
1527                    ):
1528                        p = (p0 + p1) / 2
1529                        newpts[id0] = p
1530                        newpts[id1] = p
1531                        moved += [id0, id1]
1532            self.vertices = newpts
1533            cpd = vtki.new("CleanPolyData")
1534            cpd.ConvertLinesToPointsOff()
1535            cpd.ConvertPolysToLinesOff()
1536            cpd.ConvertStripsToPolysOff()
1537            cpd.SetInputData(self.dataset)
1538            cpd.Update()
1539            self._update(cpd.GetOutput())
1540
1541        self.pipeline = OperationNode(
1542            "collapse_edges",
1543            parents=[self],
1544            comment=f"#pts {self.dataset.GetNumberOfPoints()}",
1545        )
1546        return self
1547
1548    def adjacency_list(self) -> List[set]:
1549        """
1550        Computes the adjacency list for mesh edge-graph.
1551
1552        Returns: 
1553            a list with i-th entry being the set if indices of vertices connected by an edge to i-th vertex
1554        """
1555        inc = [set()] * self.nvertices
1556        for cell in self.cells:
1557            nc = len(cell)
1558            if nc > 1:
1559                for i in range(nc-1):
1560                    ci = cell[i]
1561                    inc[ci] = inc[ci].union({cell[i-1], cell[i+1]})
1562        return inc
1563
1564    def graph_ball(self, index, n: int) -> set:
1565        """
1566        Computes the ball of radius `n` in the mesh' edge-graph metric centred in vertex `index`.
1567
1568        Arguments:
1569            index : (int)
1570                index of the vertex
1571            n : (int)
1572                radius in the graph metric
1573
1574        Returns:
1575            the set of indices of the vertices which are at most `n` edges from vertex `index`.
1576        """
1577        if n == 0:
1578            return {index}
1579        else:
1580            al = self.adjacency_list()
1581            ball = {index}
1582            i = 0
1583            while i < n and len(ball) < self.nvertices:
1584                for v in ball:
1585                    ball = ball.union(al[v])
1586                i += 1
1587            return ball
1588
1589    def smooth(self, niter=15, pass_band=0.1, edge_angle=15, feature_angle=60, boundary=False) -> Self:
1590        """
1591        Adjust mesh point positions using the so-called "Windowed Sinc" method.
1592
1593        Arguments:
1594            niter : (int)
1595                number of iterations.
1596            pass_band : (float)
1597                set the pass_band value for the windowed sinc filter.
1598            edge_angle : (float)
1599                edge angle to control smoothing along edges (either interior or boundary).
1600            feature_angle : (float)
1601                specifies the feature angle for sharp edge identification.
1602            boundary : (bool)
1603                specify if boundary should also be smoothed or kept unmodified
1604
1605        Examples: