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