vedo.mesh

Submodule to work with polygonal meshes

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