vedo.mesh

Submodule to work with polygonal meshes

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