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