vedo.mesh

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