vedo.pointcloud

Submodule to work with point clouds.

   1#!/usr/bin/env python3
   2# -*- coding: utf-8 -*-
   3import time
   4from typing import Union, List
   5from typing_extensions import Self
   6from weakref import ref as weak_ref_to
   7import numpy as np
   8
   9import vedo.vtkclasses as vtki
  10
  11import vedo
  12from vedo import colors
  13from vedo import utils
  14from vedo.transformations import LinearTransform, NonLinearTransform
  15from vedo.core import PointAlgorithms
  16from vedo.visual import PointsVisual
  17
  18__docformat__ = "google"
  19
  20__doc__ = """
  21Submodule to work with point clouds.
  22
  23![](https://vedo.embl.es/images/basic/pca.png)
  24"""
  25
  26__all__ = [
  27    "Points",
  28    "Point",
  29    "merge",
  30    "fit_line",
  31    "fit_circle",
  32    "fit_plane",
  33    "fit_sphere",
  34    "pca_ellipse",
  35    "pca_ellipsoid",
  36]
  37
  38
  39####################################################
  40def merge(*meshs, flag=False) -> Union["vedo.Mesh", "vedo.Points", None]:
  41    """
  42    Build a new Mesh (or Points) formed by the fusion of the inputs.
  43
  44    Similar to Assembly, but in this case the input objects become a single entity.
  45
  46    To keep track of the original identities of the inputs you can set `flag=True`.
  47    In this case a `pointdata` array of ids is added to the output with name "OriginalMeshID".
  48
  49    Examples:
  50        - [warp1.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/warp1.py)
  51
  52            ![](https://vedo.embl.es/images/advanced/warp1.png)
  53
  54        - [value_iteration.py](https://github.com/marcomusy/vedo/tree/master/examples/simulations/value_iteration.py)
  55
  56    """
  57    objs = [a for a in utils.flatten(meshs) if a]
  58
  59    if not objs:
  60        return None
  61
  62    idarr = []
  63    polyapp = vtki.new("AppendPolyData")
  64    for i, ob in enumerate(objs):
  65        polyapp.AddInputData(ob.dataset)
  66        if flag:
  67            idarr += [i] * ob.dataset.GetNumberOfPoints()
  68    polyapp.Update()
  69    mpoly = polyapp.GetOutput()
  70
  71    if flag:
  72        varr = utils.numpy2vtk(idarr, dtype=np.uint16, name="OriginalMeshID")
  73        mpoly.GetPointData().AddArray(varr)
  74
  75    has_mesh = False
  76    for ob in objs:
  77        if isinstance(ob, vedo.Mesh):
  78            has_mesh = True
  79            break
  80
  81    if has_mesh:
  82        msh = vedo.Mesh(mpoly)
  83    else:
  84        msh = Points(mpoly) # type: ignore
  85
  86    msh.copy_properties_from(objs[0])
  87
  88    msh.pipeline = utils.OperationNode(
  89        "merge", parents=objs, comment=f"#pts {msh.dataset.GetNumberOfPoints()}"
  90    )
  91    return msh
  92
  93
  94def _rotate_points(points, n0=None, n1=(0, 0, 1)) -> Union[np.ndarray, tuple]:
  95    # Rotate a set of 3D points from direction n0 to direction n1.
  96    # Return the rotated points and the normal to the fitting plane (if n0 is None).
  97    # The pointing direction of the normal in this case is arbitrary.
  98    points = np.asarray(points)
  99
 100    if points.ndim == 1:
 101        points = points[np.newaxis, :]
 102
 103    if len(points[0]) == 2:
 104        return points, (0, 0, 1)
 105
 106    if n0 is None:  # fit plane
 107        datamean = points.mean(axis=0)
 108        vv = np.linalg.svd(points - datamean)[2]
 109        n0 = np.cross(vv[0], vv[1])
 110
 111    n0 = n0 / np.linalg.norm(n0)
 112    n1 = n1 / np.linalg.norm(n1)
 113    k = np.cross(n0, n1)
 114    l = np.linalg.norm(k)
 115    if not l:
 116        k = n0
 117    k /= np.linalg.norm(k)
 118
 119    ct = np.dot(n0, n1)
 120    theta = np.arccos(ct)
 121    st = np.sin(theta)
 122    v = k * (1 - ct)
 123
 124    rpoints = []
 125    for p in points:
 126        a = p * ct
 127        b = np.cross(k, p) * st
 128        c = v * np.dot(k, p)
 129        rpoints.append(a + b + c)
 130
 131    return np.array(rpoints), n0
 132
 133
 134def fit_line(points: Union[np.ndarray, "vedo.Points"]) -> "vedo.shapes.Line":
 135    """
 136    Fits a line through points.
 137
 138    Extra info is stored in `Line.slope`, `Line.center`, `Line.variances`.
 139
 140    Examples:
 141        - [fitline.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/fitline.py)
 142
 143            ![](https://vedo.embl.es/images/advanced/fitline.png)
 144    """
 145    if isinstance(points, Points):
 146        points = points.coordinates
 147    data = np.asarray(points)
 148    datamean = data.mean(axis=0)
 149    _, dd, vv = np.linalg.svd(data - datamean)
 150    vv = vv[0] / np.linalg.norm(vv[0])
 151    # vv contains the first principal component, i.e. the direction
 152    # vector of the best fit line in the least squares sense.
 153    xyz_min = data.min(axis=0)
 154    xyz_max = data.max(axis=0)
 155    a = np.linalg.norm(xyz_min - datamean)
 156    b = np.linalg.norm(xyz_max - datamean)
 157    p1 = datamean - a * vv
 158    p2 = datamean + b * vv
 159    line = vedo.shapes.Line(p1, p2, lw=1)
 160    line.slope = vv
 161    line.center = datamean
 162    line.variances = dd
 163    return line
 164
 165
 166def fit_circle(points: Union[np.ndarray, "vedo.Points"]) -> tuple:
 167    """
 168    Fits a circle through a set of 3D points, with a very fast non-iterative method.
 169
 170    Returns the tuple `(center, radius, normal_to_circle)`.
 171
 172    .. warning::
 173        trying to fit s-shaped points will inevitably lead to instabilities and
 174        circles of small radius.
 175
 176    References:
 177        *J.F. Crawford, Nucl. Instr. Meth. 211, 1983, 223-225.*
 178    """
 179    if isinstance(points, Points):
 180        points = points.coordinates
 181    data = np.asarray(points)
 182
 183    offs = data.mean(axis=0)
 184    data, n0 = _rotate_points(data - offs)
 185
 186    xi = data[:, 0]
 187    yi = data[:, 1]
 188
 189    x = sum(xi)
 190    xi2 = xi * xi
 191    xx = sum(xi2)
 192    xxx = sum(xi2 * xi)
 193
 194    y = sum(yi)
 195    yi2 = yi * yi
 196    yy = sum(yi2)
 197    yyy = sum(yi2 * yi)
 198
 199    xiyi = xi * yi
 200    xy = sum(xiyi)
 201    xyy = sum(xiyi * yi)
 202    xxy = sum(xi * xiyi)
 203
 204    N = len(xi)
 205    k = (xx + yy) / N
 206
 207    a1 = xx - x * x / N
 208    b1 = xy - x * y / N
 209    c1 = 0.5 * (xxx + xyy - x * k)
 210
 211    a2 = xy - x * y / N
 212    b2 = yy - y * y / N
 213    c2 = 0.5 * (xxy + yyy - y * k)
 214
 215    d = a2 * b1 - a1 * b2
 216    if not d:
 217        return offs, 0, n0
 218    x0 = (b1 * c2 - b2 * c1) / d
 219    y0 = (c1 - a1 * x0) / b1
 220
 221    R = np.sqrt(x0 * x0 + y0 * y0 - 1 / N * (2 * x0 * x + 2 * y0 * y - xx - yy))
 222
 223    c, _ = _rotate_points([x0, y0, 0], (0, 0, 1), n0)
 224
 225    return c[0] + offs, R, n0
 226
 227
 228def fit_plane(points: Union[np.ndarray, "vedo.Points"], signed=False) -> "vedo.shapes.Plane":
 229    """
 230    Fits a plane to a set of points.
 231
 232    Extra info is stored in `Plane.normal`, `Plane.center`, `Plane.variance`.
 233
 234    Arguments:
 235        signed : (bool)
 236            if True flip sign of the normal based on the ordering of the points
 237
 238    Examples:
 239        - [fitline.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/fitline.py)
 240
 241            ![](https://vedo.embl.es/images/advanced/fitline.png)
 242    """
 243    if isinstance(points, Points):
 244        points = points.coordinates
 245    data = np.asarray(points)
 246    datamean = data.mean(axis=0)
 247    pts = data - datamean
 248    res = np.linalg.svd(pts)
 249    dd, vv = res[1], res[2]
 250    n = np.cross(vv[0], vv[1])
 251    if signed:
 252        v = np.zeros_like(pts)
 253        for i in range(len(pts) - 1):
 254            vi = np.cross(pts[i], pts[i + 1])
 255            v[i] = vi / np.linalg.norm(vi)
 256        ns = np.mean(v, axis=0)  # normal to the points plane
 257        if np.dot(n, ns) < 0:
 258            n = -n
 259    xyz_min = data.min(axis=0)
 260    xyz_max = data.max(axis=0)
 261    s = np.linalg.norm(xyz_max - xyz_min)
 262    pla = vedo.shapes.Plane(datamean, n, s=[s, s])
 263    pla.variance = dd[2]
 264    pla.name = "FitPlane"
 265    return pla
 266
 267
 268def fit_sphere(coords: Union[np.ndarray, "vedo.Points"]) -> "vedo.shapes.Sphere":
 269    """
 270    Fits a sphere to a set of points.
 271
 272    Extra info is stored in `Sphere.radius`, `Sphere.center`, `Sphere.residue`.
 273
 274    Examples:
 275        - [fitspheres1.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/fitspheres1.py)
 276
 277            ![](https://vedo.embl.es/images/advanced/fitspheres1.jpg)
 278    """
 279    if isinstance(coords, Points):
 280        coords = coords.coordinates
 281    coords = np.array(coords)
 282    n = len(coords)
 283    A = np.zeros((n, 4))
 284    A[:, :-1] = coords * 2
 285    A[:, 3] = 1
 286    f = np.zeros((n, 1))
 287    x = coords[:, 0]
 288    y = coords[:, 1]
 289    z = coords[:, 2]
 290    f[:, 0] = x * x + y * y + z * z
 291    try:
 292        C, residue, rank, _ = np.linalg.lstsq(A, f, rcond=-1)  # solve AC=f
 293    except:
 294        C, residue, rank, _ = np.linalg.lstsq(A, f)  # solve AC=f
 295    if rank < 4:
 296        return None
 297    t = (C[0] * C[0]) + (C[1] * C[1]) + (C[2] * C[2]) + C[3]
 298    radius = np.sqrt(t)[0]
 299    center = np.array([C[0][0], C[1][0], C[2][0]])
 300    if len(residue) > 0:
 301        residue = np.sqrt(residue[0]) / n
 302    else:
 303        residue = 0
 304    sph = vedo.shapes.Sphere(center, radius, c=(1, 0, 0)).wireframe(1)
 305    sph.radius = radius
 306    sph.center = center
 307    sph.residue = residue
 308    sph.name = "FitSphere"
 309    return sph
 310
 311
 312def pca_ellipse(points: Union[np.ndarray, "vedo.Points"], pvalue=0.673, res=60) -> Union["vedo.shapes.Circle", None]:
 313    """
 314    Create the oriented 2D ellipse that contains the fraction `pvalue` of points.
 315    PCA (Principal Component Analysis) is used to compute the ellipse orientation.
 316
 317    Parameter `pvalue` sets the specified fraction of points inside the ellipse.
 318    Normalized directions are stored in `ellipse.axis1`, `ellipse.axis2`.
 319    Axes sizes are stored in `ellipse.va`, `ellipse.vb`
 320
 321    Arguments:
 322        pvalue : (float)
 323            ellipse will include this fraction of points
 324        res : (int)
 325            resolution of the ellipse
 326
 327    Examples:
 328        - [pca_ellipse.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/pca_ellipse.py)
 329        - [histo_pca.py](https://github.com/marcomusy/vedo/tree/master/examples/pyplot/histo_pca.py)
 330
 331            ![](https://vedo.embl.es/images/pyplot/histo_pca.png)
 332    """
 333    from scipy.stats import f
 334
 335    if isinstance(points, Points):
 336        coords = points.coordinates
 337    else:
 338        coords = points
 339    if len(coords) < 4:
 340        vedo.logger.warning("in pca_ellipse(), there are not enough points!")
 341        return None
 342
 343    P = np.array(coords, dtype=float)[:, (0, 1)]
 344    cov = np.cov(P, rowvar=0)      # type: ignore
 345    _, s, R = np.linalg.svd(cov)   # singular value decomposition
 346    p, n = s.size, P.shape[0]
 347    fppf = f.ppf(pvalue, p, n - p) # f % point function
 348    u = np.sqrt(s * fppf / 2) * 2  # semi-axes (largest first)
 349    ua, ub = u
 350    center = utils.make3d(np.mean(P, axis=0)) # centroid of the ellipse
 351
 352    t = LinearTransform(R.T * u).translate(center)
 353    elli = vedo.shapes.Circle(alpha=0.75, res=res)
 354    elli.apply_transform(t)
 355    elli.properties.LightingOff()
 356
 357    elli.pvalue = pvalue
 358    elli.center = np.array([center[0], center[1], 0])
 359    elli.nr_of_points = n
 360    elli.va = ua
 361    elli.vb = ub
 362    
 363    # we subtract center because it's in t
 364    elli.axis1 = t.move([1, 0, 0]) - center
 365    elli.axis2 = t.move([0, 1, 0]) - center
 366
 367    elli.axis1 /= np.linalg.norm(elli.axis1)
 368    elli.axis2 /= np.linalg.norm(elli.axis2)
 369    elli.name = "PCAEllipse"
 370    return elli
 371
 372
 373def pca_ellipsoid(points: Union[np.ndarray, "vedo.Points"], pvalue=0.673, res=24) -> Union["vedo.shapes.Ellipsoid", None]:
 374    """
 375    Create the oriented ellipsoid that contains the fraction `pvalue` of points.
 376    PCA (Principal Component Analysis) is used to compute the ellipsoid orientation.
 377
 378    Axes sizes can be accessed in `ellips.va`, `ellips.vb`, `ellips.vc`,
 379    normalized directions are stored in `ellips.axis1`, `ellips.axis2` and `ellips.axis3`.
 380    Center of mass is stored in `ellips.center`.
 381
 382    Asphericity can be accessed in `ellips.asphericity()` and ellips.asphericity_error().
 383    A value of 0 means a perfect sphere.
 384
 385    Arguments:
 386        pvalue : (float)
 387            ellipsoid will include this fraction of points
 388   
 389    Examples:
 390        [pca_ellipsoid.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/pca_ellipsoid.py)
 391
 392            ![](https://vedo.embl.es/images/basic/pca.png)
 393    
 394    See also:
 395        `pca_ellipse()` for a 2D ellipse.
 396    """
 397    from scipy.stats import f
 398
 399    if isinstance(points, Points):
 400        coords = points.coordinates
 401    else:
 402        coords = points
 403    if len(coords) < 4:
 404        vedo.logger.warning("in pca_ellipsoid(), not enough input points!")
 405        return None
 406
 407    P = np.array(coords, ndmin=2, dtype=float)
 408    cov = np.cov(P, rowvar=0)     # type: ignore
 409    _, s, R = np.linalg.svd(cov)  # singular value decomposition
 410    p, n = s.size, P.shape[0]
 411    fppf = f.ppf(pvalue, p, n-p)*(n-1)*p*(n+1)/n/(n-p)  # f % point function
 412    u = np.sqrt(s*fppf)
 413    ua, ub, uc = u                # semi-axes (largest first)
 414    center = np.mean(P, axis=0)   # centroid of the hyperellipsoid
 415
 416    t = LinearTransform(R.T * u).translate(center)
 417    elli = vedo.shapes.Ellipsoid((0,0,0), (1,0,0), (0,1,0), (0,0,1), res=res)
 418    elli.apply_transform(t)
 419    elli.alpha(0.25)
 420    elli.properties.LightingOff()
 421
 422    elli.pvalue = pvalue
 423    elli.nr_of_points = n
 424    elli.center = center
 425    elli.va = ua
 426    elli.vb = ub
 427    elli.vc = uc
 428    # we subtract center because it's in t
 429    elli.axis1 = np.array(t.move([1, 0, 0])) - center
 430    elli.axis2 = np.array(t.move([0, 1, 0])) - center
 431    elli.axis3 = np.array(t.move([0, 0, 1])) - center
 432    elli.axis1 /= np.linalg.norm(elli.axis1)
 433    elli.axis2 /= np.linalg.norm(elli.axis2)
 434    elli.axis3 /= np.linalg.norm(elli.axis3)
 435    elli.name = "PCAEllipsoid"
 436    return elli
 437
 438
 439###################################################
 440def Point(pos=(0, 0, 0), r=12, c="red", alpha=1.0) -> Self:
 441    """
 442    Create a simple point in space.
 443
 444    .. note:: if you are creating many points you should use class `Points` instead!
 445    """
 446    pt = Points([[0,0,0]], r, c, alpha).pos(pos)
 447    pt.name = "Point"
 448    return pt
 449
 450
 451###################################################
 452class Points(PointsVisual, PointAlgorithms):
 453    """Work with point clouds."""
 454
 455    def __init__(self, inputobj=None, r=4, c=(0.2, 0.2, 0.2), alpha=1):
 456        """
 457        Build an object made of only vertex points for a list of 2D/3D points.
 458        Both shapes (N, 3) or (3, N) are accepted as input, if N>3.
 459
 460        Arguments:
 461            inputobj : (list, tuple)
 462            r : (int)
 463                Point radius in units of pixels.
 464            c : (str, list)
 465                Color name or rgb tuple.
 466            alpha : (float)
 467                Transparency in range [0,1].
 468
 469        Example:
 470            ```python
 471            from vedo import *
 472
 473            def fibonacci_sphere(n):
 474                s = np.linspace(0, n, num=n, endpoint=False)
 475                theta = s * 2.399963229728653
 476                y = 1 - s * (2/(n-1))
 477                r = np.sqrt(1 - y * y)
 478                x = np.cos(theta) * r
 479                z = np.sin(theta) * r
 480                return np._c[x,y,z]
 481
 482            Points(fibonacci_sphere(1000)).show(axes=1).close()
 483            ```
 484            ![](https://vedo.embl.es/images/feats/fibonacci.png)
 485        """
 486        # print("INIT POINTS")
 487        super().__init__()
 488
 489        self.name = ""
 490        self.filename = ""
 491        self.file_size = ""
 492
 493        self.info = {}
 494        self.time = time.time()
 495        
 496        self.transform = LinearTransform()
 497        self.point_locator = None
 498        self.cell_locator = None
 499        self.line_locator = None
 500
 501        self.actor = vtki.vtkActor()
 502        self.properties = self.actor.GetProperty()
 503        self.properties_backface = self.actor.GetBackfaceProperty()
 504        self.mapper = vtki.new("PolyDataMapper")
 505        self.dataset = vtki.vtkPolyData()
 506        
 507        # Create weakref so actor can access this object (eg to pick/remove):
 508        self.actor.retrieve_object = weak_ref_to(self)
 509
 510        try:
 511            self.properties.RenderPointsAsSpheresOn()
 512        except AttributeError:
 513            pass
 514
 515        if inputobj is None:  ####################
 516            return
 517        ##########################################
 518
 519        self.name = "Points"
 520
 521        ######
 522        if isinstance(inputobj, vtki.vtkActor):
 523            self.dataset.DeepCopy(inputobj.GetMapper().GetInput())
 524            pr = vtki.vtkProperty()
 525            pr.DeepCopy(inputobj.GetProperty())
 526            self.actor.SetProperty(pr)
 527            self.properties = pr
 528            self.mapper.SetScalarVisibility(inputobj.GetMapper().GetScalarVisibility())
 529
 530        elif isinstance(inputobj, vtki.vtkPolyData):
 531            self.dataset = inputobj
 532            if self.dataset.GetNumberOfCells() == 0:
 533                carr = vtki.vtkCellArray()
 534                for i in range(self.dataset.GetNumberOfPoints()):
 535                    carr.InsertNextCell(1)
 536                    carr.InsertCellPoint(i)
 537                self.dataset.SetVerts(carr)
 538
 539        elif isinstance(inputobj, Points):
 540            self.dataset = inputobj.dataset
 541            self.copy_properties_from(inputobj)
 542
 543        elif utils.is_sequence(inputobj):  # passing point coords
 544            self.dataset = utils.buildPolyData(utils.make3d(inputobj))
 545
 546        elif isinstance(inputobj, str):
 547            verts = vedo.file_io.load(inputobj)
 548            self.filename = inputobj
 549            self.dataset = verts.dataset
 550
 551        elif "meshlib" in str(type(inputobj)):
 552            from meshlib import mrmeshnumpy as mn
 553            self.dataset = utils.buildPolyData(mn.toNumpyArray(inputobj.points))
 554
 555        else:
 556            # try to extract the points from a generic VTK input data object
 557            if hasattr(inputobj, "dataset"):
 558                inputobj = inputobj.dataset
 559            try:
 560                vvpts = inputobj.GetPoints()
 561                self.dataset = vtki.vtkPolyData()
 562                self.dataset.SetPoints(vvpts)
 563                for i in range(inputobj.GetPointData().GetNumberOfArrays()):
 564                    arr = inputobj.GetPointData().GetArray(i)
 565                    self.dataset.GetPointData().AddArray(arr)
 566                carr = vtki.vtkCellArray()
 567                for i in range(self.dataset.GetNumberOfPoints()):
 568                    carr.InsertNextCell(1)
 569                    carr.InsertCellPoint(i)
 570                self.dataset.SetVerts(carr)
 571            except:
 572                vedo.logger.error(f"cannot build Points from type {type(inputobj)}")
 573                raise RuntimeError()
 574
 575        self.actor.SetMapper(self.mapper)
 576        self.mapper.SetInputData(self.dataset)
 577
 578        self.properties.SetColor(colors.get_color(c))
 579        self.properties.SetOpacity(alpha)
 580        self.properties.SetRepresentationToPoints()
 581        self.properties.SetPointSize(r)
 582        self.properties.LightingOff()
 583
 584        self.pipeline = utils.OperationNode(
 585            self, parents=[], comment=f"#pts {self.dataset.GetNumberOfPoints()}"
 586        )
 587
 588    def _update(self, polydata, reset_locators=True) -> Self:
 589        """Overwrite the polygonal dataset with a new vtkPolyData."""
 590        self.dataset = polydata
 591        self.mapper.SetInputData(self.dataset)
 592        self.mapper.Modified()
 593        if reset_locators:
 594            self.point_locator = None
 595            self.line_locator = None
 596            self.cell_locator = None
 597        return self
 598
 599    def __str__(self):
 600        """Print a description of the Points/Mesh."""
 601        module = self.__class__.__module__
 602        name = self.__class__.__name__
 603        out = vedo.printc(
 604            f"{module}.{name} at ({hex(self.memory_address())})".ljust(75),
 605            c="g", bold=True, invert=True, return_string=True,
 606        )
 607        out += "\x1b[0m\x1b[32;1m"
 608
 609        if self.name:
 610            out += "name".ljust(14) + ": " + self.name
 611            if "legend" in self.info.keys() and self.info["legend"]:
 612                out+= f", legend='{self.info['legend']}'"
 613            out += "\n"
 614 
 615        if self.filename:
 616            out+= "file name".ljust(14) + ": " + self.filename + "\n"
 617
 618        if not self.mapper.GetScalarVisibility():
 619            col = utils.precision(self.properties.GetColor(), 3)
 620            cname = vedo.colors.get_color_name(self.properties.GetColor())
 621            out+= "color".ljust(14) + ": " + cname 
 622            out+= f", rgb={col}, alpha={self.properties.GetOpacity()}\n"
 623            if self.actor.GetBackfaceProperty():
 624                bcol = self.actor.GetBackfaceProperty().GetDiffuseColor()
 625                cname = vedo.colors.get_color_name(bcol)
 626                out+= "backface color".ljust(14) + ": " 
 627                out+= f"{cname}, rgb={utils.precision(bcol,3)}\n"
 628
 629        npt = self.dataset.GetNumberOfPoints()
 630        npo, nln = self.dataset.GetNumberOfPolys(), self.dataset.GetNumberOfLines()
 631        out+= "elements".ljust(14) + f": vertices={npt:,} polygons={npo:,} lines={nln:,}"
 632        if self.dataset.GetNumberOfStrips():
 633            out+= f", strips={self.dataset.GetNumberOfStrips():,}"
 634        out+= "\n"
 635        if self.dataset.GetNumberOfPieces() > 1:
 636            out+= "pieces".ljust(14) + ": " + str(self.dataset.GetNumberOfPieces()) + "\n"
 637
 638        out+= "position".ljust(14) + ": " + f"{utils.precision(self.pos(), 6)}\n"
 639        try:
 640            sc = self.transform.get_scale()
 641            out+= "scaling".ljust(14)  + ": "
 642            out+= utils.precision(sc, 6) + "\n"
 643        except AttributeError:
 644            pass
 645
 646        if self.npoints:
 647            out+="size".ljust(14)+ ": average=" + utils.precision(self.average_size(),6)
 648            out+=", diagonal="+ utils.precision(self.diagonal_size(), 6)+ "\n"
 649            out+="center of mass".ljust(14) + ": " + utils.precision(self.center_of_mass(),6)+"\n"
 650
 651        bnds = self.bounds()
 652        bx1, bx2 = utils.precision(bnds[0], 3), utils.precision(bnds[1], 3)
 653        by1, by2 = utils.precision(bnds[2], 3), utils.precision(bnds[3], 3)
 654        bz1, bz2 = utils.precision(bnds[4], 3), utils.precision(bnds[5], 3)
 655        out+= "bounds".ljust(14) + ":"
 656        out+= " x=(" + bx1 + ", " + bx2 + "),"
 657        out+= " y=(" + by1 + ", " + by2 + "),"
 658        out+= " z=(" + bz1 + ", " + bz2 + ")\n"
 659
 660        for key in self.pointdata.keys():
 661            arr = self.pointdata[key]
 662            dim = arr.shape[1] if arr.ndim > 1 else 1
 663            mark_active = "pointdata"
 664            a_scalars = self.dataset.GetPointData().GetScalars()
 665            a_vectors = self.dataset.GetPointData().GetVectors()
 666            a_tensors = self.dataset.GetPointData().GetTensors()
 667            if   a_scalars and a_scalars.GetName() == key:
 668                mark_active += " *"
 669            elif a_vectors and a_vectors.GetName() == key:
 670                mark_active += " **"
 671            elif a_tensors and a_tensors.GetName() == key:
 672                mark_active += " ***"
 673            out += mark_active.ljust(14) + f': "{key}" ({arr.dtype}), dim={dim}'
 674            if dim == 1 and len(arr):
 675                rng = utils.precision(arr.min(), 3) + ", " + utils.precision(arr.max(), 3)
 676                out += f", range=({rng})\n"
 677            else:
 678                out += "\n"
 679
 680        for key in self.celldata.keys():
 681            arr = self.celldata[key]
 682            dim = arr.shape[1] if arr.ndim > 1 else 1
 683            mark_active = "celldata"
 684            a_scalars = self.dataset.GetCellData().GetScalars()
 685            a_vectors = self.dataset.GetCellData().GetVectors()
 686            a_tensors = self.dataset.GetCellData().GetTensors()
 687            if   a_scalars and a_scalars.GetName() == key:
 688                mark_active += " *"
 689            elif a_vectors and a_vectors.GetName() == key:
 690                mark_active += " **"
 691            elif a_tensors and a_tensors.GetName() == key:
 692                mark_active += " ***"
 693            out += mark_active.ljust(14) + f': "{key}" ({arr.dtype}), dim={dim}'
 694            if dim == 1 and len(arr):
 695                rng = utils.precision(arr.min(), 3) + ", " + utils.precision(arr.max(), 3)
 696                out += f", range=({rng})\n"
 697            else:
 698                out += "\n"
 699
 700        for key in self.metadata.keys():
 701            arr = self.metadata[key]
 702            if len(arr) > 3:
 703                out+= "metadata".ljust(14) + ": " + f'"{key}" ({len(arr)} values)\n'
 704            else:
 705                out+= "metadata".ljust(14) + ": " + f'"{key}" = {arr}\n'
 706
 707        if self.picked3d is not None:
 708            idp = self.closest_point(self.picked3d, return_point_id=True)
 709            idc = self.closest_point(self.picked3d, return_cell_id=True)
 710            out+= "clicked point".ljust(14) + ": " + utils.precision(self.picked3d, 6)
 711            out+= f", pointID={idp}, cellID={idc}\n"
 712
 713        return out.rstrip() + "\x1b[0m"
 714
 715    def _repr_html_(self):
 716        """
 717        HTML representation of the Point cloud object for Jupyter Notebooks.
 718
 719        Returns:
 720            HTML text with the image and some properties.
 721        """
 722        import io
 723        import base64
 724        from PIL import Image
 725
 726        library_name = "vedo.pointcloud.Points"
 727        help_url = "https://vedo.embl.es/docs/vedo/pointcloud.html#Points"
 728
 729        arr = self.thumbnail()
 730        im = Image.fromarray(arr)
 731        buffered = io.BytesIO()
 732        im.save(buffered, format="PNG", quality=100)
 733        encoded = base64.b64encode(buffered.getvalue()).decode("utf-8")
 734        url = "data:image/png;base64," + encoded
 735        image = f"<img src='{url}'></img>"
 736
 737        bounds = "<br/>".join(
 738            [
 739                utils.precision(min_x, 4) + " ... " + utils.precision(max_x, 4)
 740                for min_x, max_x in zip(self.bounds()[::2], self.bounds()[1::2])
 741            ]
 742        )
 743        average_size = "{size:.3f}".format(size=self.average_size())
 744
 745        help_text = ""
 746        if self.name:
 747            help_text += f"<b> {self.name}: &nbsp&nbsp</b>"
 748        help_text += '<b><a href="' + help_url + '" target="_blank">' + library_name + "</a></b>"
 749        if self.filename:
 750            dots = ""
 751            if len(self.filename) > 30:
 752                dots = "..."
 753            help_text += f"<br/><code><i>({dots}{self.filename[-30:]})</i></code>"
 754
 755        pdata = ""
 756        if self.dataset.GetPointData().GetScalars():
 757            if self.dataset.GetPointData().GetScalars().GetName():
 758                name = self.dataset.GetPointData().GetScalars().GetName()
 759                pdata = "<tr><td><b> point data array </b></td><td>" + name + "</td></tr>"
 760
 761        cdata = ""
 762        if self.dataset.GetCellData().GetScalars():
 763            if self.dataset.GetCellData().GetScalars().GetName():
 764                name = self.dataset.GetCellData().GetScalars().GetName()
 765                cdata = "<tr><td><b> cell data array </b></td><td>" + name + "</td></tr>"
 766
 767        allt = [
 768            "<table>",
 769            "<tr>",
 770            "<td>",
 771            image,
 772            "</td>",
 773            "<td style='text-align: center; vertical-align: center;'><br/>",
 774            help_text,
 775            "<table>",
 776            "<tr><td><b> bounds </b> <br/> (x/y/z) </td><td>" + str(bounds) + "</td></tr>",
 777            "<tr><td><b> center of mass </b></td><td>"
 778            + utils.precision(self.center_of_mass(), 3)
 779            + "</td></tr>",
 780            "<tr><td><b> average size </b></td><td>" + str(average_size) + "</td></tr>",
 781            "<tr><td><b> nr. points </b></td><td>" + str(self.npoints) + "</td></tr>",
 782            pdata,
 783            cdata,
 784            "</table>",
 785            "</table>",
 786        ]
 787        return "\n".join(allt)
 788
 789    ##################################################################################
 790    def __add__(self, meshs):
 791        """
 792        Add two meshes or a list of meshes together to form an `Assembly` object.
 793        """
 794        if isinstance(meshs, list):
 795            alist = [self]
 796            for l in meshs:
 797                if isinstance(l, vedo.Assembly):
 798                    alist += l.unpack()
 799                else:
 800                    alist += l
 801            return vedo.assembly.Assembly(alist)
 802
 803        if isinstance(meshs, vedo.Assembly):
 804            return meshs + self  # use Assembly.__add__
 805
 806        return vedo.assembly.Assembly([self, meshs])
 807
 808    def polydata(self, **kwargs):
 809        """
 810        Obsolete. Use property `.dataset` instead.
 811        Returns the underlying `vtkPolyData` object.
 812        """
 813        colors.printc(
 814            "WARNING: call to .polydata() is obsolete, use property .dataset instead.",
 815            c="y")
 816        return self.dataset
 817
 818    def __copy__(self):
 819        return self.clone(deep=False)
 820
 821    def __deepcopy__(self, memo):
 822        return self.clone(deep=memo)
 823    
 824    def copy(self, deep=True) -> Self:
 825        """Return a copy of the object. Alias of `clone()`."""
 826        return self.clone(deep=deep)
 827
 828    def clone(self, deep=True) -> Self:
 829        """
 830        Clone a `PointCloud` or `Mesh` object to make an exact copy of it.
 831        Alias of `copy()`.
 832
 833        Arguments:
 834            deep : (bool)
 835                if False return a shallow copy of the mesh without copying the points array.
 836
 837        Examples:
 838            - [mirror.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/mirror.py)
 839
 840               ![](https://vedo.embl.es/images/basic/mirror.png)
 841        """
 842        poly = vtki.vtkPolyData()
 843        if deep or isinstance(deep, dict): # if a memo object is passed this checks as True
 844            poly.DeepCopy(self.dataset)
 845        else:
 846            poly.ShallowCopy(self.dataset)
 847
 848        if isinstance(self, vedo.Mesh):
 849            cloned = vedo.Mesh(poly)
 850        else:
 851            cloned = Points(poly)
 852        # print([self], self.__class__)
 853        # cloned = self.__class__(poly)
 854
 855        cloned.transform = self.transform.clone()
 856
 857        cloned.copy_properties_from(self)
 858
 859        cloned.name = str(self.name)
 860        cloned.filename = str(self.filename)
 861        cloned.info = dict(self.info)
 862        cloned.pipeline = utils.OperationNode("clone", parents=[self], shape="diamond", c="#edede9")
 863
 864        if isinstance(deep, dict):
 865            deep[id(self)] = cloned
 866
 867        return cloned
 868
 869    def compute_normals_with_pca(self, n=20, orientation_point=None, invert=False) -> Self:
 870        """
 871        Generate point normals using PCA (principal component analysis).
 872        This algorithm estimates a local tangent plane around each sample point p
 873        by considering a small neighborhood of points around p, and fitting a plane
 874        to the neighborhood (via PCA).
 875
 876        Arguments:
 877            n : (int)
 878                neighborhood size to calculate the normal
 879            orientation_point : (list)
 880                adjust the +/- sign of the normals so that
 881                the normals all point towards a specified point. If None, perform a traversal
 882                of the point cloud and flip neighboring normals so that they are mutually consistent.
 883            invert : (bool)
 884                flip all normals
 885        """
 886        poly = self.dataset
 887        pcan = vtki.new("PCANormalEstimation")
 888        pcan.SetInputData(poly)
 889        pcan.SetSampleSize(n)
 890
 891        if orientation_point is not None:
 892            pcan.SetNormalOrientationToPoint()
 893            pcan.SetOrientationPoint(orientation_point)
 894        else:
 895            pcan.SetNormalOrientationToGraphTraversal()
 896
 897        if invert:
 898            pcan.FlipNormalsOn()
 899        pcan.Update()
 900
 901        varr = pcan.GetOutput().GetPointData().GetNormals()
 902        varr.SetName("Normals")
 903        self.dataset.GetPointData().SetNormals(varr)
 904        self.dataset.GetPointData().Modified()
 905        return self
 906
 907    def compute_acoplanarity(self, n=25, radius=None, on="points") -> Self:
 908        """
 909        Compute acoplanarity which is a measure of how much a local region of the mesh
 910        differs from a plane.
 911        
 912        The information is stored in a `pointdata` or `celldata` array with name 'Acoplanarity'.
 913        
 914        Either `n` (number of neighbour points) or `radius` (radius of local search) can be specified.
 915        If a radius value is given and not enough points fall inside it, then a -1 is stored.
 916
 917        Example:
 918            ```python
 919            from vedo import *
 920            msh = ParametricShape('RandomHills')
 921            msh.compute_acoplanarity(radius=0.1, on='cells')
 922            msh.cmap("coolwarm", on='cells').add_scalarbar()
 923            msh.show(axes=1).close()
 924            ```
 925            ![](https://vedo.embl.es/images/feats/acoplanarity.jpg)
 926        """
 927        acoplanarities = []
 928        if "point" in on:
 929            pts = self.coordinates
 930        elif "cell" in on:
 931            pts = self.cell_centers().coordinates
 932        else:
 933            raise ValueError(f"In compute_acoplanarity() set on to either 'cells' or 'points', not {on}")
 934
 935        for p in utils.progressbar(pts, delay=5, width=15, title=f"{on} acoplanarity"):
 936            if n:
 937                data = self.closest_point(p, n=n)
 938                npts = n
 939            elif radius:
 940                data = self.closest_point(p, radius=radius)
 941                npts = len(data)
 942
 943            try:
 944                center = data.mean(axis=0)
 945                res = np.linalg.svd(data - center)
 946                acoplanarities.append(res[1][2] / npts)
 947            except:
 948                acoplanarities.append(-1.0)
 949
 950        if "point" in on:
 951            self.pointdata["Acoplanarity"] = np.array(acoplanarities, dtype=float)
 952        else:
 953            self.celldata["Acoplanarity"] = np.array(acoplanarities, dtype=float)
 954        return self
 955
 956    def distance_to(self, pcloud, signed=False, invert=False, name="Distance") -> np.ndarray:
 957        """
 958        Computes the distance from one point cloud or mesh to another point cloud or mesh.
 959        This new `pointdata` array is saved with default name "Distance".
 960
 961        Keywords `signed` and `invert` are used to compute signed distance,
 962        but the mesh in that case must have polygonal faces (not a simple point cloud),
 963        and normals must also be computed.
 964
 965        Examples:
 966            - [distance2mesh.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/distance2mesh.py)
 967
 968                ![](https://vedo.embl.es/images/basic/distance2mesh.png)
 969        """
 970        if pcloud.dataset.GetNumberOfPolys():
 971
 972            poly1 = self.dataset
 973            poly2 = pcloud.dataset
 974            df = vtki.new("DistancePolyDataFilter")
 975            df.ComputeSecondDistanceOff()
 976            df.SetInputData(0, poly1)
 977            df.SetInputData(1, poly2)
 978            df.SetSignedDistance(signed)
 979            df.SetNegateDistance(invert)
 980            df.Update()
 981            scals = df.GetOutput().GetPointData().GetScalars()
 982            dists = utils.vtk2numpy(scals)
 983
 984        else:  # has no polygons
 985
 986            if signed:
 987                vedo.logger.warning("distance_to() called with signed=True but input object has no polygons")
 988
 989            if not pcloud.point_locator:
 990                pcloud.point_locator = vtki.new("PointLocator")
 991                pcloud.point_locator.SetDataSet(pcloud.dataset)
 992                pcloud.point_locator.BuildLocator()
 993
 994            ids = []
 995            ps1 = self.coordinates
 996            ps2 = pcloud.coordinates
 997            for p in ps1:
 998                pid = pcloud.point_locator.FindClosestPoint(p)
 999                ids.append(pid)
1000
1001            deltas = ps2[ids] - ps1
1002            dists = np.linalg.norm(deltas, axis=1).astype(np.float32)
1003            scals = utils.numpy2vtk(dists)
1004
1005        scals.SetName(name)
1006        self.dataset.GetPointData().AddArray(scals)
1007        self.dataset.GetPointData().SetActiveScalars(scals.GetName())
1008        rng = scals.GetRange()
1009        self.mapper.SetScalarRange(rng[0], rng[1])
1010        self.mapper.ScalarVisibilityOn()
1011
1012        self.pipeline = utils.OperationNode(
1013            "distance_to",
1014            parents=[self, pcloud],
1015            shape="cylinder",
1016            comment=f"#pts {self.dataset.GetNumberOfPoints()}",
1017        )
1018        return dists
1019
1020    def clean(self) -> Self:
1021        """Clean pointcloud or mesh by removing coincident points."""
1022        cpd = vtki.new("CleanPolyData")
1023        cpd.PointMergingOn()
1024        cpd.ConvertLinesToPointsOff()
1025        cpd.ConvertPolysToLinesOff()
1026        cpd.ConvertStripsToPolysOff()
1027        cpd.SetInputData(self.dataset)
1028        cpd.Update()
1029        self._update(cpd.GetOutput())
1030        self.pipeline = utils.OperationNode(
1031            "clean", parents=[self], comment=f"#pts {self.dataset.GetNumberOfPoints()}"
1032        )
1033        return self
1034
1035    def subsample(self, fraction: float, absolute=False) -> Self:
1036        """
1037        Subsample a point cloud by requiring that the points
1038        or vertices are far apart at least by the specified fraction of the object size.
1039        If a Mesh is passed the polygonal faces are not removed
1040        but holes can appear as their vertices are removed.
1041
1042        Examples:
1043            - [moving_least_squares1D.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/moving_least_squares1D.py)
1044
1045                ![](https://vedo.embl.es/images/advanced/moving_least_squares1D.png)
1046
1047            - [recosurface.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/recosurface.py)
1048
1049                ![](https://vedo.embl.es/images/advanced/recosurface.png)
1050        """
1051        if not absolute:
1052            if fraction > 1:
1053                vedo.logger.warning(
1054                    f"subsample(fraction=...), fraction must be < 1, but is {fraction}"
1055                )
1056            if fraction <= 0:
1057                return self
1058
1059        cpd = vtki.new("CleanPolyData")
1060        cpd.PointMergingOn()
1061        cpd.ConvertLinesToPointsOn()
1062        cpd.ConvertPolysToLinesOn()
1063        cpd.ConvertStripsToPolysOn()
1064        cpd.SetInputData(self.dataset)
1065        if absolute:
1066            cpd.SetTolerance(fraction / self.diagonal_size())
1067            # cpd.SetToleranceIsAbsolute(absolute)
1068        else:
1069            cpd.SetTolerance(fraction)
1070        cpd.Update()
1071
1072        ps = 2
1073        if self.properties.GetRepresentation() == 0:
1074            ps = self.properties.GetPointSize()
1075
1076        self._update(cpd.GetOutput())
1077        self.ps(ps)
1078
1079        self.pipeline = utils.OperationNode(
1080            "subsample", parents=[self], comment=f"#pts {self.dataset.GetNumberOfPoints()}"
1081        )
1082        return self
1083
1084    def threshold(self, scalars: str, above=None, below=None, on="points") -> Self:
1085        """
1086        Extracts cells where scalar value satisfies threshold criterion.
1087
1088        Arguments:
1089            scalars : (str)
1090                name of the scalars array.
1091            above : (float)
1092                minimum value of the scalar
1093            below : (float)
1094                maximum value of the scalar
1095            on : (str)
1096                if 'cells' assume array of scalars refers to cell data.
1097
1098        Examples:
1099            - [mesh_threshold.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/mesh_threshold.py)
1100        """
1101        thres = vtki.new("Threshold")
1102        thres.SetInputData(self.dataset)
1103
1104        if on.startswith("c"):
1105            asso = vtki.vtkDataObject.FIELD_ASSOCIATION_CELLS
1106        else:
1107            asso = vtki.vtkDataObject.FIELD_ASSOCIATION_POINTS
1108
1109        thres.SetInputArrayToProcess(0, 0, 0, asso, scalars)
1110
1111        if above is None and below is not None:
1112            try:  # vtk 9.2
1113                thres.ThresholdByLower(below)
1114            except AttributeError:  # vtk 9.3
1115                thres.SetUpperThreshold(below)
1116
1117        elif below is None and above is not None:
1118            try:
1119                thres.ThresholdByUpper(above)
1120            except AttributeError:
1121                thres.SetLowerThreshold(above)
1122        else:
1123            try:
1124                thres.ThresholdBetween(above, below)
1125            except AttributeError:
1126                thres.SetUpperThreshold(below)
1127                thres.SetLowerThreshold(above)
1128
1129        thres.Update()
1130
1131        gf = vtki.new("GeometryFilter")
1132        gf.SetInputData(thres.GetOutput())
1133        gf.Update()
1134        self._update(gf.GetOutput())
1135        self.pipeline = utils.OperationNode("threshold", parents=[self])
1136        return self
1137
1138    def quantize(self, value: float) -> Self:
1139        """
1140        The user should input a value and all {x,y,z} coordinates
1141        will be quantized to that absolute grain size.
1142        """
1143        qp = vtki.new("QuantizePolyDataPoints")
1144        qp.SetInputData(self.dataset)
1145        qp.SetQFactor(value)
1146        qp.Update()
1147        self._update(qp.GetOutput())
1148        self.pipeline = utils.OperationNode("quantize", parents=[self])
1149        return self
1150
1151    @property
1152    def vertex_normals(self) -> np.ndarray:
1153        """
1154        Retrieve vertex normals as a numpy array. Same as `point_normals`.
1155        Check out also `compute_normals()` and `compute_normals_with_pca()`.
1156        """
1157        vtknormals = self.dataset.GetPointData().GetNormals()
1158        return utils.vtk2numpy(vtknormals)
1159
1160    @property
1161    def point_normals(self) -> np.ndarray:
1162        """
1163        Retrieve vertex normals as a numpy array. Same as `vertex_normals`.
1164        Check out also `compute_normals()` and `compute_normals_with_pca()`.
1165        """
1166        vtknormals = self.dataset.GetPointData().GetNormals()
1167        return utils.vtk2numpy(vtknormals)
1168
1169    def align_to(self, target, iters=100, rigid=False, invert=False, use_centroids=False) -> Self:
1170        """
1171        Aligned to target mesh through the `Iterative Closest Point` algorithm.
1172
1173        The core of the algorithm is to match each vertex in one surface with
1174        the closest surface point on the other, then apply the transformation
1175        that modify one surface to best match the other (in the least-square sense).
1176
1177        Arguments:
1178            rigid : (bool)
1179                if True do not allow scaling
1180            invert : (bool)
1181                if True start by aligning the target to the source but
1182                invert the transformation finally. Useful when the target is smaller
1183                than the source.
1184            use_centroids : (bool)
1185                start by matching the centroids of the two objects.
1186
1187        Examples:
1188            - [align1.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/align1.py)
1189
1190                ![](https://vedo.embl.es/images/basic/align1.png)
1191
1192            - [align2.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/align2.py)
1193
1194                ![](https://vedo.embl.es/images/basic/align2.png)
1195        """
1196        icp = vtki.new("IterativeClosestPointTransform")
1197        icp.SetSource(self.dataset)
1198        icp.SetTarget(target.dataset)
1199        if invert:
1200            icp.Inverse()
1201        icp.SetMaximumNumberOfIterations(iters)
1202        if rigid:
1203            icp.GetLandmarkTransform().SetModeToRigidBody()
1204        icp.SetStartByMatchingCentroids(use_centroids)
1205        icp.Update()
1206
1207        self.apply_transform(icp.GetMatrix())
1208
1209        self.pipeline = utils.OperationNode(
1210            "align_to", parents=[self, target], comment=f"rigid = {rigid}"
1211        )
1212        return self
1213
1214    def align_to_bounding_box(self, msh, rigid=False) -> Self:
1215        """
1216        Align the current object's bounding box to the bounding box
1217        of the input object.
1218
1219        Use `rigid=True` to disable scaling.
1220
1221        Example:
1222            [align6.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/align6.py)
1223        """
1224        lmt = vtki.vtkLandmarkTransform()
1225        ss = vtki.vtkPoints()
1226        xss0, xss1, yss0, yss1, zss0, zss1 = self.bounds()
1227        for p in [
1228            [xss0, yss0, zss0],
1229            [xss1, yss0, zss0],
1230            [xss1, yss1, zss0],
1231            [xss0, yss1, zss0],
1232            [xss0, yss0, zss1],
1233            [xss1, yss0, zss1],
1234            [xss1, yss1, zss1],
1235            [xss0, yss1, zss1],
1236        ]:
1237            ss.InsertNextPoint(p)
1238        st = vtki.vtkPoints()
1239        xst0, xst1, yst0, yst1, zst0, zst1 = msh.bounds()
1240        for p in [
1241            [xst0, yst0, zst0],
1242            [xst1, yst0, zst0],
1243            [xst1, yst1, zst0],
1244            [xst0, yst1, zst0],
1245            [xst0, yst0, zst1],
1246            [xst1, yst0, zst1],
1247            [xst1, yst1, zst1],
1248            [xst0, yst1, zst1],
1249        ]:
1250            st.InsertNextPoint(p)
1251
1252        lmt.SetSourceLandmarks(ss)
1253        lmt.SetTargetLandmarks(st)
1254        lmt.SetModeToAffine()
1255        if rigid:
1256            lmt.SetModeToRigidBody()
1257        lmt.Update()
1258
1259        LT = LinearTransform(lmt)
1260        self.apply_transform(LT)
1261        return self
1262
1263    def align_with_landmarks(
1264        self,
1265        source_landmarks,
1266        target_landmarks,
1267        rigid=False,
1268        affine=False,
1269        least_squares=False,
1270    ) -> Self:
1271        """
1272        Transform mesh orientation and position based on a set of landmarks points.
1273        The algorithm finds the best matching of source points to target points
1274        in the mean least square sense, in one single step.
1275
1276        If `affine` is True the x, y and z axes can scale independently but stay collinear.
1277        With least_squares they can vary orientation.
1278
1279        Examples:
1280            - [align5.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/align5.py)
1281
1282                ![](https://vedo.embl.es/images/basic/align5.png)
1283        """
1284
1285        if utils.is_sequence(source_landmarks):
1286            ss = vtki.vtkPoints()
1287            for p in source_landmarks:
1288                ss.InsertNextPoint(p)
1289        else:
1290            ss = source_landmarks.dataset.GetPoints()
1291            if least_squares:
1292                source_landmarks = source_landmarks.coordinates
1293
1294        if utils.is_sequence(target_landmarks):
1295            st = vtki.vtkPoints()
1296            for p in target_landmarks:
1297                st.InsertNextPoint(p)
1298        else:
1299            st = target_landmarks.GetPoints()
1300            if least_squares:
1301                target_landmarks = target_landmarks.coordinates
1302
1303        if ss.GetNumberOfPoints() != st.GetNumberOfPoints():
1304            n1 = ss.GetNumberOfPoints()
1305            n2 = st.GetNumberOfPoints()
1306            vedo.logger.error(f"source and target have different nr of points {n1} vs {n2}")
1307            raise RuntimeError()
1308
1309        if int(rigid) + int(affine) + int(least_squares) > 1:
1310            vedo.logger.error(
1311                "only one of rigid, affine, least_squares can be True at a time"
1312            )
1313            raise RuntimeError()
1314
1315        lmt = vtki.vtkLandmarkTransform()
1316        lmt.SetSourceLandmarks(ss)
1317        lmt.SetTargetLandmarks(st)
1318        lmt.SetModeToSimilarity()
1319
1320        if rigid:
1321            lmt.SetModeToRigidBody()
1322            lmt.Update()
1323
1324        elif affine:
1325            lmt.SetModeToAffine()
1326            lmt.Update()
1327
1328        elif least_squares:
1329            cms = source_landmarks.mean(axis=0)
1330            cmt = target_landmarks.mean(axis=0)
1331            m = np.linalg.lstsq(source_landmarks - cms, target_landmarks - cmt, rcond=None)[0]
1332            M = vtki.vtkMatrix4x4()
1333            for i in range(3):
1334                for j in range(3):
1335                    M.SetElement(j, i, m[i][j])
1336            lmt = vtki.vtkTransform()
1337            lmt.Translate(cmt)
1338            lmt.Concatenate(M)
1339            lmt.Translate(-cms)
1340
1341        else:
1342            lmt.Update()
1343
1344        self.apply_transform(lmt)
1345        self.pipeline = utils.OperationNode("transform_with_landmarks", parents=[self])
1346        return self
1347
1348    def normalize(self) -> Self:
1349        """Scale average size to unit. The scaling is performed around the center of mass."""
1350        coords = self.coordinates
1351        if not coords.shape[0]:
1352            return self
1353        cm = np.mean(coords, axis=0)
1354        pts = coords - cm
1355        xyz2 = np.sum(pts * pts, axis=0)
1356        scale = 1 / np.sqrt(np.sum(xyz2) / len(pts))
1357        self.scale(scale, origin=cm)
1358        self.pipeline = utils.OperationNode("normalize", parents=[self])
1359        return self
1360
1361    def mirror(self, axis="x", origin=True) -> Self:
1362        """
1363        Mirror reflect along one of the cartesian axes
1364
1365        Arguments:
1366            axis : (str)
1367                axis to use for mirroring, must be set to `x, y, z`.
1368                Or any combination of those.
1369            origin : (list)
1370                use this point as the origin of the mirroring transformation.
1371
1372        Examples:
1373            - [mirror.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/mirror.py)
1374
1375                ![](https://vedo.embl.es/images/basic/mirror.png)
1376        """
1377        sx, sy, sz = 1, 1, 1
1378        if "x" in axis.lower(): sx = -1
1379        if "y" in axis.lower(): sy = -1
1380        if "z" in axis.lower(): sz = -1
1381
1382        self.scale([sx, sy, sz], origin=origin)
1383
1384        self.pipeline = utils.OperationNode(
1385            "mirror", comment=f"axis = {axis}", parents=[self])
1386
1387        if sx * sy * sz < 0:
1388            if hasattr(self, "reverse"):
1389                self.reverse()
1390        return self
1391
1392    def flip_normals(self) -> Self:
1393        """Flip all normals orientation."""
1394        rs = vtki.new("ReverseSense")
1395        rs.SetInputData(self.dataset)
1396        rs.ReverseCellsOff()
1397        rs.ReverseNormalsOn()
1398        rs.Update()
1399        self._update(rs.GetOutput())
1400        self.pipeline = utils.OperationNode("flip_normals", parents=[self])
1401        return self
1402
1403    def add_gaussian_noise(self, sigma=1.0) -> Self:
1404        """
1405        Add gaussian noise to point positions.
1406        An extra array is added named "GaussianNoise" with the displacements.
1407
1408        Arguments:
1409            sigma : (float)
1410                nr. of standard deviations, expressed in percent of the diagonal size of mesh.
1411                Can also be a list `[sigma_x, sigma_y, sigma_z]`.
1412
1413        Example:
1414            ```python
1415            from vedo import Sphere
1416            Sphere().add_gaussian_noise(1.0).point_size(8).show().close()
1417            ```
1418        """
1419        sz = self.diagonal_size()
1420        pts = self.coordinates
1421        n = len(pts)
1422        ns = (np.random.randn(n, 3) * sigma) * (sz / 100)
1423        vpts = vtki.vtkPoints()
1424        vpts.SetNumberOfPoints(n)
1425        vpts.SetData(utils.numpy2vtk(pts + ns, dtype=np.float32))
1426        self.dataset.SetPoints(vpts)
1427        self.dataset.GetPoints().Modified()
1428        self.pointdata["GaussianNoise"] = -ns
1429        self.pipeline = utils.OperationNode(
1430            "gaussian_noise", parents=[self], shape="egg", comment=f"sigma = {sigma}"
1431        )
1432        return self
1433
1434    def closest_point(
1435        self, pt, n=1, radius=None, return_point_id=False, return_cell_id=False
1436    ) -> Union[List[int], int, np.ndarray]:
1437        """
1438        Find the closest point(s) on a mesh given from the input point `pt`.
1439
1440        Arguments:
1441            n : (int)
1442                if greater than 1, return a list of n ordered closest points
1443            radius : (float)
1444                if given, get all points within that radius. Then n is ignored.
1445            return_point_id : (bool)
1446                return point ID instead of coordinates
1447            return_cell_id : (bool)
1448                return cell ID in which the closest point sits
1449
1450        Examples:
1451            - [align1.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/align1.py)
1452            - [fitplanes.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/fitplanes.py)
1453            - [quadratic_morphing.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/quadratic_morphing.py)
1454
1455        .. note::
1456            The appropriate tree search locator is built on the fly and cached for speed.
1457
1458            If you want to reset it use `mymesh.point_locator=None`
1459            and / or `mymesh.cell_locator=None`.
1460        """
1461        if len(pt) != 3:
1462            pt = [pt[0], pt[1], 0]
1463
1464        # NB: every time the mesh moves or is warped the locators are set to None
1465        if ((n > 1 or radius) or (n == 1 and return_point_id)) and not return_cell_id:
1466            poly = None
1467            if not self.point_locator:
1468                poly = self.dataset
1469                self.point_locator = vtki.new("StaticPointLocator")
1470                self.point_locator.SetDataSet(poly)
1471                self.point_locator.BuildLocator()
1472
1473            ##########
1474            if radius:
1475                vtklist = vtki.vtkIdList()
1476                self.point_locator.FindPointsWithinRadius(radius, pt, vtklist)
1477            elif n > 1:
1478                vtklist = vtki.vtkIdList()
1479                self.point_locator.FindClosestNPoints(n, pt, vtklist)
1480            else:  # n==1 hence return_point_id==True
1481                ########
1482                return self.point_locator.FindClosestPoint(pt)
1483                ########
1484
1485            if return_point_id:
1486                ########
1487                return utils.vtk2numpy(vtklist)
1488                ########
1489
1490            if not poly:
1491                poly = self.dataset
1492            trgp = []
1493            for i in range(vtklist.GetNumberOfIds()):
1494                trgp_ = [0, 0, 0]
1495                vi = vtklist.GetId(i)
1496                poly.GetPoints().GetPoint(vi, trgp_)
1497                trgp.append(trgp_)
1498            ########
1499            return np.array(trgp)
1500            ########
1501
1502        else:
1503
1504            if not self.cell_locator:
1505                poly = self.dataset
1506
1507                # As per Miquel example with limbs the vtkStaticCellLocator doesnt work !!
1508                # https://discourse.vtk.org/t/vtkstaticcelllocator-problem-vtk9-0-3/7854/4
1509                if vedo.vtk_version[0] >= 9 and vedo.vtk_version[1] > 0:
1510                    self.cell_locator = vtki.new("StaticCellLocator")
1511                else:
1512                    self.cell_locator = vtki.new("CellLocator")
1513
1514                self.cell_locator.SetDataSet(poly)
1515                self.cell_locator.BuildLocator()
1516
1517            if radius is not None:
1518                vedo.printc("Warning: closest_point() with radius is not implemented for cells.", c='r')   
1519 
1520            if n != 1:
1521                vedo.printc("Warning: closest_point() with n>1 is not implemented for cells.", c='r')   
1522 
1523            trgp = [0, 0, 0]
1524            cid = vtki.mutable(0)
1525            dist2 = vtki.mutable(0)
1526            subid = vtki.mutable(0)
1527            self.cell_locator.FindClosestPoint(pt, trgp, cid, subid, dist2)
1528
1529            if return_cell_id:
1530                return int(cid)
1531
1532            return np.array(trgp)
1533
1534    def auto_distance(self) -> np.ndarray:
1535        """
1536        Calculate the distance to the closest point in the same cloud of points.
1537        The output is stored in a new pointdata array called "AutoDistance",
1538        and it is also returned by the function.
1539        """
1540        points = self.coordinates
1541        if not self.point_locator:
1542            self.point_locator = vtki.new("StaticPointLocator")
1543            self.point_locator.SetDataSet(self.dataset)
1544            self.point_locator.BuildLocator()
1545        qs = []
1546        vtklist = vtki.vtkIdList()
1547        vtkpoints = self.dataset.GetPoints()
1548        for p in points:
1549            self.point_locator.FindClosestNPoints(2, p, vtklist)
1550            q = [0, 0, 0]
1551            pid = vtklist.GetId(1)
1552            vtkpoints.GetPoint(pid, q)
1553            qs.append(q)
1554        dists = np.linalg.norm(points - np.array(qs), axis=1)
1555        self.pointdata["AutoDistance"] = dists
1556        return dists
1557
1558    def hausdorff_distance(self, points) -> float:
1559        """
1560        Compute the Hausdorff distance to the input point set.
1561        Returns a single `float`.
1562
1563        Example:
1564            ```python
1565            from vedo import *
1566            t = np.linspace(0, 2*np.pi, 100)
1567            x = 4/3 * sin(t)**3
1568            y = cos(t) - cos(2*t)/3 - cos(3*t)/6 - cos(4*t)/12
1569            pol1 = Line(np.c_[x,y], closed=True).triangulate()
1570            pol2 = Polygon(nsides=5).pos(2,2)
1571            d12 = pol1.distance_to(pol2)
1572            d21 = pol2.distance_to(pol1)
1573            pol1.lw(0).cmap("viridis")
1574            pol2.lw(0).cmap("viridis")
1575            print("distance d12, d21 :", min(d12), min(d21))
1576            print("hausdorff distance:", pol1.hausdorff_distance(pol2))
1577            print("chamfer distance  :", pol1.chamfer_distance(pol2))
1578            show(pol1, pol2, axes=1)
1579            ```
1580            ![](https://vedo.embl.es/images/feats/heart.png)
1581        """
1582        hp = vtki.new("HausdorffDistancePointSetFilter")
1583        hp.SetInputData(0, self.dataset)
1584        hp.SetInputData(1, points.dataset)
1585        hp.SetTargetDistanceMethodToPointToCell()
1586        hp.Update()
1587        return hp.GetHausdorffDistance()
1588
1589    def chamfer_distance(self, pcloud) -> float:
1590        """
1591        Compute the Chamfer distance to the input point set.
1592
1593        Example:
1594            ```python
1595            from vedo import *
1596            cloud1 = np.random.randn(1000, 3)
1597            cloud2 = np.random.randn(1000, 3) + [1, 2, 3]
1598            c1 = Points(cloud1, r=5, c="red")
1599            c2 = Points(cloud2, r=5, c="green")
1600            d = c1.chamfer_distance(c2)
1601            show(f"Chamfer distance = {d}", c1, c2, axes=1).close()
1602            ```
1603        """
1604        # Definition of Chamfer distance may vary, here we use the average
1605        if not pcloud.point_locator:
1606            pcloud.point_locator = vtki.new("PointLocator")
1607            pcloud.point_locator.SetDataSet(pcloud.dataset)
1608            pcloud.point_locator.BuildLocator()
1609        if not self.point_locator:
1610            self.point_locator = vtki.new("PointLocator")
1611            self.point_locator.SetDataSet(self.dataset)
1612            self.point_locator.BuildLocator()
1613
1614        ps1 = self.coordinates
1615        ps2 = pcloud.coordinates
1616
1617        ids12 = []
1618        for p in ps1:
1619            pid12 = pcloud.point_locator.FindClosestPoint(p)
1620            ids12.append(pid12)
1621        deltav = ps2[ids12] - ps1
1622        da = np.mean(np.linalg.norm(deltav, axis=1))
1623
1624        ids21 = []
1625        for p in ps2:
1626            pid21 = self.point_locator.FindClosestPoint(p)
1627            ids21.append(pid21)
1628        deltav = ps1[ids21] - ps2
1629        db = np.mean(np.linalg.norm(deltav, axis=1))
1630        return (da + db) / 2
1631
1632    def remove_outliers(self, radius: float, neighbors=5) -> Self:
1633        """
1634        Remove outliers from a cloud of points within the specified `radius` search.
1635
1636        Arguments:
1637            radius : (float)
1638                Specify the local search radius.
1639            neighbors : (int)
1640                Specify the number of neighbors that a point must have,
1641                within the specified radius, for the point to not be considered isolated.
1642
1643        Examples:
1644            - [clustering.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/clustering.py)
1645
1646                ![](https://vedo.embl.es/images/basic/clustering.png)
1647        """
1648        removal = vtki.new("RadiusOutlierRemoval")
1649        removal.SetInputData(self.dataset)
1650        removal.SetRadius(radius)
1651        removal.SetNumberOfNeighbors(neighbors)
1652        removal.GenerateOutliersOff()
1653        removal.Update()
1654        inputobj = removal.GetOutput()
1655        if inputobj.GetNumberOfCells() == 0:
1656            carr = vtki.vtkCellArray()
1657            for i in range(inputobj.GetNumberOfPoints()):
1658                carr.InsertNextCell(1)
1659                carr.InsertCellPoint(i)
1660            inputobj.SetVerts(carr)
1661        self._update(removal.GetOutput())
1662        self.pipeline = utils.OperationNode("remove_outliers", parents=[self])
1663        return self
1664
1665    def relax_point_positions(
1666            self, 
1667            n=10,
1668            iters=10,
1669            sub_iters=10,
1670            packing_factor=1,
1671            max_step=0,
1672            constraints=(),
1673        ) -> Self:
1674        """
1675        Smooth mesh or points with a 
1676        [Laplacian algorithm](https://vtk.org/doc/nightly/html/classvtkPointSmoothingFilter.html)
1677        variant. This modifies the coordinates of the input points by adjusting their positions
1678        to create a smooth distribution (and thereby form a pleasing packing of the points).
1679        Smoothing is performed by considering the effects of neighboring points on one another
1680        it uses a cubic cutoff function to produce repulsive forces between close points
1681        and attractive forces that are a little further away.
1682        
1683        In general, the larger the neighborhood size, the greater the reduction in high frequency
1684        information. The memory and computational requirements of the algorithm may also
1685        significantly increase.
1686
1687        The algorithm incrementally adjusts the point positions through an iterative process.
1688        Basically points are moved due to the influence of neighboring points. 
1689        
1690        As points move, both the local connectivity and data attributes associated with each point
1691        must be updated. Rather than performing these expensive operations after every iteration,
1692        a number of sub-iterations can be specified. If so, then the neighborhood and attribute
1693        value updates occur only every sub iteration, which can improve performance significantly.
1694        
1695        Arguments:
1696            n : (int)
1697                neighborhood size to calculate the Laplacian.
1698            iters : (int)
1699                number of iterations.
1700            sub_iters : (int)
1701                number of sub-iterations, i.e. the number of times the neighborhood and attribute
1702                value updates occur during each iteration.
1703            packing_factor : (float)
1704                adjust convergence speed.
1705            max_step : (float)
1706                Specify the maximum smoothing step size for each smoothing iteration.
1707                This limits the the distance over which a point can move in each iteration.
1708                As in all iterative methods, the stability of the process is sensitive to this parameter.
1709                In general, small step size and large numbers of iterations are more stable than a larger
1710                step size and a smaller numbers of iterations.
1711            constraints : (dict)
1712                dictionary of constraints.
1713                Point constraints are used to prevent points from moving,
1714                or to move only on a plane. This can prevent shrinking or growing point clouds.
1715                If enabled, a local topological analysis is performed to determine whether a point
1716                should be marked as fixed" i.e., never moves, or the point only moves on a plane,
1717                or the point can move freely.
1718                If all points in the neighborhood surrounding a point are in the cone defined by
1719                `fixed_angle`, then the point is classified as fixed.
1720                If all points in the neighborhood surrounding a point are in the cone defined by
1721                `boundary_angle`, then the point is classified as lying on a plane.
1722                Angles are expressed in degrees.
1723        
1724        Example:
1725            ```py
1726            import numpy as np
1727            from vedo import Points, show
1728            from vedo.pyplot import histogram
1729
1730            vpts1 = Points(np.random.rand(10_000, 3))
1731            dists = vpts1.auto_distance()
1732            h1 = histogram(dists, xlim=(0,0.08)).clone2d()
1733
1734            vpts2 = vpts1.clone().relax_point_positions(n=100, iters=20, sub_iters=10)
1735            dists = vpts2.auto_distance()
1736            h2 = histogram(dists, xlim=(0,0.08)).clone2d()
1737
1738            show([[vpts1, h1], [vpts2, h2]], N=2).close()
1739            ```
1740        """
1741        smooth = vtki.new("PointSmoothingFilter")
1742        smooth.SetInputData(self.dataset)
1743        smooth.SetSmoothingModeToUniform()
1744        smooth.SetNumberOfIterations(iters)
1745        smooth.SetNumberOfSubIterations(sub_iters)
1746        smooth.SetPackingFactor(packing_factor)
1747        if self.point_locator:
1748            smooth.SetLocator(self.point_locator)
1749        if not max_step:
1750            max_step = self.diagonal_size() / 100
1751        smooth.SetMaximumStepSize(max_step)
1752        smooth.SetNeighborhoodSize(n)
1753        if constraints:
1754            fixed_angle = constraints.get("fixed_angle", 45)
1755            boundary_angle = constraints.get("boundary_angle", 110)
1756            smooth.EnableConstraintsOn()
1757            smooth.SetFixedAngle(fixed_angle)
1758            smooth.SetBoundaryAngle(boundary_angle)
1759            smooth.GenerateConstraintScalarsOn()
1760            smooth.GenerateConstraintNormalsOn()
1761        smooth.Update()
1762        self._update(smooth.GetOutput())
1763        self.metadata["PackingRadius"] = smooth.GetPackingRadius()
1764        self.pipeline = utils.OperationNode("relax_point_positions", parents=[self])
1765        return self
1766
1767    def smooth_mls_1d(self, f=0.2, radius=None, n=0) -> Self:
1768        """
1769        Smooth mesh or points with a `Moving Least Squares` variant.
1770        The point data array "Variances" will contain the residue calculated for each point.
1771
1772        Arguments:
1773            f : (float)
1774                smoothing factor - typical range is [0,2].
1775            radius : (float)
1776                radius search in absolute units.
1777                If set then `f` is ignored.
1778            n : (int)
1779                number of neighbours to be used for the fit.
1780                If set then `f` and `radius` are ignored.
1781
1782        Examples:
1783            - [moving_least_squares1D.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/moving_least_squares1D.py)
1784            - [skeletonize.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/skeletonize.py)
1785
1786            ![](https://vedo.embl.es/images/advanced/moving_least_squares1D.png)
1787        """
1788        coords = self.coordinates
1789        ncoords = len(coords)
1790
1791        if n:
1792            Ncp = n
1793        elif radius:
1794            Ncp = 1
1795        else:
1796            Ncp = int(ncoords * f / 10)
1797            if Ncp < 5:
1798                vedo.logger.warning(f"Please choose a fraction higher than {f}")
1799                Ncp = 5
1800
1801        variances, newline = [], []
1802        for p in coords:
1803            points = self.closest_point(p, n=Ncp, radius=radius)
1804            if len(points) < 4:
1805                continue
1806
1807            points = np.array(points)
1808            pointsmean = points.mean(axis=0)  # plane center
1809            _, dd, vv = np.linalg.svd(points - pointsmean)
1810            newp = np.dot(p - pointsmean, vv[0]) * vv[0] + pointsmean
1811            variances.append(dd[1] + dd[2])
1812            newline.append(newp)
1813
1814        self.pointdata["Variances"] = np.array(variances).astype(np.float32)
1815        self.coordinates = newline
1816        self.pipeline = utils.OperationNode("smooth_mls_1d", parents=[self])
1817        return self
1818
1819    def smooth_mls_2d(self, f=0.2, radius=None, n=0) -> Self:
1820        """
1821        Smooth mesh or points with a `Moving Least Squares` algorithm variant.
1822
1823        The `mesh.pointdata['MLSVariance']` array will contain the residue calculated for each point.
1824        When a radius is specified, points that are isolated will not be moved and will get
1825        a 0 entry in array `mesh.pointdata['MLSValidPoint']`.
1826
1827        Arguments:
1828            f : (float)
1829                smoothing factor - typical range is [0, 2].
1830            radius : (float | array)
1831                radius search in absolute units. Can be single value (float) or sequence
1832                for adaptive smoothing. If set then `f` is ignored.
1833            n : (int)
1834                number of neighbours to be used for the fit.
1835                If set then `f` and `radius` are ignored.
1836
1837        Examples:
1838            - [moving_least_squares2D.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/moving_least_squares2D.py)
1839            - [recosurface.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/recosurface.py)
1840
1841                ![](https://vedo.embl.es/images/advanced/recosurface.png)
1842        """
1843        coords = self.coordinates
1844        ncoords = len(coords)
1845
1846        if n:
1847            Ncp = n
1848            radius = None
1849        elif radius is not None:
1850            Ncp = 1
1851        else:
1852            Ncp = int(ncoords * f / 100)
1853            if Ncp < 4:
1854                vedo.logger.error(f"please choose a f-value higher than {f}")
1855                Ncp = 4
1856
1857        variances, newpts, valid = [], [], []
1858        radius_is_sequence = utils.is_sequence(radius)
1859
1860        pb = None
1861        if ncoords > 10000:
1862            pb = utils.ProgressBar(0, ncoords, delay=3)
1863
1864        for i, p in enumerate(coords):
1865            if pb:
1866                pb.print("smooth_mls_2d working ...")
1867            
1868            # if a radius was provided for each point
1869            if radius_is_sequence:
1870                pts = self.closest_point(p, n=Ncp, radius=radius[i])
1871            else:
1872                pts = self.closest_point(p, n=Ncp, radius=radius)
1873
1874            if len(pts) > 3:
1875                ptsmean = pts.mean(axis=0)  # plane center
1876                _, dd, vv = np.linalg.svd(pts - ptsmean)
1877                cv = np.cross(vv[0], vv[1])
1878                t = (np.dot(cv, ptsmean) - np.dot(cv, p)) / np.dot(cv, cv)
1879                newpts.append(p + cv * t)
1880                variances.append(dd[2])
1881                if radius is not None:
1882                    valid.append(1)
1883            else:
1884                newpts.append(p)
1885                variances.append(0)
1886                if radius is not None:
1887                    valid.append(0)
1888
1889        if radius is not None:
1890            self.pointdata["MLSValidPoint"] = np.array(valid).astype(np.uint8)
1891        self.pointdata["MLSVariance"] = np.array(variances).astype(np.float32)
1892
1893        self.coordinates = newpts
1894
1895        self.pipeline = utils.OperationNode("smooth_mls_2d", parents=[self])
1896        return self
1897
1898    def smooth_lloyd_2d(self, iterations=2, bounds=None, options="Qbb Qc Qx") -> Self:
1899        """
1900        Lloyd relaxation of a 2D pointcloud.
1901        
1902        Arguments:
1903            iterations : (int)
1904                number of iterations.
1905            bounds : (list)
1906                bounding box of the domain.
1907            options : (str)
1908                options for the Qhull algorithm.
1909        """
1910        # Credits: https://hatarilabs.com/ih-en/
1911        # tutorial-to-create-a-geospatial-voronoi-sh-mesh-with-python-scipy-and-geopandas
1912        from scipy.spatial import Voronoi as scipy_voronoi
1913
1914        def _constrain_points(points):
1915            # Update any points that have drifted beyond the boundaries of this space
1916            if bounds is not None:
1917                for point in points:
1918                    if point[0] < bounds[0]: point[0] = bounds[0]
1919                    if point[0] > bounds[1]: point[0] = bounds[1]
1920                    if point[1] < bounds[2]: point[1] = bounds[2]
1921                    if point[1] > bounds[3]: point[1] = bounds[3]
1922            return points
1923
1924        def _find_centroid(vertices):
1925            # The equation for the method used here to find the centroid of a
1926            # 2D polygon is given here: https://en.wikipedia.org/wiki/Centroid#Of_a_polygon
1927            area = 0
1928            centroid_x = 0
1929            centroid_y = 0
1930            for i in range(len(vertices) - 1):
1931                step = (vertices[i, 0] * vertices[i + 1, 1]) - (vertices[i + 1, 0] * vertices[i, 1])
1932                centroid_x += (vertices[i, 0] + vertices[i + 1, 0]) * step
1933                centroid_y += (vertices[i, 1] + vertices[i + 1, 1]) * step
1934                area += step
1935            if area:
1936                centroid_x = (1.0 / (3.0 * area)) * centroid_x
1937                centroid_y = (1.0 / (3.0 * area)) * centroid_y
1938            # prevent centroids from escaping bounding box
1939            return _constrain_points([[centroid_x, centroid_y]])[0]
1940
1941        def _relax(voron):
1942            # Moves each point to the centroid of its cell in the voronoi
1943            # map to "relax" the points (i.e. jitter the points so as
1944            # to spread them out within the space).
1945            centroids = []
1946            for idx in voron.point_region:
1947                # the region is a series of indices into voronoi.vertices
1948                # remove point at infinity, designated by index -1
1949                region = [i for i in voron.regions[idx] if i != -1]
1950                # enclose the polygon
1951                region = region + [region[0]]
1952                verts = voron.vertices[region]
1953                # find the centroid of those vertices
1954                centroids.append(_find_centroid(verts))
1955            return _constrain_points(centroids)
1956
1957        if bounds is None:
1958            bounds = self.bounds()
1959
1960        pts = self.vertices[:, (0, 1)]
1961        for i in range(iterations):
1962            vor = scipy_voronoi(pts, qhull_options=options)
1963            _constrain_points(vor.vertices)
1964            pts = _relax(vor)
1965        out = Points(pts)
1966        out.name = "MeshSmoothLloyd2D"
1967        out.pipeline = utils.OperationNode("smooth_lloyd", parents=[self])
1968        return out
1969
1970    def project_on_plane(self, plane="z", point=None, direction=None) -> Self:
1971        """
1972        Project the mesh on one of the Cartesian planes.
1973
1974        Arguments:
1975            plane : (str, Plane)
1976                if plane is `str`, plane can be one of ['x', 'y', 'z'],
1977                represents x-plane, y-plane and z-plane, respectively.
1978                Otherwise, plane should be an instance of `vedo.shapes.Plane`.
1979            point : (float, array)
1980                if plane is `str`, point should be a float represents the intercept.
1981                Otherwise, point is the camera point of perspective projection
1982            direction : (array)
1983                direction of oblique projection
1984
1985        Note:
1986            Parameters `point` and `direction` are only used if the given plane
1987            is an instance of `vedo.shapes.Plane`. And one of these two params
1988            should be left as `None` to specify the projection type.
1989
1990        Example:
1991            ```python
1992            s.project_on_plane(plane='z') # project to z-plane
1993            plane = Plane(pos=(4, 8, -4), normal=(-1, 0, 1), s=(5,5))
1994            s.project_on_plane(plane=plane)                       # orthogonal projection
1995            s.project_on_plane(plane=plane, point=(6, 6, 6))      # perspective projection
1996            s.project_on_plane(plane=plane, direction=(1, 2, -1)) # oblique projection
1997            ```
1998
1999        Examples:
2000            - [silhouette2.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/silhouette2.py)
2001
2002                ![](https://vedo.embl.es/images/basic/silhouette2.png)
2003        """
2004        coords = self.coordinates
2005
2006        if plane == "x":
2007            coords[:, 0] = self.transform.position[0]
2008            intercept = self.xbounds()[0] if point is None else point
2009            self.x(intercept)
2010        elif plane == "y":
2011            coords[:, 1] = self.transform.position[1]
2012            intercept = self.ybounds()[0] if point is None else point
2013            self.y(intercept)
2014        elif plane == "z":
2015            coords[:, 2] = self.transform.position[2]
2016            intercept = self.zbounds()[0] if point is None else point
2017            self.z(intercept)
2018
2019        elif isinstance(plane, vedo.shapes.Plane):
2020            normal = plane.normal / np.linalg.norm(plane.normal)
2021            pl = np.hstack((normal, -np.dot(plane.pos(), normal))).reshape(4, 1)
2022            if direction is None and point is None:
2023                # orthogonal projection
2024                pt = np.hstack((normal, [0])).reshape(4, 1)
2025                # proj_mat = pt.T @ pl * np.eye(4) - pt @ pl.T # python3 only
2026                proj_mat = np.matmul(pt.T, pl) * np.eye(4) - np.matmul(pt, pl.T)
2027
2028            elif direction is None:
2029                # perspective projection
2030                pt = np.hstack((np.array(point), [1])).reshape(4, 1)
2031                # proj_mat = pt.T @ pl * np.eye(4) - pt @ pl.T
2032                proj_mat = np.matmul(pt.T, pl) * np.eye(4) - np.matmul(pt, pl.T)
2033
2034            elif point is None:
2035                # oblique projection
2036                pt = np.hstack((np.array(direction), [0])).reshape(4, 1)
2037                # proj_mat = pt.T @ pl * np.eye(4) - pt @ pl.T
2038                proj_mat = np.matmul(pt.T, pl) * np.eye(4) - np.matmul(pt, pl.T)
2039
2040            coords = np.concatenate([coords, np.ones((coords.shape[:-1] + (1,)))], axis=-1)
2041            # coords = coords @ proj_mat.T
2042            coords = np.matmul(coords, proj_mat.T)
2043            coords = coords[:, :3] / coords[:, 3:]
2044
2045        else:
2046            vedo.logger.error(f"unknown plane {plane}")
2047            raise RuntimeError()
2048
2049        self.alpha(0.1)
2050        self.coordinates = coords
2051        return self
2052
2053    def warp(self, source, target, sigma=1.0, mode="3d") -> Self:
2054        """
2055        "Thin Plate Spline" transformations describe a nonlinear warp transform defined by a set
2056        of source and target landmarks. Any point on the mesh close to a source landmark will
2057        be moved to a place close to the corresponding target landmark.
2058        The points in between are interpolated smoothly using
2059        Bookstein's Thin Plate Spline algorithm.
2060
2061        Transformation object can be accessed with `mesh.transform`.
2062
2063        Arguments:
2064            sigma : (float)
2065                specify the 'stiffness' of the spline.
2066            mode : (str)
2067                set the basis function to either abs(R) (for 3d) or R2LogR (for 2d meshes)
2068
2069        Examples:
2070            - [interpolate_field.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/interpolate_field.py)
2071            - [warp1.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/warp1.py)
2072            - [warp2.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/warp2.py)
2073            - [warp3.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/warp3.py)
2074            - [warp4a.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/warp4a.py)
2075            - [warp4b.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/warp4b.py)
2076            - [warp6.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/warp6.py)
2077
2078            ![](https://vedo.embl.es/images/advanced/warp2.png)
2079        """
2080        parents = [self]
2081
2082        try:
2083            source = source.coordinates
2084            parents.append(source)
2085        except AttributeError:
2086            source = utils.make3d(source)
2087        
2088        try:
2089            target = target.coordinates
2090            parents.append(target)
2091        except AttributeError:
2092            target = utils.make3d(target)
2093
2094        ns = len(source)
2095        nt = len(target)
2096        if ns != nt:
2097            vedo.logger.error(f"#source {ns} != {nt} #target points")
2098            raise RuntimeError()
2099
2100        NLT = NonLinearTransform()
2101        NLT.source_points = source
2102        NLT.target_points = target
2103        self.apply_transform(NLT)
2104
2105        self.pipeline = utils.OperationNode("warp", parents=parents)
2106        return self
2107
2108    def cut_with_plane(
2109            self,
2110            origin=(0, 0, 0),
2111            normal=(1, 0, 0),
2112            invert=False,
2113            # generate_ids=False,
2114    ) -> Self:
2115        """
2116        Cut the mesh with the plane defined by a point and a normal.
2117
2118        Arguments:
2119            origin : (array)
2120                the cutting plane goes through this point
2121            normal : (array)
2122                normal of the cutting plane
2123            invert : (bool)
2124                select which side of the plane to keep
2125
2126        Example:
2127            ```python
2128            from vedo import Cube
2129            cube = Cube().cut_with_plane(normal=(1,1,1))
2130            cube.back_color('pink').show().close()
2131            ```
2132            ![](https://vedo.embl.es/images/feats/cut_with_plane_cube.png)
2133
2134        Examples:
2135            - [trail.py](https://github.com/marcomusy/vedo/blob/master/examples/simulations/trail.py)
2136
2137                ![](https://vedo.embl.es/images/simulations/trail.gif)
2138
2139        Check out also:
2140            `cut_with_box()`, `cut_with_cylinder()`, `cut_with_sphere()`.
2141        """
2142        s = str(normal)
2143        if "x" in s:
2144            normal = (1, 0, 0)
2145            if "-" in s:
2146                normal = -np.array(normal)
2147        elif "y" in s:
2148            normal = (0, 1, 0)
2149            if "-" in s:
2150                normal = -np.array(normal)
2151        elif "z" in s:
2152            normal = (0, 0, 1)
2153            if "-" in s:
2154                normal = -np.array(normal)
2155        plane = vtki.vtkPlane()
2156        plane.SetOrigin(origin)
2157        plane.SetNormal(normal)
2158
2159        clipper = vtki.new("ClipPolyData")
2160        clipper.SetInputData(self.dataset)
2161        clipper.SetClipFunction(plane)
2162        clipper.GenerateClippedOutputOff()
2163        clipper.SetGenerateClipScalars(0)
2164        clipper.SetInsideOut(invert)
2165        clipper.SetValue(0)
2166        clipper.Update()
2167
2168        # if generate_ids:
2169        #     saved_scalars = None # otherwise the scalars are lost
2170        #     if self.dataset.GetPointData().GetScalars():
2171        #         saved_scalars = self.dataset.GetPointData().GetScalars()
2172        #     varr = clipper.GetOutput().GetPointData().GetScalars()
2173        #     if varr.GetName() is None:
2174        #         varr.SetName("DistanceToCut")
2175        #     arr = utils.vtk2numpy(varr)
2176        #     # array of original ids
2177        #     ids = np.arange(arr.shape[0]).astype(int)
2178        #     ids[arr == 0] = -1
2179        #     ids_arr = utils.numpy2vtk(ids, dtype=int)
2180        #     ids_arr.SetName("OriginalIds")
2181        #     clipper.GetOutput().GetPointData().AddArray(ids_arr)
2182        #     if saved_scalars:
2183        #         clipper.GetOutput().GetPointData().AddArray(saved_scalars)
2184
2185        self._update(clipper.GetOutput())
2186        self.pipeline = utils.OperationNode("cut_with_plane", parents=[self])
2187        return self
2188
2189    def cut_with_planes(self, origins, normals, invert=False) -> Self:
2190        """
2191        Cut the mesh with a convex set of planes defined by points and normals.
2192
2193        Arguments:
2194            origins : (array)
2195                each cutting plane goes through this point
2196            normals : (array)
2197                normal of each of the cutting planes
2198            invert : (bool)
2199                if True, cut outside instead of inside
2200
2201        Check out also:
2202            `cut_with_box()`, `cut_with_cylinder()`, `cut_with_sphere()`
2203        """
2204
2205        vpoints = vtki.vtkPoints()
2206        for p in utils.make3d(origins):
2207            vpoints.InsertNextPoint(p)
2208        normals = utils.make3d(normals)
2209
2210        planes = vtki.vtkPlanes()
2211        planes.SetPoints(vpoints)
2212        planes.SetNormals(utils.numpy2vtk(normals, dtype=float))
2213
2214        clipper = vtki.new("ClipPolyData")
2215        clipper.SetInputData(self.dataset)
2216        clipper.SetInsideOut(invert)
2217        clipper.SetClipFunction(planes)
2218        clipper.GenerateClippedOutputOff()
2219        clipper.GenerateClipScalarsOff()
2220        clipper.SetValue(0)
2221        clipper.Update()
2222
2223        self._update(clipper.GetOutput())
2224
2225        self.pipeline = utils.OperationNode("cut_with_planes", parents=[self])
2226        return self
2227
2228    def cut_with_box(self, bounds, invert=False) -> Self:
2229        """
2230        Cut the current mesh with a box or a set of boxes.
2231        This is much faster than `cut_with_mesh()`.
2232
2233        Input `bounds` can be either:
2234        - a Mesh or Points object
2235        - a list of 6 number representing a bounding box `[xmin,xmax, ymin,ymax, zmin,zmax]`
2236        - a list of bounding boxes like the above: `[[xmin1,...], [xmin2,...], ...]`
2237
2238        Example:
2239            ```python
2240            from vedo import Sphere, Cube, show
2241            mesh = Sphere(r=1, res=50)
2242            box  = Cube(side=1.5).wireframe()
2243            mesh.cut_with_box(box)
2244            show(mesh, box, axes=1).close()
2245            ```
2246            ![](https://vedo.embl.es/images/feats/cut_with_box_cube.png)
2247
2248        Check out also:
2249            `cut_with_line()`, `cut_with_plane()`, `cut_with_cylinder()`
2250        """
2251        if isinstance(bounds, Points):
2252            bounds = bounds.bounds()
2253
2254        box = vtki.new("Box")
2255        if utils.is_sequence(bounds[0]):
2256            for bs in bounds:
2257                box.AddBounds(bs)
2258        else:
2259            box.SetBounds(bounds)
2260
2261        clipper = vtki.new("ClipPolyData")
2262        clipper.SetInputData(self.dataset)
2263        clipper.SetClipFunction(box)
2264        clipper.SetInsideOut(not invert)
2265        clipper.GenerateClippedOutputOff()
2266        clipper.GenerateClipScalarsOff()
2267        clipper.SetValue(0)
2268        clipper.Update()
2269        self._update(clipper.GetOutput())
2270
2271        self.pipeline = utils.OperationNode("cut_with_box", parents=[self])
2272        return self
2273
2274    def cut_with_line(self, points, invert=False, closed=True) -> Self:
2275        """
2276        Cut the current mesh with a line vertically in the z-axis direction like a cookie cutter.
2277        The polyline is defined by a set of points (z-coordinates are ignored).
2278        This is much faster than `cut_with_mesh()`.
2279
2280        Check out also:
2281            `cut_with_box()`, `cut_with_plane()`, `cut_with_sphere()`
2282        """
2283        pplane = vtki.new("PolyPlane")
2284        if isinstance(points, Points):
2285            points = points.coordinates.tolist()
2286
2287        if closed:
2288            if isinstance(points, np.ndarray):
2289                points = points.tolist()
2290            points.append(points[0])
2291
2292        vpoints = vtki.vtkPoints()
2293        for p in points:
2294            if len(p) == 2:
2295                p = [p[0], p[1], 0.0]
2296            vpoints.InsertNextPoint(p)
2297
2298        n = len(points)
2299        polyline = vtki.new("PolyLine")
2300        polyline.Initialize(n, vpoints)
2301        polyline.GetPointIds().SetNumberOfIds(n)
2302        for i in range(n):
2303            polyline.GetPointIds().SetId(i, i)
2304        pplane.SetPolyLine(polyline)
2305
2306        clipper = vtki.new("ClipPolyData")
2307        clipper.SetInputData(self.dataset)
2308        clipper.SetClipFunction(pplane)
2309        clipper.SetInsideOut(invert)
2310        clipper.GenerateClippedOutputOff()
2311        clipper.GenerateClipScalarsOff()
2312        clipper.SetValue(0)
2313        clipper.Update()
2314        self._update(clipper.GetOutput())
2315
2316        self.pipeline = utils.OperationNode("cut_with_line", parents=[self])
2317        return self
2318
2319    def cut_with_cookiecutter(self, lines) -> Self:
2320        """
2321        Cut the current mesh with a single line or a set of lines.
2322
2323        Input `lines` can be either:
2324        - a `Mesh` or `Points` object
2325        - a list of 3D points: `[(x1,y1,z1), (x2,y2,z2), ...]`
2326        - a list of 2D points: `[(x1,y1), (x2,y2), ...]`
2327
2328        Example:
2329            ```python
2330            from vedo import *
2331            grid = Mesh(dataurl + "dolfin_fine.vtk")
2332            grid.compute_quality().cmap("Greens")
2333            pols = merge(
2334                Polygon(nsides=10, r=0.3).pos(0.7, 0.3),
2335                Polygon(nsides=10, r=0.2).pos(0.3, 0.7),
2336            )
2337            lines = pols.boundaries()
2338            cgrid = grid.clone().cut_with_cookiecutter(lines)
2339            grid.alpha(0.1).wireframe()
2340            show(grid, cgrid, lines, axes=8, bg='blackboard').close()
2341            ```
2342            ![](https://vedo.embl.es/images/feats/cookiecutter.png)
2343
2344        Check out also:
2345            `cut_with_line()` and `cut_with_point_loop()`
2346
2347        Note:
2348            In case of a warning message like:
2349                "Mesh and trim loop point data attributes are different"
2350            consider interpolating the mesh point data to the loop points,
2351            Eg. (in the above example):
2352            ```python
2353            lines = pols.boundaries().interpolate_data_from(grid, n=2)
2354            ```
2355
2356        Note:
2357            trying to invert the selection by reversing the loop order
2358            will have no effect in this method, hence it does not have
2359            the `invert` option.
2360        """
2361        if utils.is_sequence(lines):
2362            lines = utils.make3d(lines)
2363            iline = list(range(len(lines))) + [0]
2364            poly = utils.buildPolyData(lines, lines=[iline])
2365        else:
2366            poly = lines.dataset
2367
2368        # if invert: # not working
2369        #     rev = vtki.new("ReverseSense")
2370        #     rev.ReverseCellsOn()
2371        #     rev.SetInputData(poly)
2372        #     rev.Update()
2373        #     poly = rev.GetOutput()
2374
2375        # Build loops from the polyline
2376        build_loops = vtki.new("ContourLoopExtraction")
2377        build_loops.SetGlobalWarningDisplay(0)
2378        build_loops.SetInputData(poly)
2379        build_loops.Update()
2380        boundary_poly = build_loops.GetOutput()
2381
2382        ccut = vtki.new("CookieCutter")
2383        ccut.SetInputData(self.dataset)
2384        ccut.SetLoopsData(boundary_poly)
2385        ccut.SetPointInterpolationToMeshEdges()
2386        # ccut.SetPointInterpolationToLoopEdges()
2387        ccut.PassCellDataOn()
2388        ccut.PassPointDataOn()
2389        ccut.Update()
2390        self._update(ccut.GetOutput())
2391
2392        self.pipeline = utils.OperationNode("cut_with_cookiecutter", parents=[self])
2393        return self
2394
2395    def cut_with_cylinder(self, center=(0, 0, 0), axis=(0, 0, 1), r=1, invert=False) -> Self:
2396        """
2397        Cut the current mesh with an infinite cylinder.
2398        This is much faster than `cut_with_mesh()`.
2399
2400        Arguments:
2401            center : (array)
2402                the center of the cylinder
2403            normal : (array)
2404                direction of the cylinder axis
2405            r : (float)
2406                radius of the cylinder
2407
2408        Example:
2409            ```python
2410            from vedo import Disc, show
2411            disc = Disc(r1=1, r2=1.2)
2412            mesh = disc.extrude(3, res=50).linewidth(1)
2413            mesh.cut_with_cylinder([0,0,2], r=0.4, axis='y', invert=True)
2414            show(mesh, axes=1).close()
2415            ```
2416            ![](https://vedo.embl.es/images/feats/cut_with_cylinder.png)
2417
2418        Examples:
2419            - [optics_main1.py](https://github.com/marcomusy/vedo/blob/master/examples/simulations/optics_main1.py)
2420
2421        Check out also:
2422            `cut_with_box()`, `cut_with_plane()`, `cut_with_sphere()`
2423        """
2424        s = str(axis)
2425        if "x" in s:
2426            axis = (1, 0, 0)
2427        elif "y" in s:
2428            axis = (0, 1, 0)
2429        elif "z" in s:
2430            axis = (0, 0, 1)
2431        cyl = vtki.new("Cylinder")
2432        cyl.SetCenter(center)
2433        cyl.SetAxis(axis[0], axis[1], axis[2])
2434        cyl.SetRadius(r)
2435
2436        clipper = vtki.new("ClipPolyData")
2437        clipper.SetInputData(self.dataset)
2438        clipper.SetClipFunction(cyl)
2439        clipper.SetInsideOut(not invert)
2440        clipper.GenerateClippedOutputOff()
2441        clipper.GenerateClipScalarsOff()
2442        clipper.SetValue(0)
2443        clipper.Update()
2444        self._update(clipper.GetOutput())
2445
2446        self.pipeline = utils.OperationNode("cut_with_cylinder", parents=[self])
2447        return self
2448
2449    def cut_with_sphere(self, center=(0, 0, 0), r=1.0, invert=False) -> Self:
2450        """
2451        Cut the current mesh with an sphere.
2452        This is much faster than `cut_with_mesh()`.
2453
2454        Arguments:
2455            center : (array)
2456                the center of the sphere
2457            r : (float)
2458                radius of the sphere
2459
2460        Example:
2461            ```python
2462            from vedo import Disc, show
2463            disc = Disc(r1=1, r2=1.2)
2464            mesh = disc.extrude(3, res=50).linewidth(1)
2465            mesh.cut_with_sphere([1,-0.7,2], r=1.5, invert=True)
2466            show(mesh, axes=1).close()
2467            ```
2468            ![](https://vedo.embl.es/images/feats/cut_with_sphere.png)
2469
2470        Check out also:
2471            `cut_with_box()`, `cut_with_plane()`, `cut_with_cylinder()`
2472        """
2473        sph = vtki.new("Sphere")
2474        sph.SetCenter(center)
2475        sph.SetRadius(r)
2476
2477        clipper = vtki.new("ClipPolyData")
2478        clipper.SetInputData(self.dataset)
2479        clipper.SetClipFunction(sph)
2480        clipper.SetInsideOut(not invert)
2481        clipper.GenerateClippedOutputOff()
2482        clipper.GenerateClipScalarsOff()
2483        clipper.SetValue(0)
2484        clipper.Update()
2485        self._update(clipper.GetOutput())
2486        self.pipeline = utils.OperationNode("cut_with_sphere", parents=[self])
2487        return self
2488
2489    def cut_with_mesh(self, mesh, invert=False, keep=False) -> Union[Self, "vedo.Assembly"]:
2490        """
2491        Cut an `Mesh` mesh with another `Mesh`.
2492
2493        Use `invert` to invert the selection.
2494
2495        Use `keep` to keep the cutoff part, in this case an `Assembly` is returned:
2496        the "cut" object and the "discarded" part of the original object.
2497        You can access both via `assembly.unpack()` method.
2498
2499        Example:
2500        ```python
2501        from vedo import *
2502        arr = np.random.randn(100000, 3)/2
2503        pts = Points(arr).c('red3').pos(5,0,0)
2504        cube = Cube().pos(4,0.5,0)
2505        assem = pts.cut_with_mesh(cube, keep=True)
2506        show(assem.unpack(), axes=1).close()
2507        ```
2508        ![](https://vedo.embl.es/images/feats/cut_with_mesh.png)
2509
2510       Check out also:
2511            `cut_with_box()`, `cut_with_plane()`, `cut_with_cylinder()`
2512       """
2513        polymesh = mesh.dataset
2514        poly = self.dataset
2515
2516        # Create an array to hold distance information
2517        signed_distances = vtki.vtkFloatArray()
2518        signed_distances.SetNumberOfComponents(1)
2519        signed_distances.SetName("SignedDistances")
2520
2521        # implicit function that will be used to slice the mesh
2522        ippd = vtki.new("ImplicitPolyDataDistance")
2523        ippd.SetInput(polymesh)
2524
2525        # Evaluate the signed distance function at all of the grid points
2526        for pointId in range(poly.GetNumberOfPoints()):
2527            p = poly.GetPoint(pointId)
2528            signed_distance = ippd.EvaluateFunction(p)
2529            signed_distances.InsertNextValue(signed_distance)
2530
2531        currentscals = poly.GetPointData().GetScalars()
2532        if currentscals:
2533            currentscals = currentscals.GetName()
2534
2535        poly.GetPointData().AddArray(signed_distances)
2536        poly.GetPointData().SetActiveScalars("SignedDistances")
2537
2538        clipper = vtki.new("ClipPolyData")
2539        clipper.SetInputData(poly)
2540        clipper.SetInsideOut(not invert)
2541        clipper.SetGenerateClippedOutput(keep)
2542        clipper.SetValue(0.0)
2543        clipper.Update()
2544        cpoly = clipper.GetOutput()
2545
2546        if keep:
2547            kpoly = clipper.GetOutput(1)
2548
2549        vis = False
2550        if currentscals:
2551            cpoly.GetPointData().SetActiveScalars(currentscals)
2552            vis = self.mapper.GetScalarVisibility()
2553
2554        self._update(cpoly)
2555
2556        self.pointdata.remove("SignedDistances")
2557        self.mapper.SetScalarVisibility(vis)
2558        if keep:
2559            if isinstance(self, vedo.Mesh):
2560                cutoff = vedo.Mesh(kpoly)
2561            else:
2562                cutoff = vedo.Points(kpoly)
2563            # cutoff = self.__class__(kpoly) # this does not work properly
2564            cutoff.properties = vtki.vtkProperty()
2565            cutoff.properties.DeepCopy(self.properties)
2566            cutoff.actor.SetProperty(cutoff.properties)
2567            cutoff.c("k5").alpha(0.2)
2568            return vedo.Assembly([self, cutoff])
2569
2570        self.pipeline = utils.OperationNode("cut_with_mesh", parents=[self, mesh])
2571        return self
2572
2573    def cut_with_point_loop(
2574        self, points, invert=False, on="points", include_boundary=False
2575    ) -> Self:
2576        """
2577        Cut an `Mesh` object with a set of points forming a closed loop.
2578
2579        Arguments:
2580            invert : (bool)
2581                invert selection (inside-out)
2582            on : (str)
2583                if 'cells' will extract the whole cells lying inside (or outside) the point loop
2584            include_boundary : (bool)
2585                include cells lying exactly on the boundary line. Only relevant on 'cells' mode
2586
2587        Examples:
2588            - [cut_with_points1.py](https://github.com/marcomusy/vedo/blob/master/examples/advanced/cut_with_points1.py)
2589
2590                ![](https://vedo.embl.es/images/advanced/cutWithPoints1.png)
2591
2592            - [cut_with_points2.py](https://github.com/marcomusy/vedo/blob/master/examples/advanced/cut_with_points2.py)
2593
2594                ![](https://vedo.embl.es/images/advanced/cutWithPoints2.png)
2595        """
2596        if isinstance(points, Points):
2597            parents = [points]
2598            vpts = points.dataset.GetPoints()
2599            points = points.coordinates
2600        else:
2601            parents = [self]
2602            vpts = vtki.vtkPoints()
2603            points = utils.make3d(points)
2604            for p in points:
2605                vpts.InsertNextPoint(p)
2606
2607        if "cell" in on:
2608            ippd = vtki.new("ImplicitSelectionLoop")
2609            ippd.SetLoop(vpts)
2610            ippd.AutomaticNormalGenerationOn()
2611            clipper = vtki.new("ExtractPolyDataGeometry")
2612            clipper.SetInputData(self.dataset)
2613            clipper.SetImplicitFunction(ippd)
2614            clipper.SetExtractInside(not invert)
2615            clipper.SetExtractBoundaryCells(include_boundary)
2616        else:
2617            spol = vtki.new("SelectPolyData")
2618            spol.SetLoop(vpts)
2619            spol.GenerateSelectionScalarsOn()
2620            spol.GenerateUnselectedOutputOff()
2621            spol.SetInputData(self.dataset)
2622            spol.Update()
2623            clipper = vtki.new("ClipPolyData")
2624            clipper.SetInputData(spol.GetOutput())
2625            clipper.SetInsideOut(not invert)
2626            clipper.SetValue(0.0)
2627        clipper.Update()
2628        self._update(clipper.GetOutput())
2629
2630        self.pipeline = utils.OperationNode("cut_with_pointloop", parents=parents)
2631        return self
2632
2633    def cut_with_scalar(self, value: float, name="", invert=False) -> Self:
2634        """
2635        Cut a mesh or point cloud with some input scalar point-data.
2636
2637        Arguments:
2638            value : (float)
2639                cutting value
2640            name : (str)
2641                array name of the scalars to be used
2642            invert : (bool)
2643                flip selection
2644
2645        Example:
2646            ```python
2647            from vedo import *
2648            s = Sphere().lw(1)
2649            pts = s.points
2650            scalars = np.sin(3*pts[:,2]) + pts[:,0]
2651            s.pointdata["somevalues"] = scalars
2652            s.cut_with_scalar(0.3)
2653            s.cmap("Spectral", "somevalues").add_scalarbar()
2654            s.show(axes=1).close()
2655            ```
2656            ![](https://vedo.embl.es/images/feats/cut_with_scalars.png)
2657        """
2658        if name:
2659            self.pointdata.select(name)
2660        clipper = vtki.new("ClipPolyData")
2661        clipper.SetInputData(self.dataset)
2662        clipper.SetValue(value)
2663        clipper.GenerateClippedOutputOff()
2664        clipper.SetInsideOut(not invert)
2665        clipper.Update()
2666        self._update(clipper.GetOutput())
2667        self.pipeline = utils.OperationNode("cut_with_scalar", parents=[self])
2668        return self
2669
2670    def crop(self,
2671             top=None, bottom=None, right=None, left=None, front=None, back=None,
2672             bounds=()) -> Self:
2673        """
2674        Crop an `Mesh` object.
2675
2676        Arguments:
2677            top : (float)
2678                fraction to crop from the top plane (positive z)
2679            bottom : (float)
2680                fraction to crop from the bottom plane (negative z)
2681            front : (float)
2682                fraction to crop from the front plane (positive y)
2683            back : (float)
2684                fraction to crop from the back plane (negative y)
2685            right : (float)
2686                fraction to crop from the right plane (positive x)
2687            left : (float)
2688                fraction to crop from the left plane (negative x)
2689            bounds : (list)
2690                bounding box of the crop region as `[x0,x1, y0,y1, z0,z1]`
2691
2692        Example:
2693            ```python
2694            from vedo import Sphere
2695            Sphere().crop(right=0.3, left=0.1).show()
2696            ```
2697            ![](https://user-images.githubusercontent.com/32848391/57081955-0ef1e800-6cf6-11e9-99de-b45220939bc9.png)
2698        """
2699        if not len(bounds):
2700            pos = np.array(self.pos())
2701            x0, x1, y0, y1, z0, z1 = self.bounds()
2702            x0, y0, z0 = [x0, y0, z0] - pos
2703            x1, y1, z1 = [x1, y1, z1] - pos
2704
2705            dx, dy, dz = x1 - x0, y1 - y0, z1 - z0
2706            if top:
2707                z1 = z1 - top * dz
2708            if bottom:
2709                z0 = z0 + bottom * dz
2710            if front:
2711                y1 = y1 - front * dy
2712            if back:
2713                y0 = y0 + back * dy
2714            if right:
2715                x1 = x1 - right * dx
2716            if left:
2717                x0 = x0 + left * dx
2718            bounds = (x0, x1, y0, y1, z0, z1)
2719
2720        cu = vtki.new("Box")
2721        cu.SetBounds(bounds)
2722
2723        clipper = vtki.new("ClipPolyData")
2724        clipper.SetInputData(self.dataset)
2725        clipper.SetClipFunction(cu)
2726        clipper.InsideOutOn()
2727        clipper.GenerateClippedOutputOff()
2728        clipper.GenerateClipScalarsOff()
2729        clipper.SetValue(0)
2730        clipper.Update()
2731        self._update(clipper.GetOutput())
2732
2733        self.pipeline = utils.OperationNode(
2734            "crop", parents=[self], comment=f"#pts {self.dataset.GetNumberOfPoints()}"
2735        )
2736        return self
2737
2738    def generate_surface_halo(
2739            self, 
2740            distance=0.05,
2741            res=(50, 50, 50),
2742            bounds=(),
2743            maxdist=None,
2744    ) -> "vedo.Mesh":
2745        """
2746        Generate the surface halo which sits at the specified distance from the input one.
2747
2748        Arguments:
2749            distance : (float)
2750                distance from the input surface
2751            res : (int)
2752                resolution of the surface
2753            bounds : (list)
2754                bounding box of the surface
2755            maxdist : (float)
2756                maximum distance to generate the surface
2757        """
2758        if not bounds:
2759            bounds = self.bounds()
2760
2761        if not maxdist:
2762            maxdist = self.diagonal_size() / 2
2763
2764        imp = vtki.new("ImplicitModeller")
2765        imp.SetInputData(self.dataset)
2766        imp.SetSampleDimensions(res)
2767        if maxdist:
2768            imp.SetMaximumDistance(maxdist)
2769        if len(bounds) == 6:
2770            imp.SetModelBounds(bounds)
2771        contour = vtki.new("ContourFilter")
2772        contour.SetInputConnection(imp.GetOutputPort())
2773        contour.SetValue(0, distance)
2774        contour.Update()
2775        out = vedo.Mesh(contour.GetOutput())
2776        out.c("lightblue").alpha(0.25).lighting("off")
2777        out.pipeline = utils.OperationNode("generate_surface_halo", parents=[self])
2778        return out
2779
2780    def generate_mesh(
2781        self,
2782        line_resolution=None,
2783        mesh_resolution=None,
2784        smooth=0.0,
2785        jitter=0.001,
2786        grid=None,
2787        quads=False,
2788        invert=False,
2789    ) -> Self:
2790        """
2791        Generate a polygonal Mesh from a closed contour line.
2792        If line is not closed it will be closed with a straight segment.
2793
2794        Check also `generate_delaunay2d()`.
2795
2796        Arguments:
2797            line_resolution : (int)
2798                resolution of the contour line. The default is None, in this case
2799                the contour is not resampled.
2800            mesh_resolution : (int)
2801                resolution of the internal triangles not touching the boundary.
2802            smooth : (float)
2803                smoothing of the contour before meshing.
2804            jitter : (float)
2805                add a small noise to the internal points.
2806            grid : (Grid)
2807                manually pass a Grid object. The default is True.
2808            quads : (bool)
2809                generate a mesh of quads instead of triangles.
2810            invert : (bool)
2811                flip the line orientation. The default is False.
2812
2813        Examples:
2814            - [line2mesh_tri.py](https://github.com/marcomusy/vedo/blob/master/examples/advanced/line2mesh_tri.py)
2815
2816                ![](https://vedo.embl.es/images/advanced/line2mesh_tri.jpg)
2817
2818            - [line2mesh_quads.py](https://github.com/marcomusy/vedo/blob/master/examples/advanced/line2mesh_quads.py)
2819
2820                ![](https://vedo.embl.es/images/advanced/line2mesh_quads.png)
2821        """
2822        if line_resolution is None:
2823            contour = vedo.shapes.Line(self.coordinates)
2824        else:
2825            contour = vedo.shapes.Spline(self.coordinates, smooth=smooth, res=line_resolution)
2826        contour.clean()
2827
2828        length = contour.length()
2829        density = length / contour.npoints
2830        # print(f"tomesh():\n\tline length = {length}")
2831        # print(f"\tdensity = {density} length/pt_separation")
2832
2833        x0, x1 = contour.xbounds()
2834        y0, y1 = contour.ybounds()
2835
2836        if grid is None:
2837            if mesh_resolution is None:
2838                resx = int((x1 - x0) / density + 0.5)
2839                resy = int((y1 - y0) / density + 0.5)
2840                # print(f"tmesh_resolution = {[resx, resy]}")
2841            else:
2842                if utils.is_sequence(mesh_resolution):
2843                    resx, resy = mesh_resolution
2844                else:
2845                    resx, resy = mesh_resolution, mesh_resolution
2846            grid = vedo.shapes.Grid(
2847                [(x0 + x1) / 2, (y0 + y1) / 2, 0],
2848                s=((x1 - x0) * 1.025, (y1 - y0) * 1.025),
2849                res=(resx, resy),
2850            )
2851        else:
2852            grid = grid.clone()
2853
2854        cpts = contour.coordinates
2855
2856        # make sure it's closed
2857        p0, p1 = cpts[0], cpts[-1]
2858        nj = max(2, int(utils.mag(p1 - p0) / density + 0.5))
2859        joinline = vedo.shapes.Line(p1, p0, res=nj)
2860        contour = vedo.merge(contour, joinline).subsample(0.0001)
2861
2862        ####################################### quads
2863        if quads:
2864            cmesh = grid.clone().cut_with_point_loop(contour, on="cells", invert=invert)
2865            cmesh.wireframe(False).lw(0.5)
2866            cmesh.pipeline = utils.OperationNode(
2867                "generate_mesh",
2868                parents=[self, contour],
2869                comment=f"#quads {cmesh.dataset.GetNumberOfCells()}",
2870            )
2871            return cmesh
2872        #############################################
2873
2874        grid_tmp = grid.coordinates.copy()
2875
2876        if jitter:
2877            np.random.seed(0)
2878            sigma = 1.0 / np.sqrt(grid.npoints) * grid.diagonal_size() * jitter
2879            # print(f"\tsigma jittering = {sigma}")
2880            grid_tmp += np.random.rand(grid.npoints, 3) * sigma
2881            grid_tmp[:, 2] = 0.0
2882
2883        todel = []
2884        density /= np.sqrt(3)
2885        vgrid_tmp = Points(grid_tmp)
2886
2887        for p in contour.coordinates:
2888            out = vgrid_tmp.closest_point(p, radius=density, return_point_id=True)
2889            todel += out.tolist()
2890
2891        grid_tmp = grid_tmp.tolist()
2892        for index in sorted(list(set(todel)), reverse=True):
2893            del grid_tmp[index]
2894
2895        points = contour.coordinates.tolist() + grid_tmp
2896        if invert:
2897            boundary = list(reversed(range(contour.npoints)))
2898        else:
2899            boundary = list(range(contour.npoints))
2900
2901        dln = Points(points).generate_delaunay2d(mode="xy", boundaries=[boundary])
2902        dln.compute_normals(points=False)  # fixes reversd faces
2903        dln.lw(1)
2904
2905        dln.pipeline = utils.OperationNode(
2906            "generate_mesh",
2907            parents=[self, contour],
2908            comment=f"#cells {dln.dataset.GetNumberOfCells()}",
2909        )
2910        return dln
2911
2912    def reconstruct_surface(
2913        self,
2914        dims=(100, 100, 100),
2915        radius=None,
2916        sample_size=None,
2917        hole_filling=True,
2918        bounds=(),
2919        padding=0.05,
2920    ) -> "vedo.Mesh":
2921        """
2922        Surface reconstruction from a scattered cloud of points.
2923
2924        Arguments:
2925            dims : (int)
2926                number of voxels in x, y and z to control precision.
2927            radius : (float)
2928                radius of influence of each point.
2929                Smaller values generally improve performance markedly.
2930                Note that after the signed distance function is computed,
2931                any voxel taking on the value >= radius
2932                is presumed to be "unseen" or uninitialized.
2933            sample_size : (int)
2934                if normals are not present
2935                they will be calculated using this sample size per point.
2936            hole_filling : (bool)
2937                enables hole filling, this generates
2938                separating surfaces between the empty and unseen portions of the volume.
2939            bounds : (list)
2940                region in space in which to perform the sampling
2941                in format (xmin,xmax, ymin,ymax, zim, zmax)
2942            padding : (float)
2943                increase by this fraction the bounding box
2944
2945        Examples:
2946            - [recosurface.py](https://github.com/marcomusy/vedo/blob/master/examples/advanced/recosurface.py)
2947
2948                ![](https://vedo.embl.es/images/advanced/recosurface.png)
2949        """
2950        if not utils.is_sequence(dims):
2951            dims = (dims, dims, dims)
2952
2953        sdf = vtki.new("SignedDistance")
2954
2955        if len(bounds) == 6:
2956            sdf.SetBounds(bounds)
2957        else:
2958            x0, x1, y0, y1, z0, z1 = self.bounds()
2959            sdf.SetBounds(
2960                x0 - (x1 - x0) * padding,
2961                x1 + (x1 - x0) * padding,
2962                y0 - (y1 - y0) * padding,
2963                y1 + (y1 - y0) * padding,
2964                z0 - (z1 - z0) * padding,
2965                z1 + (z1 - z0) * padding,
2966            )
2967        
2968        bb = sdf.GetBounds()
2969        if bb[0]==bb[1]:
2970            vedo.logger.warning("reconstruct_surface(): zero x-range")
2971        if bb[2]==bb[3]:
2972            vedo.logger.warning("reconstruct_surface(): zero y-range")
2973        if bb[4]==bb[5]:
2974            vedo.logger.warning("reconstruct_surface(): zero z-range")
2975
2976        pd = self.dataset
2977
2978        if pd.GetPointData().GetNormals():
2979            sdf.SetInputData(pd)
2980        else:
2981            normals = vtki.new("PCANormalEstimation")
2982            normals.SetInputData(pd)
2983            if not sample_size:
2984                sample_size = int(pd.GetNumberOfPoints() / 50)
2985            normals.SetSampleSize(sample_size)
2986            normals.SetNormalOrientationToGraphTraversal()
2987            sdf.SetInputConnection(normals.GetOutputPort())
2988            # print("Recalculating normals with sample size =", sample_size)
2989
2990        if radius is None:
2991            radius = self.diagonal_size() / (sum(dims) / 3) * 5
2992            # print("Calculating mesh from points with radius =", radius)
2993
2994        sdf.SetRadius(radius)
2995        sdf.SetDimensions(dims)
2996        sdf.Update()
2997
2998        surface = vtki.new("ExtractSurface")
2999        surface.SetRadius(radius * 0.99)
3000        surface.SetHoleFilling(hole_filling)
3001        surface.ComputeNormalsOff()
3002        surface.ComputeGradientsOff()
3003        surface.SetInputConnection(sdf.GetOutputPort())
3004        surface.Update()
3005        m = vedo.mesh.Mesh(surface.GetOutput(), c=self.color())
3006
3007        m.pipeline = utils.OperationNode(
3008            "reconstruct_surface",
3009            parents=[self],
3010            comment=f"#pts {m.dataset.GetNumberOfPoints()}",
3011        )
3012        return m
3013
3014    def compute_clustering(self, radius: float) -> Self:
3015        """
3016        Cluster points in space. The `radius` is the radius of local search.
3017        
3018        An array named "ClusterId" is added to `pointdata`.
3019
3020        Examples:
3021            - [clustering.py](https://github.com/marcomusy/vedo/blob/master/examples/basic/clustering.py)
3022
3023                ![](https://vedo.embl.es/images/basic/clustering.png)
3024        """
3025        cluster = vtki.new("EuclideanClusterExtraction")
3026        cluster.SetInputData(self.dataset)
3027        cluster.SetExtractionModeToAllClusters()
3028        cluster.SetRadius(radius)
3029        cluster.ColorClustersOn()
3030        cluster.Update()
3031        idsarr = cluster.GetOutput().GetPointData().GetArray("ClusterId")
3032        self.dataset.GetPointData().AddArray(idsarr)
3033        self.pipeline = utils.OperationNode(
3034            "compute_clustering", parents=[self], comment=f"radius = {radius}"
3035        )
3036        return self
3037
3038    def compute_connections(self, radius, mode=0, regions=(), vrange=(0, 1), seeds=(), angle=0.0) -> Self:
3039        """
3040        Extracts and/or segments points from a point cloud based on geometric distance measures
3041        (e.g., proximity, normal alignments, etc.) and optional measures such as scalar range.
3042        The default operation is to segment the points into "connected" regions where the connection
3043        is determined by an appropriate distance measure. Each region is given a region id.
3044
3045        Optionally, the filter can output the largest connected region of points; a particular region
3046        (via id specification); those regions that are seeded using a list of input point ids;
3047        or the region of points closest to a specified position.
3048
3049        The key parameter of this filter is the radius defining a sphere around each point which defines
3050        a local neighborhood: any other points in the local neighborhood are assumed connected to the point.
3051        Note that the radius is defined in absolute terms.
3052
3053        Other parameters are used to further qualify what it means to be a neighboring point.
3054        For example, scalar range and/or point normals can be used to further constrain the neighborhood.
3055        Also the extraction mode defines how the filter operates.
3056        By default, all regions are extracted but it is possible to extract particular regions;
3057        the region closest to a seed point; seeded regions; or the largest region found while processing.
3058        By default, all regions are extracted.
3059
3060        On output, all points are labeled with a region number.
3061        However note that the number of input and output points may not be the same:
3062        if not extracting all regions then the output size may be less than the input size.
3063
3064        Arguments:
3065            radius : (float)
3066                variable specifying a local sphere used to define local point neighborhood
3067            mode : (int)
3068                - 0,  Extract all regions
3069                - 1,  Extract point seeded regions
3070                - 2,  Extract largest region
3071                - 3,  Test specified regions
3072                - 4,  Extract all regions with scalar connectivity
3073                - 5,  Extract point seeded regions
3074            regions : (list)
3075                a list of non-negative regions id to extract
3076            vrange : (list)
3077                scalar range to use to extract points based on scalar connectivity
3078            seeds : (list)
3079                a list of non-negative point seed ids
3080            angle : (list)
3081                points are connected if the angle between their normals is
3082                within this angle threshold (expressed in degrees).
3083        """
3084        # https://vtk.org/doc/nightly/html/classvtkConnectedPointsFilter.html
3085        cpf = vtki.new("ConnectedPointsFilter")
3086        cpf.SetInputData(self.dataset)
3087        cpf.SetRadius(radius)
3088        if mode == 0:  # Extract all regions
3089            pass
3090
3091        elif mode == 1:  # Extract point seeded regions
3092            cpf.SetExtractionModeToPointSeededRegions()
3093            for s in seeds:
3094                cpf.AddSeed(s)
3095
3096        elif mode == 2:  # Test largest region
3097            cpf.SetExtractionModeToLargestRegion()
3098
3099        elif mode == 3:  # Test specified regions
3100            cpf.SetExtractionModeToSpecifiedRegions()
3101            for r in regions:
3102                cpf.AddSpecifiedRegion(r)
3103
3104        elif mode == 4:  # Extract all regions with scalar connectivity
3105            cpf.SetExtractionModeToLargestRegion()
3106            cpf.ScalarConnectivityOn()
3107            cpf.SetScalarRange(vrange[0], vrange[1])
3108
3109        elif mode == 5:  # Extract point seeded regions
3110            cpf.SetExtractionModeToLargestRegion()
3111            cpf.ScalarConnectivityOn()
3112            cpf.SetScalarRange(vrange[0], vrange[1])
3113            cpf.AlignedNormalsOn()
3114            cpf.SetNormalAngle(angle)
3115
3116        cpf.Update()
3117        self._update(cpf.GetOutput(), reset_locators=False)
3118        return self
3119
3120    def compute_camera_distance(self) -> np.ndarray:
3121        """
3122        Calculate the distance from points to the camera.
3123        
3124        A pointdata array is created with name 'DistanceToCamera' and returned.
3125        """
3126        if vedo.plotter_instance and vedo.plotter_instance.renderer:
3127            poly = self.dataset
3128            dc = vtki.new("DistanceToCamera")
3129            dc.SetInputData(poly)
3130            dc.SetRenderer(vedo.plotter_instance.renderer)
3131            dc.Update()
3132            self._update(dc.GetOutput(), reset_locators=False)
3133            return self.pointdata["DistanceToCamera"]
3134        return np.array([])
3135
3136    def densify(self, target_distance=0.1, nclosest=6, radius=None, niter=1, nmax=None) -> Self:
3137        """
3138        Return a copy of the cloud with new added points.
3139        The new points are created in such a way that all points in any local neighborhood are
3140        within a target distance of one another.
3141
3142        For each input point, the distance to all points in its neighborhood is computed.
3143        If any of its neighbors is further than the target distance,
3144        the edge connecting the point and its neighbor is bisected and
3145        a new point is inserted at the bisection point.
3146        A single pass is completed once all the input points are visited.
3147        Then the process repeats to the number of iterations.
3148
3149        Examples:
3150            - [densifycloud.py](https://github.com/marcomusy/vedo/blob/master/examples/volumetric/densifycloud.py)
3151
3152                ![](https://vedo.embl.es/images/volumetric/densifycloud.png)
3153
3154        .. note::
3155            Points will be created in an iterative fashion until all points in their
3156            local neighborhood are the target distance apart or less.
3157            Note that the process may terminate early due to the
3158            number of iterations. By default the target distance is set to 0.5.
3159            Note that the target_distance should be less than the radius
3160            or nothing will change on output.
3161
3162        .. warning::
3163            This class can generate a lot of points very quickly.
3164            The maximum number of iterations is by default set to =1.0 for this reason.
3165            Increase the number of iterations very carefully.
3166            Also, `nmax` can be set to limit the explosion of points.
3167            It is also recommended that a N closest neighborhood is used.
3168
3169        """
3170        src = vtki.new("ProgrammableSource")
3171        opts = self.coordinates
3172        # zeros = np.zeros(3)
3173
3174        def _read_points():
3175            output = src.GetPolyDataOutput()
3176            points = vtki.vtkPoints()
3177            for p in opts:
3178                # print(p)
3179                # if not np.array_equal(p, zeros):
3180                points.InsertNextPoint(p)
3181            output.SetPoints(points)
3182
3183        src.SetExecuteMethod(_read_points)
3184
3185        dens = vtki.new("DensifyPointCloudFilter")
3186        dens.SetInputConnection(src.GetOutputPort())
3187        # dens.SetInputData(self.dataset) # this does not work
3188        dens.InterpolateAttributeDataOn()
3189        dens.SetTargetDistance(target_distance)
3190        dens.SetMaximumNumberOfIterations(niter)
3191        if nmax:
3192            dens.SetMaximumNumberOfPoints(nmax)
3193
3194        if radius:
3195            dens.SetNeighborhoodTypeToRadius()
3196            dens.SetRadius(radius)
3197        elif nclosest:
3198            dens.SetNeighborhoodTypeToNClosest()
3199            dens.SetNumberOfClosestPoints(nclosest)
3200        else:
3201            vedo.logger.error("set either radius or nclosest")
3202            raise RuntimeError()
3203        dens.Update()
3204
3205        cld = Points(dens.GetOutput())
3206        cld.copy_properties_from(self)
3207        cld.interpolate_data_from(self, n=nclosest, radius=radius)
3208        cld.name = "DensifiedCloud"
3209        cld.pipeline = utils.OperationNode(
3210            "densify",
3211            parents=[self],
3212            c="#e9c46a:",
3213            comment=f"#pts {cld.dataset.GetNumberOfPoints()}",
3214        )
3215        return cld
3216
3217    ###############################################################################
3218    ## stuff returning a Volume
3219    ###############################################################################
3220
3221    def density(
3222        self, dims=(40, 40, 40), bounds=None, radius=None, compute_gradient=False, locator=None
3223    ) -> "vedo.Volume":
3224        """
3225        Generate a density field from a point cloud. Input can also be a set of 3D coordinates.
3226        Output is a `Volume`.
3227
3228        The local neighborhood is specified as the `radius` around each sample position (each voxel).
3229        If left to None, the radius is automatically computed as the diagonal of the bounding box
3230        and can be accessed via `vol.metadata["radius"]`.
3231        The density is expressed as the number of counts in the radius search.
3232
3233        Arguments:
3234            dims : (int, list)
3235                number of voxels in x, y and z of the output Volume.
3236            compute_gradient : (bool)
3237                Turn on/off the generation of the gradient vector,
3238                gradient magnitude scalar, and function classification scalar.
3239                By default this is off. Note that this will increase execution time
3240                and the size of the output. (The names of these point data arrays are:
3241                "Gradient", "Gradient Magnitude", and "Classification")
3242            locator : (vtkPointLocator)
3243                can be assigned from a previous call for speed (access it via `object.point_locator`).
3244
3245        Examples:
3246            - [plot_density3d.py](https://github.com/marcomusy/vedo/blob/master/examples/pyplot/plot_density3d.py)
3247
3248                ![](https://vedo.embl.es/images/pyplot/plot_density3d.png)
3249        """
3250        pdf = vtki.new("PointDensityFilter")
3251        pdf.SetInputData(self.dataset)
3252
3253        if not utils.is_sequence(dims):
3254            dims = [dims, dims, dims]
3255
3256        if bounds is None:
3257            bounds = list(self.bounds())
3258        elif len(bounds) == 4:
3259            bounds = [*bounds, 0, 0]
3260
3261        if bounds[5] - bounds[4] == 0 or len(dims) == 2:  # its 2D
3262            dims = list(dims)
3263            dims = [dims[0], dims[1], 2]
3264            diag = self.diagonal_size()
3265            bounds[5] = bounds[4] + diag / 1000
3266        pdf.SetModelBounds(bounds)
3267
3268        pdf.SetSampleDimensions(dims)
3269
3270        if locator:
3271            pdf.SetLocator(locator)
3272
3273        pdf.SetDensityEstimateToFixedRadius()
3274        if radius is None:
3275            radius = self.diagonal_size() / 20
3276        pdf.SetRadius(radius)
3277        pdf.SetComputeGradient(compute_gradient)
3278        pdf.Update()
3279
3280        vol = vedo.Volume(pdf.GetOutput()).mode(1)
3281        vol.name = "PointDensity"
3282        vol.metadata["radius"] = radius
3283        vol.locator = pdf.GetLocator()
3284        vol.pipeline = utils.OperationNode(
3285            "density", parents=[self], comment=f"dims={tuple(vol.dimensions())}"
3286        )
3287        return vol
3288
3289
3290    def tovolume(
3291        self,
3292        kernel="shepard",
3293        radius=None,
3294        n=None,
3295        bounds=None,
3296        null_value=None,
3297        dims=(25, 25, 25),
3298    ) -> "vedo.Volume":
3299        """
3300        Generate a `Volume` by interpolating a scalar
3301        or vector field which is only known on a scattered set of points or mesh.
3302        Available interpolation kernels are: shepard, gaussian, or linear.
3303
3304        Arguments:
3305            kernel : (str)
3306                interpolation kernel type [shepard]
3307            radius : (float)
3308                radius of the local search
3309            n : (int)
3310                number of point to use for interpolation
3311            bounds : (list)
3312                bounding box of the output Volume object
3313            dims : (list)
3314                dimensions of the output Volume object
3315            null_value : (float)
3316                value to be assigned to invalid points
3317
3318        Examples:
3319            - [interpolate_volume.py](https://github.com/marcomusy/vedo/blob/master/examples/volumetric/interpolate_volume.py)
3320
3321                ![](https://vedo.embl.es/images/volumetric/59095175-1ec5a300-8918-11e9-8bc0-fd35c8981e2b.jpg)
3322        """
3323        if radius is None and not n:
3324            vedo.logger.error("please set either radius or n")
3325            raise RuntimeError
3326
3327        poly = self.dataset
3328
3329        # Create a probe volume
3330        probe = vtki.vtkImageData()
3331        probe.SetDimensions(dims)
3332        if bounds is None:
3333            bounds = self.bounds()
3334        probe.SetOrigin(bounds[0], bounds[2], bounds[4])
3335        probe.SetSpacing(
3336            (bounds[1] - bounds[0]) / dims[0],
3337            (bounds[3] - bounds[2]) / dims[1],
3338            (bounds[5] - bounds[4]) / dims[2],
3339        )
3340
3341        if not self.point_locator:
3342            self.point_locator = vtki.new("PointLocator")
3343            self.point_locator.SetDataSet(poly)
3344            self.point_locator.BuildLocator()
3345
3346        if kernel == "shepard":
3347            kern = vtki.new("ShepardKernel")
3348            kern.SetPowerParameter(2)
3349        elif kernel == "gaussian":
3350            kern = vtki.new("GaussianKernel")
3351        elif kernel == "linear":
3352            kern = vtki.new("LinearKernel")
3353        else:
3354            vedo.logger.error("Error in tovolume(), available kernels are:")
3355            vedo.logger.error(" [shepard, gaussian, linear]")
3356            raise RuntimeError()
3357
3358        if radius:
3359            kern.SetRadius(radius)
3360
3361        interpolator = vtki.new("PointInterpolator")
3362        interpolator.SetInputData(probe)
3363        interpolator.SetSourceData(poly)
3364        interpolator.SetKernel(kern)
3365        interpolator.SetLocator(self.point_locator)
3366
3367        if n:
3368            kern.SetNumberOfPoints(n)
3369            kern.SetKernelFootprintToNClosest()
3370        else:
3371            kern.SetRadius(radius)
3372
3373        if null_value is not None:
3374            interpolator.SetNullValue(null_value)
3375        else:
3376            interpolator.SetNullPointsStrategyToClosestPoint()
3377        interpolator.Update()
3378
3379        vol = vedo.Volume(interpolator.GetOutput())
3380
3381        vol.pipeline = utils.OperationNode(
3382            "signed_distance",
3383            parents=[self],
3384            comment=f"dims={tuple(vol.dimensions())}",
3385            c="#e9c46a:#0096c7",
3386        )
3387        return vol
3388
3389    #################################################################################    
3390    def generate_segments(self, istart=0, rmax=1e30, niter=3) -> "vedo.shapes.Lines":
3391        """
3392        Generate a line segments from a set of points.
3393        The algorithm is based on the closest point search.
3394
3395        Returns a `Line` object.
3396        This object contains the a metadata array of used vertex counts in "UsedVertexCount"
3397        and the sum of the length of the segments in "SegmentsLengthSum".
3398
3399        Arguments:
3400            istart : (int)
3401                index of the starting point
3402            rmax : (float)
3403                maximum length of a segment
3404            niter : (int)
3405                number of iterations or passes through the points
3406
3407        Examples:
3408            - [moving_least_squares1D.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/moving_least_squares1D.py)
3409        """
3410        points = self.coordinates
3411        segments = []
3412        dists = []
3413        n = len(points)
3414        used = np.zeros(n, dtype=int)
3415        for _ in range(niter):
3416            i = istart
3417            for _ in range(n):
3418                p = points[i]
3419                ids = self.closest_point(p, n=4, return_point_id=True)
3420                j = ids[1]
3421                if used[j] > 1 or [j, i] in segments:
3422                    j = ids[2]
3423                if used[j] > 1:
3424                    j = ids[3]
3425                d = np.linalg.norm(p - points[j])
3426                if used[j] > 1 or used[i] > 1 or d > rmax:
3427                    i += 1
3428                    if i >= n:
3429                        i = 0
3430                    continue
3431                used[i] += 1
3432                used[j] += 1
3433                segments.append([i, j])
3434                dists.append(d)
3435                i = j
3436        segments = np.array(segments, dtype=int)
3437
3438        lines = vedo.shapes.Lines(points[segments], c="k", lw=3)
3439        lines.metadata["UsedVertexCount"] = used
3440        lines.metadata["SegmentsLengthSum"] = np.sum(dists)
3441        lines.pipeline = utils.OperationNode("generate_segments", parents=[self])
3442        lines.name = "Segments"
3443        return lines
3444
3445    def generate_delaunay2d(
3446        self,
3447        mode="scipy",
3448        boundaries=(),
3449        tol=None,
3450        alpha=0.0,
3451        offset=0.0,
3452        transform=None,
3453    ) -> "vedo.mesh.Mesh":
3454        """
3455        Create a mesh from points in the XY plane.
3456        If `mode='fit'` then the filter computes a best fitting
3457        plane and projects the points onto it.
3458
3459        Check also `generate_mesh()`.
3460
3461        Arguments:
3462            tol : (float)
3463                specify a tolerance to control discarding of closely spaced points.
3464                This tolerance is specified as a fraction of the diagonal length of the bounding box of the points.
3465            alpha : (float)
3466                for a non-zero alpha value, only edges or triangles contained
3467                within a sphere centered at mesh vertices will be output.
3468                Otherwise, only triangles will be output.
3469            offset : (float)
3470                multiplier to control the size of the initial, bounding Delaunay triangulation.
3471            transform: (LinearTransform, NonLinearTransform)
3472                a transformation which is applied to points to generate a 2D problem.
3473                This maps a 3D dataset into a 2D dataset where triangulation can be done on the XY plane.
3474                The points are transformed and triangulated.
3475                The topology of triangulated points is used as the output topology.
3476
3477        Examples:
3478            - [delaunay2d.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/delaunay2d.py)
3479
3480                ![](https://vedo.embl.es/images/basic/delaunay2d.png)
3481        """
3482        plist = self.coordinates.copy()
3483
3484        #########################################################
3485        if mode == "scipy":
3486            from scipy.spatial import Delaunay as scipy_delaunay
3487
3488            tri = scipy_delaunay(plist[:, 0:2])
3489            return vedo.mesh.Mesh([plist, tri.simplices])
3490        ##########################################################
3491
3492        pd = vtki.vtkPolyData()
3493        vpts = vtki.vtkPoints()
3494        vpts.SetData(utils.numpy2vtk(plist, dtype=np.float32))
3495        pd.SetPoints(vpts)
3496
3497        delny = vtki.new("Delaunay2D")
3498        delny.SetInputData(pd)
3499        if tol:
3500            delny.SetTolerance(tol)
3501        delny.SetAlpha(alpha)
3502        delny.SetOffset(offset)
3503
3504        if transform:
3505            delny.SetTransform(transform.T)
3506        elif mode == "fit":
3507            delny.SetProjectionPlaneMode(vtki.get_class("VTK_BEST_FITTING_PLANE"))
3508        elif mode == "xy" and boundaries:
3509            boundary = vtki.vtkPolyData()
3510            boundary.SetPoints(vpts)
3511            cell_array = vtki.vtkCellArray()
3512            for b in boundaries:
3513                cpolygon = vtki.vtkPolygon()
3514                for idd in b:
3515                    cpolygon.GetPointIds().InsertNextId(idd)
3516                cell_array.InsertNextCell(cpolygon)
3517            boundary.SetPolys(cell_array)
3518            delny.SetSourceData(boundary)
3519
3520        delny.Update()
3521
3522        msh = vedo.mesh.Mesh(delny.GetOutput())
3523        msh.name = "Delaunay2D"
3524        msh.clean().lighting("off")
3525        msh.pipeline = utils.OperationNode(
3526            "delaunay2d",
3527            parents=[self],
3528            comment=f"#cells {msh.dataset.GetNumberOfCells()}",
3529        )
3530        return msh
3531
3532    def generate_voronoi(self, padding=0.0, fit=False, method="vtk") -> "vedo.Mesh":
3533        """
3534        Generate the 2D Voronoi convex tiling of the input points (z is ignored).
3535        The points are assumed to lie in a plane. The output is a Mesh. Each output cell is a convex polygon.
3536
3537        A cell array named "VoronoiID" is added to the output Mesh.
3538
3539        The 2D Voronoi tessellation is a tiling of space, where each Voronoi tile represents the region nearest
3540        to one of the input points. Voronoi tessellations are important in computational geometry
3541        (and many other fields), and are the dual of Delaunay triangulations.
3542
3543        Thus the triangulation is constructed in the x-y plane, and the z coordinate is ignored
3544        (although carried through to the output).
3545        If you desire to triangulate in a different plane, you can use fit=True.
3546
3547        A brief summary is as follows. Each (generating) input point is associated with
3548        an initial Voronoi tile, which is simply the bounding box of the point set.
3549        A locator is then used to identify nearby points: each neighbor in turn generates a
3550        clipping line positioned halfway between the generating point and the neighboring point,
3551        and orthogonal to the line connecting them. Clips are readily performed by evaluationg the
3552        vertices of the convex Voronoi tile as being on either side (inside,outside) of the clip line.
3553        If two intersections of the Voronoi tile are found, the portion of the tile "outside" the clip
3554        line is discarded, resulting in a new convex, Voronoi tile. As each clip occurs,
3555        the Voronoi "Flower" error metric (the union of error spheres) is compared to the extent of the region
3556        containing the neighboring clip points. The clip region (along with the points contained in it) is grown
3557        by careful expansion (e.g., outward spiraling iterator over all candidate clip points).
3558        When the Voronoi Flower is contained within the clip region, the algorithm terminates and the Voronoi
3559        tile is output. Once complete, it is possible to construct the Delaunay triangulation from the Voronoi
3560        tessellation. Note that topological and geometric information is used to generate a valid triangulation
3561        (e.g., merging points and validating topology).
3562
3563        Arguments:
3564            pts : (list)
3565                list of input points.
3566            padding : (float)
3567                padding distance. The default is 0.
3568            fit : (bool)
3569                detect automatically the best fitting plane. The default is False.
3570
3571        Examples:
3572            - [voronoi1.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/voronoi1.py)
3573
3574                ![](https://vedo.embl.es/images/basic/voronoi1.png)
3575
3576            - [voronoi2.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/voronoi2.py)
3577
3578                ![](https://vedo.embl.es/images/advanced/voronoi2.png)
3579        """
3580        pts = self.coordinates
3581
3582        if method == "scipy":
3583            from scipy.spatial import Voronoi as scipy_voronoi
3584
3585            pts = np.asarray(pts)[:, (0, 1)]
3586            vor = scipy_voronoi(pts)
3587            regs = []  # filter out invalid indices
3588            for r in vor.regions:
3589                flag = True
3590                for x in r:
3591                    if x < 0:
3592                        flag = False
3593                        break
3594                if flag and len(r) > 0:
3595                    regs.append(r)
3596
3597            m = vedo.Mesh([vor.vertices, regs])
3598            m.celldata["VoronoiID"] = np.array(list(range(len(regs)))).astype(int)
3599
3600        elif method == "vtk":
3601            vor = vtki.new("Voronoi2D")
3602            if isinstance(pts, Points):
3603                vor.SetInputData(pts)
3604            else:
3605                pts = np.asarray(pts)
3606                if pts.shape[1] == 2:
3607                    pts = np.c_[pts, np.zeros(len(pts))]
3608                pd = vtki.vtkPolyData()
3609                vpts = vtki.vtkPoints()
3610                vpts.SetData(utils.numpy2vtk(pts, dtype=np.float32))
3611                pd.SetPoints(vpts)
3612                vor.SetInputData(pd)
3613            vor.SetPadding(padding)
3614            vor.SetGenerateScalarsToPointIds()
3615            if fit:
3616                vor.SetProjectionPlaneModeToBestFittingPlane()
3617            else:
3618                vor.SetProjectionPlaneModeToXYPlane()
3619            vor.Update()
3620            poly = vor.GetOutput()
3621            arr = poly.GetCellData().GetArray(0)
3622            if arr:
3623                arr.SetName("VoronoiID")
3624            m = vedo.Mesh(poly, c="orange5")
3625
3626        else:
3627            vedo.logger.error(f"Unknown method {method} in voronoi()")
3628            raise RuntimeError
3629
3630        m.lw(2).lighting("off").wireframe()
3631        m.name = "Voronoi"
3632        return m
3633
3634    ##########################################################################
3635    def generate_delaunay3d(self, radius=0, tol=None) -> "vedo.TetMesh":
3636        """
3637        Create 3D Delaunay triangulation of input points.
3638
3639        Arguments:
3640            radius : (float)
3641                specify distance (or "alpha") value to control output.
3642                For a non-zero values, only tetra contained within the circumsphere
3643                will be output.
3644            tol : (float)
3645                Specify a tolerance to control discarding of closely spaced points.
3646                This tolerance is specified as a fraction of the diagonal length of
3647                the bounding box of the points.
3648        """
3649        deln = vtki.new("Delaunay3D")
3650        deln.SetInputData(self.dataset)
3651        deln.SetAlpha(radius)
3652        deln.AlphaTetsOn()
3653        deln.AlphaTrisOff()
3654        deln.AlphaLinesOff()
3655        deln.AlphaVertsOff()
3656        deln.BoundingTriangulationOff()
3657        if tol:
3658            deln.SetTolerance(tol)
3659        deln.Update()
3660        m = vedo.TetMesh(deln.GetOutput())
3661        m.pipeline = utils.OperationNode(
3662            "generate_delaunay3d", c="#e9c46a:#edabab", parents=[self],
3663        )
3664        m.name = "Delaunay3D"
3665        return m
3666
3667    ####################################################
3668    def visible_points(self, area=(), tol=None, invert=False) -> Union[Self, None]:
3669        """
3670        Extract points based on whether they are visible or not.
3671        Visibility is determined by accessing the z-buffer of a rendering window.
3672        The position of each input point is converted into display coordinates,
3673        and then the z-value at that point is obtained.
3674        If within the user-specified tolerance, the point is considered visible.
3675        Associated data attributes are passed to the output as well.
3676
3677        This filter also allows you to specify a rectangular window in display (pixel)
3678        coordinates in which the visible points must lie.
3679
3680        Arguments:
3681            area : (list)
3682                specify a rectangular region as (xmin,xmax,ymin,ymax)
3683            tol : (float)
3684                a tolerance in normalized display coordinate system
3685            invert : (bool)
3686                select invisible points instead.
3687
3688        Example:
3689            ```python
3690            from vedo import Ellipsoid, show
3691            s = Ellipsoid().rotate_y(30)
3692
3693            # Camera options: pos, focal_point, viewup, distance
3694            camopts = dict(pos=(0,0,25), focal_point=(0,0,0))
3695            show(s, camera=camopts, offscreen=True)
3696
3697            m = s.visible_points()
3698            # print('visible pts:', m.vertices)  # numpy array
3699            show(m, new=True, axes=1).close() # optionally draw result in a new window
3700            ```
3701            ![](https://vedo.embl.es/images/feats/visible_points.png)
3702        """
3703        svp = vtki.new("SelectVisiblePoints")
3704        svp.SetInputData(self.dataset)
3705
3706        ren = None
3707        if vedo.plotter_instance:
3708            if vedo.plotter_instance.renderer:
3709                ren = vedo.plotter_instance.renderer
3710                svp.SetRenderer(ren)
3711        if not ren:
3712            vedo.logger.warning(
3713                "visible_points() can only be used after a rendering step"
3714            )
3715            return None
3716
3717        if len(area) == 2:
3718            area = utils.flatten(area)
3719        if len(area) == 4:
3720            # specify a rectangular region
3721            svp.SetSelection(area[0], area[1], area[2], area[3])
3722        if tol is not None:
3723            svp.SetTolerance(tol)
3724        if invert:
3725            svp.SelectInvisibleOn()
3726        svp.Update()
3727
3728        m = Points(svp.GetOutput())
3729        m.name = "VisiblePoints"
3730        return m
 453class Points(PointsVisual, PointAlgorithms):
 454    """Work with point clouds."""
 455
 456    def __init__(self, inputobj=None, r=4, c=(0.2, 0.2, 0.2), alpha=1):
 457        """
 458        Build an object made of only vertex points for a list of 2D/3D points.
 459        Both shapes (N, 3) or (3, N) are accepted as input, if N>3.
 460
 461        Arguments:
 462            inputobj : (list, tuple)
 463            r : (int)
 464                Point radius in units of pixels.
 465            c : (str, list)
 466                Color name or rgb tuple.
 467            alpha : (float)
 468                Transparency in range [0,1].
 469
 470        Example:
 471            ```python
 472            from vedo import *
 473
 474            def fibonacci_sphere(n):
 475                s = np.linspace(0, n, num=n, endpoint=False)
 476                theta = s * 2.399963229728653
 477                y = 1 - s * (2/(n-1))
 478                r = np.sqrt(1 - y * y)
 479                x = np.cos(theta) * r
 480                z = np.sin(theta) * r
 481                return np._c[x,y,z]
 482
 483            Points(fibonacci_sphere(1000)).show(axes=1).close()
 484            ```
 485            ![](https://vedo.embl.es/images/feats/fibonacci.png)
 486        """
 487        # print("INIT POINTS")
 488        super().__init__()
 489
 490        self.name = ""
 491        self.filename = ""
 492        self.file_size = ""
 493
 494        self.info = {}
 495        self.time = time.time()
 496        
 497        self.transform = LinearTransform()
 498        self.point_locator = None
 499        self.cell_locator = None
 500        self.line_locator = None
 501
 502        self.actor = vtki.vtkActor()
 503        self.properties = self.actor.GetProperty()
 504        self.properties_backface = self.actor.GetBackfaceProperty()
 505        self.mapper = vtki.new("PolyDataMapper")
 506        self.dataset = vtki.vtkPolyData()
 507        
 508        # Create weakref so actor can access this object (eg to pick/remove):
 509        self.actor.retrieve_object = weak_ref_to(self)
 510
 511        try:
 512            self.properties.RenderPointsAsSpheresOn()
 513        except AttributeError:
 514            pass
 515
 516        if inputobj is None:  ####################
 517            return
 518        ##########################################
 519
 520        self.name = "Points"
 521
 522        ######
 523        if isinstance(inputobj, vtki.vtkActor):
 524            self.dataset.DeepCopy(inputobj.GetMapper().GetInput())
 525            pr = vtki.vtkProperty()
 526            pr.DeepCopy(inputobj.GetProperty())
 527            self.actor.SetProperty(pr)
 528            self.properties = pr
 529            self.mapper.SetScalarVisibility(inputobj.GetMapper().GetScalarVisibility())
 530
 531        elif isinstance(inputobj, vtki.vtkPolyData):
 532            self.dataset = inputobj
 533            if self.dataset.GetNumberOfCells() == 0:
 534                carr = vtki.vtkCellArray()
 535                for i in range(self.dataset.GetNumberOfPoints()):
 536                    carr.InsertNextCell(1)
 537                    carr.InsertCellPoint(i)
 538                self.dataset.SetVerts(carr)
 539
 540        elif isinstance(inputobj, Points):
 541            self.dataset = inputobj.dataset
 542            self.copy_properties_from(inputobj)
 543
 544        elif utils.is_sequence(inputobj):  # passing point coords
 545            self.dataset = utils.buildPolyData(utils.make3d(inputobj))
 546
 547        elif isinstance(inputobj, str):
 548            verts = vedo.file_io.load(inputobj)
 549            self.filename = inputobj
 550            self.dataset = verts.dataset
 551
 552        elif "meshlib" in str(type(inputobj)):
 553            from meshlib import mrmeshnumpy as mn
 554            self.dataset = utils.buildPolyData(mn.toNumpyArray(inputobj.points))
 555
 556        else:
 557            # try to extract the points from a generic VTK input data object
 558            if hasattr(inputobj, "dataset"):
 559                inputobj = inputobj.dataset
 560            try:
 561                vvpts = inputobj.GetPoints()
 562                self.dataset = vtki.vtkPolyData()
 563                self.dataset.SetPoints(vvpts)
 564                for i in range(inputobj.GetPointData().GetNumberOfArrays()):
 565                    arr = inputobj.GetPointData().GetArray(i)
 566                    self.dataset.GetPointData().AddArray(arr)
 567                carr = vtki.vtkCellArray()
 568                for i in range(self.dataset.GetNumberOfPoints()):
 569                    carr.InsertNextCell(1)
 570                    carr.InsertCellPoint(i)
 571                self.dataset.SetVerts(carr)
 572            except:
 573                vedo.logger.error(f"cannot build Points from type {type(inputobj)}")
 574                raise RuntimeError()
 575
 576        self.actor.SetMapper(self.mapper)
 577        self.mapper.SetInputData(self.dataset)
 578
 579        self.properties.SetColor(colors.get_color(c))
 580        self.properties.SetOpacity(alpha)
 581        self.properties.SetRepresentationToPoints()
 582        self.properties.SetPointSize(r)
 583        self.properties.LightingOff()
 584
 585        self.pipeline = utils.OperationNode(
 586            self, parents=[], comment=f"#pts {self.dataset.GetNumberOfPoints()}"
 587        )
 588
 589    def _update(self, polydata, reset_locators=True) -> Self:
 590        """Overwrite the polygonal dataset with a new vtkPolyData."""
 591        self.dataset = polydata
 592        self.mapper.SetInputData(self.dataset)
 593        self.mapper.Modified()
 594        if reset_locators:
 595            self.point_locator = None
 596            self.line_locator = None
 597            self.cell_locator = None
 598        return self
 599
 600    def __str__(self):
 601        """Print a description of the Points/Mesh."""
 602        module = self.__class__.__module__
 603        name = self.__class__.__name__
 604        out = vedo.printc(
 605            f"{module}.{name} at ({hex(self.memory_address())})".ljust(75),
 606            c="g", bold=True, invert=True, return_string=True,
 607        )
 608        out += "\x1b[0m\x1b[32;1m"
 609
 610        if self.name:
 611            out += "name".ljust(14) + ": " + self.name
 612            if "legend" in self.info.keys() and self.info["legend"]:
 613                out+= f", legend='{self.info['legend']}'"
 614            out += "\n"
 615 
 616        if self.filename:
 617            out+= "file name".ljust(14) + ": " + self.filename + "\n"
 618
 619        if not self.mapper.GetScalarVisibility():
 620            col = utils.precision(self.properties.GetColor(), 3)
 621            cname = vedo.colors.get_color_name(self.properties.GetColor())
 622            out+= "color".ljust(14) + ": " + cname 
 623            out+= f", rgb={col}, alpha={self.properties.GetOpacity()}\n"
 624            if self.actor.GetBackfaceProperty():
 625                bcol = self.actor.GetBackfaceProperty().GetDiffuseColor()
 626                cname = vedo.colors.get_color_name(bcol)
 627                out+= "backface color".ljust(14) + ": " 
 628                out+= f"{cname}, rgb={utils.precision(bcol,3)}\n"
 629
 630        npt = self.dataset.GetNumberOfPoints()
 631        npo, nln = self.dataset.GetNumberOfPolys(), self.dataset.GetNumberOfLines()
 632        out+= "elements".ljust(14) + f": vertices={npt:,} polygons={npo:,} lines={nln:,}"
 633        if self.dataset.GetNumberOfStrips():
 634            out+= f", strips={self.dataset.GetNumberOfStrips():,}"
 635        out+= "\n"
 636        if self.dataset.GetNumberOfPieces() > 1:
 637            out+= "pieces".ljust(14) + ": " + str(self.dataset.GetNumberOfPieces()) + "\n"
 638
 639        out+= "position".ljust(14) + ": " + f"{utils.precision(self.pos(), 6)}\n"
 640        try:
 641            sc = self.transform.get_scale()
 642            out+= "scaling".ljust(14)  + ": "
 643            out+= utils.precision(sc, 6) + "\n"
 644        except AttributeError:
 645            pass
 646
 647        if self.npoints:
 648            out+="size".ljust(14)+ ": average=" + utils.precision(self.average_size(),6)
 649            out+=", diagonal="+ utils.precision(self.diagonal_size(), 6)+ "\n"
 650            out+="center of mass".ljust(14) + ": " + utils.precision(self.center_of_mass(),6)+"\n"
 651
 652        bnds = self.bounds()
 653        bx1, bx2 = utils.precision(bnds[0], 3), utils.precision(bnds[1], 3)
 654        by1, by2 = utils.precision(bnds[2], 3), utils.precision(bnds[3], 3)
 655        bz1, bz2 = utils.precision(bnds[4], 3), utils.precision(bnds[5], 3)
 656        out+= "bounds".ljust(14) + ":"
 657        out+= " x=(" + bx1 + ", " + bx2 + "),"
 658        out+= " y=(" + by1 + ", " + by2 + "),"
 659        out+= " z=(" + bz1 + ", " + bz2 + ")\n"
 660
 661        for key in self.pointdata.keys():
 662            arr = self.pointdata[key]
 663            dim = arr.shape[1] if arr.ndim > 1 else 1
 664            mark_active = "pointdata"
 665            a_scalars = self.dataset.GetPointData().GetScalars()
 666            a_vectors = self.dataset.GetPointData().GetVectors()
 667            a_tensors = self.dataset.GetPointData().GetTensors()
 668            if   a_scalars and a_scalars.GetName() == key:
 669                mark_active += " *"
 670            elif a_vectors and a_vectors.GetName() == key:
 671                mark_active += " **"
 672            elif a_tensors and a_tensors.GetName() == key:
 673                mark_active += " ***"
 674            out += mark_active.ljust(14) + f': "{key}" ({arr.dtype}), dim={dim}'
 675            if dim == 1 and len(arr):
 676                rng = utils.precision(arr.min(), 3) + ", " + utils.precision(arr.max(), 3)
 677                out += f", range=({rng})\n"
 678            else:
 679                out += "\n"
 680
 681        for key in self.celldata.keys():
 682            arr = self.celldata[key]
 683            dim = arr.shape[1] if arr.ndim > 1 else 1
 684            mark_active = "celldata"
 685            a_scalars = self.dataset.GetCellData().GetScalars()
 686            a_vectors = self.dataset.GetCellData().GetVectors()
 687            a_tensors = self.dataset.GetCellData().GetTensors()
 688            if   a_scalars and a_scalars.GetName() == key:
 689                mark_active += " *"
 690            elif a_vectors and a_vectors.GetName() == key:
 691                mark_active += " **"
 692            elif a_tensors and a_tensors.GetName() == key:
 693                mark_active += " ***"
 694            out += mark_active.ljust(14) + f': "{key}" ({arr.dtype}), dim={dim}'
 695            if dim == 1 and len(arr):
 696                rng = utils.precision(arr.min(), 3) + ", " + utils.precision(arr.max(), 3)
 697                out += f", range=({rng})\n"
 698            else:
 699                out += "\n"
 700
 701        for key in self.metadata.keys():
 702            arr = self.metadata[key]
 703            if len(arr) > 3:
 704                out+= "metadata".ljust(14) + ": " + f'"{key}" ({len(arr)} values)\n'
 705            else:
 706                out+= "metadata".ljust(14) + ": " + f'"{key}" = {arr}\n'
 707
 708        if self.picked3d is not None:
 709            idp = self.closest_point(self.picked3d, return_point_id=True)
 710            idc = self.closest_point(self.picked3d, return_cell_id=True)
 711            out+= "clicked point".ljust(14) + ": " + utils.precision(self.picked3d, 6)
 712            out+= f", pointID={idp}, cellID={idc}\n"
 713
 714        return out.rstrip() + "\x1b[0m"
 715
 716    def _repr_html_(self):
 717        """
 718        HTML representation of the Point cloud object for Jupyter Notebooks.
 719
 720        Returns:
 721            HTML text with the image and some properties.
 722        """
 723        import io
 724        import base64
 725        from PIL import Image
 726
 727        library_name = "vedo.pointcloud.Points"
 728        help_url = "https://vedo.embl.es/docs/vedo/pointcloud.html#Points"
 729
 730        arr = self.thumbnail()
 731        im = Image.fromarray(arr)
 732        buffered = io.BytesIO()
 733        im.save(buffered, format="PNG", quality=100)
 734        encoded = base64.b64encode(buffered.getvalue()).decode("utf-8")
 735        url = "data:image/png;base64," + encoded
 736        image = f"<img src='{url}'></img>"
 737
 738        bounds = "<br/>".join(
 739            [
 740                utils.precision(min_x, 4) + " ... " + utils.precision(max_x, 4)
 741                for min_x, max_x in zip(self.bounds()[::2], self.bounds()[1::2])
 742            ]
 743        )
 744        average_size = "{size:.3f}".format(size=self.average_size())
 745
 746        help_text = ""
 747        if self.name:
 748            help_text += f"<b> {self.name}: &nbsp&nbsp</b>"
 749        help_text += '<b><a href="' + help_url + '" target="_blank">' + library_name + "</a></b>"
 750        if self.filename:
 751            dots = ""
 752            if len(self.filename) > 30:
 753                dots = "..."
 754            help_text += f"<br/><code><i>({dots}{self.filename[-30:]})</i></code>"
 755
 756        pdata = ""
 757        if self.dataset.GetPointData().GetScalars():
 758            if self.dataset.GetPointData().GetScalars().GetName():
 759                name = self.dataset.GetPointData().GetScalars().GetName()
 760                pdata = "<tr><td><b> point data array </b></td><td>" + name + "</td></tr>"
 761
 762        cdata = ""
 763        if self.dataset.GetCellData().GetScalars():
 764            if self.dataset.GetCellData().GetScalars().GetName():
 765                name = self.dataset.GetCellData().GetScalars().GetName()
 766                cdata = "<tr><td><b> cell data array </b></td><td>" + name + "</td></tr>"
 767
 768        allt = [
 769            "<table>",
 770            "<tr>",
 771            "<td>",
 772            image,
 773            "</td>",
 774            "<td style='text-align: center; vertical-align: center;'><br/>",
 775            help_text,
 776            "<table>",
 777            "<tr><td><b> bounds </b> <br/> (x/y/z) </td><td>" + str(bounds) + "</td></tr>",
 778            "<tr><td><b> center of mass </b></td><td>"
 779            + utils.precision(self.center_of_mass(), 3)
 780            + "</td></tr>",
 781            "<tr><td><b> average size </b></td><td>" + str(average_size) + "</td></tr>",
 782            "<tr><td><b> nr. points </b></td><td>" + str(self.npoints) + "</td></tr>",
 783            pdata,
 784            cdata,
 785            "</table>",
 786            "</table>",
 787        ]
 788        return "\n".join(allt)
 789
 790    ##################################################################################
 791    def __add__(self, meshs):
 792        """
 793        Add two meshes or a list of meshes together to form an `Assembly` object.
 794        """
 795        if isinstance(meshs, list):
 796            alist = [self]
 797            for l in meshs:
 798                if isinstance(l, vedo.Assembly):
 799                    alist += l.unpack()
 800                else:
 801                    alist += l
 802            return vedo.assembly.Assembly(alist)
 803
 804        if isinstance(meshs, vedo.Assembly):
 805            return meshs + self  # use Assembly.__add__
 806
 807        return vedo.assembly.Assembly([self, meshs])
 808
 809    def polydata(self, **kwargs):
 810        """
 811        Obsolete. Use property `.dataset` instead.
 812        Returns the underlying `vtkPolyData` object.
 813        """
 814        colors.printc(
 815            "WARNING: call to .polydata() is obsolete, use property .dataset instead.",
 816            c="y")
 817        return self.dataset
 818
 819    def __copy__(self):
 820        return self.clone(deep=False)
 821
 822    def __deepcopy__(self, memo):
 823        return self.clone(deep=memo)
 824    
 825    def copy(self, deep=True) -> Self:
 826        """Return a copy of the object. Alias of `clone()`."""
 827        return self.clone(deep=deep)
 828
 829    def clone(self, deep=True) -> Self:
 830        """
 831        Clone a `PointCloud` or `Mesh` object to make an exact copy of it.
 832        Alias of `copy()`.
 833
 834        Arguments:
 835            deep : (bool)
 836                if False return a shallow copy of the mesh without copying the points array.
 837
 838        Examples:
 839            - [mirror.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/mirror.py)
 840
 841               ![](https://vedo.embl.es/images/basic/mirror.png)
 842        """
 843        poly = vtki.vtkPolyData()
 844        if deep or isinstance(deep, dict): # if a memo object is passed this checks as True
 845            poly.DeepCopy(self.dataset)
 846        else:
 847            poly.ShallowCopy(self.dataset)
 848
 849        if isinstance(self, vedo.Mesh):
 850            cloned = vedo.Mesh(poly)
 851        else:
 852            cloned = Points(poly)
 853        # print([self], self.__class__)
 854        # cloned = self.__class__(poly)
 855
 856        cloned.transform = self.transform.clone()
 857
 858        cloned.copy_properties_from(self)
 859
 860        cloned.name = str(self.name)
 861        cloned.filename = str(self.filename)
 862        cloned.info = dict(self.info)
 863        cloned.pipeline = utils.OperationNode("clone", parents=[self], shape="diamond", c="#edede9")
 864
 865        if isinstance(deep, dict):
 866            deep[id(self)] = cloned
 867
 868        return cloned
 869
 870    def compute_normals_with_pca(self, n=20, orientation_point=None, invert=False) -> Self:
 871        """
 872        Generate point normals using PCA (principal component analysis).
 873        This algorithm estimates a local tangent plane around each sample point p
 874        by considering a small neighborhood of points around p, and fitting a plane
 875        to the neighborhood (via PCA).
 876
 877        Arguments:
 878            n : (int)
 879                neighborhood size to calculate the normal
 880            orientation_point : (list)
 881                adjust the +/- sign of the normals so that
 882                the normals all point towards a specified point. If None, perform a traversal
 883                of the point cloud and flip neighboring normals so that they are mutually consistent.
 884            invert : (bool)
 885                flip all normals
 886        """
 887        poly = self.dataset
 888        pcan = vtki.new("PCANormalEstimation")
 889        pcan.SetInputData(poly)
 890        pcan.SetSampleSize(n)
 891
 892        if orientation_point is not None:
 893            pcan.SetNormalOrientationToPoint()
 894            pcan.SetOrientationPoint(orientation_point)
 895        else:
 896            pcan.SetNormalOrientationToGraphTraversal()
 897
 898        if invert:
 899            pcan.FlipNormalsOn()
 900        pcan.Update()
 901
 902        varr = pcan.GetOutput().GetPointData().GetNormals()
 903        varr.SetName("Normals")
 904        self.dataset.GetPointData().SetNormals(varr)
 905        self.dataset.GetPointData().Modified()
 906        return self
 907
 908    def compute_acoplanarity(self, n=25, radius=None, on="points") -> Self:
 909        """
 910        Compute acoplanarity which is a measure of how much a local region of the mesh
 911        differs from a plane.
 912        
 913        The information is stored in a `pointdata` or `celldata` array with name 'Acoplanarity'.
 914        
 915        Either `n` (number of neighbour points) or `radius` (radius of local search) can be specified.
 916        If a radius value is given and not enough points fall inside it, then a -1 is stored.
 917
 918        Example:
 919            ```python
 920            from vedo import *
 921            msh = ParametricShape('RandomHills')
 922            msh.compute_acoplanarity(radius=0.1, on='cells')
 923            msh.cmap("coolwarm", on='cells').add_scalarbar()
 924            msh.show(axes=1).close()
 925            ```
 926            ![](https://vedo.embl.es/images/feats/acoplanarity.jpg)
 927        """
 928        acoplanarities = []
 929        if "point" in on:
 930            pts = self.coordinates
 931        elif "cell" in on:
 932            pts = self.cell_centers().coordinates
 933        else:
 934            raise ValueError(f"In compute_acoplanarity() set on to either 'cells' or 'points', not {on}")
 935
 936        for p in utils.progressbar(pts, delay=5, width=15, title=f"{on} acoplanarity"):
 937            if n:
 938                data = self.closest_point(p, n=n)
 939                npts = n
 940            elif radius:
 941                data = self.closest_point(p, radius=radius)
 942                npts = len(data)
 943
 944            try:
 945                center = data.mean(axis=0)
 946                res = np.linalg.svd(data - center)
 947                acoplanarities.append(res[1][2] / npts)
 948            except:
 949                acoplanarities.append(-1.0)
 950
 951        if "point" in on:
 952            self.pointdata["Acoplanarity"] = np.array(acoplanarities, dtype=float)
 953        else:
 954            self.celldata["Acoplanarity"] = np.array(acoplanarities, dtype=float)
 955        return self
 956
 957    def distance_to(self, pcloud, signed=False, invert=False, name="Distance") -> np.ndarray:
 958        """
 959        Computes the distance from one point cloud or mesh to another point cloud or mesh.
 960        This new `pointdata` array is saved with default name "Distance".
 961
 962        Keywords `signed` and `invert` are used to compute signed distance,
 963        but the mesh in that case must have polygonal faces (not a simple point cloud),
 964        and normals must also be computed.
 965
 966        Examples:
 967            - [distance2mesh.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/distance2mesh.py)
 968
 969                ![](https://vedo.embl.es/images/basic/distance2mesh.png)
 970        """
 971        if pcloud.dataset.GetNumberOfPolys():
 972
 973            poly1 = self.dataset
 974            poly2 = pcloud.dataset
 975            df = vtki.new("DistancePolyDataFilter")
 976            df.ComputeSecondDistanceOff()
 977            df.SetInputData(0, poly1)
 978            df.SetInputData(1, poly2)
 979            df.SetSignedDistance(signed)
 980            df.SetNegateDistance(invert)
 981            df.Update()
 982            scals = df.GetOutput().GetPointData().GetScalars()
 983            dists = utils.vtk2numpy(scals)
 984
 985        else:  # has no polygons
 986
 987            if signed:
 988                vedo.logger.warning("distance_to() called with signed=True but input object has no polygons")
 989
 990            if not pcloud.point_locator:
 991                pcloud.point_locator = vtki.new("PointLocator")
 992                pcloud.point_locator.SetDataSet(pcloud.dataset)
 993                pcloud.point_locator.BuildLocator()
 994
 995            ids = []
 996            ps1 = self.coordinates
 997            ps2 = pcloud.coordinates
 998            for p in ps1:
 999                pid = pcloud.point_locator.FindClosestPoint(p)
1000                ids.append(pid)
1001
1002            deltas = ps2[ids] - ps1
1003            dists = np.linalg.norm(deltas, axis=1).astype(np.float32)
1004            scals = utils.numpy2vtk(dists)
1005
1006        scals.SetName(name)
1007        self.dataset.GetPointData().AddArray(scals)
1008        self.dataset.GetPointData().SetActiveScalars(scals.GetName())
1009        rng = scals.GetRange()
1010        self.mapper.SetScalarRange(rng[0], rng[1])
1011        self.mapper.ScalarVisibilityOn()
1012
1013        self.pipeline = utils.OperationNode(
1014            "distance_to",
1015            parents=[self, pcloud],
1016            shape="cylinder",
1017            comment=f"#pts {self.dataset.GetNumberOfPoints()}",
1018        )
1019        return dists
1020
1021    def clean(self) -> Self:
1022        """Clean pointcloud or mesh by removing coincident points."""
1023        cpd = vtki.new("CleanPolyData")
1024        cpd.PointMergingOn()
1025        cpd.ConvertLinesToPointsOff()
1026        cpd.ConvertPolysToLinesOff()
1027        cpd.ConvertStripsToPolysOff()
1028        cpd.SetInputData(self.dataset)
1029        cpd.Update()
1030        self._update(cpd.GetOutput())
1031        self.pipeline = utils.OperationNode(
1032            "clean", parents=[self], comment=f"#pts {self.dataset.GetNumberOfPoints()}"
1033        )
1034        return self
1035
1036    def subsample(self, fraction: float, absolute=False) -> Self:
1037        """
1038        Subsample a point cloud by requiring that the points
1039        or vertices are far apart at least by the specified fraction of the object size.
1040        If a Mesh is passed the polygonal faces are not removed
1041        but holes can appear as their vertices are removed.
1042
1043        Examples:
1044            - [moving_least_squares1D.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/moving_least_squares1D.py)
1045
1046                ![](https://vedo.embl.es/images/advanced/moving_least_squares1D.png)
1047
1048            - [recosurface.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/recosurface.py)
1049
1050                ![](https://vedo.embl.es/images/advanced/recosurface.png)
1051        """
1052        if not absolute:
1053            if fraction > 1:
1054                vedo.logger.warning(
1055                    f"subsample(fraction=...), fraction must be < 1, but is {fraction}"
1056                )
1057            if fraction <= 0:
1058                return self
1059
1060        cpd = vtki.new("CleanPolyData")
1061        cpd.PointMergingOn()
1062        cpd.ConvertLinesToPointsOn()
1063        cpd.ConvertPolysToLinesOn()
1064        cpd.ConvertStripsToPolysOn()
1065        cpd.SetInputData(self.dataset)
1066        if absolute:
1067            cpd.SetTolerance(fraction / self.diagonal_size())
1068            # cpd.SetToleranceIsAbsolute(absolute)
1069        else:
1070            cpd.SetTolerance(fraction)
1071        cpd.Update()
1072
1073        ps = 2
1074        if self.properties.GetRepresentation() == 0:
1075            ps = self.properties.GetPointSize()
1076
1077        self._update(cpd.GetOutput())
1078        self.ps(ps)
1079
1080        self.pipeline = utils.OperationNode(
1081            "subsample", parents=[self], comment=f"#pts {self.dataset.GetNumberOfPoints()}"
1082        )
1083        return self
1084
1085    def threshold(self, scalars: str, above=None, below=None, on="points") -> Self:
1086        """
1087        Extracts cells where scalar value satisfies threshold criterion.
1088
1089        Arguments:
1090            scalars : (str)
1091                name of the scalars array.
1092            above : (float)
1093                minimum value of the scalar
1094            below : (float)
1095                maximum value of the scalar
1096            on : (str)
1097                if 'cells' assume array of scalars refers to cell data.
1098
1099        Examples:
1100            - [mesh_threshold.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/mesh_threshold.py)
1101        """
1102        thres = vtki.new("Threshold")
1103        thres.SetInputData(self.dataset)
1104
1105        if on.startswith("c"):
1106            asso = vtki.vtkDataObject.FIELD_ASSOCIATION_CELLS
1107        else:
1108            asso = vtki.vtkDataObject.FIELD_ASSOCIATION_POINTS
1109
1110        thres.SetInputArrayToProcess(0, 0, 0, asso, scalars)
1111
1112        if above is None and below is not None:
1113            try:  # vtk 9.2
1114                thres.ThresholdByLower(below)
1115            except AttributeError:  # vtk 9.3
1116                thres.SetUpperThreshold(below)
1117