vedo.pointcloud

Submodule to work with point clouds

   1#!/usr/bin/env python3
   2# -*- coding: utf-8 -*-
   3from deprecated import deprecated
   4import numpy as np
   5
   6try:
   7    import vedo.vtkclasses as vtk
   8except ImportError:
   9    import vtkmodules.all as vtk
  10
  11import vedo
  12from vedo import colors
  13from vedo import utils
  14from vedo.base import BaseActor
  15
  16__docformat__ = "google"
  17
  18__doc__ = """
  19Submodule to work with point clouds <br>
  20
  21![](https://vedo.embl.es/images/basic/pca.png)
  22"""
  23
  24__all__ = [
  25    "Points",
  26    "Point",
  27    "merge",
  28    "visible_points",
  29    "delaunay2d",
  30    "voronoi",
  31    "fit_line",
  32    "fit_circle",
  33    "fit_plane",
  34    "fit_sphere",
  35    "pca_ellipse",
  36    "pca_ellipsoid",
  37    "pcaEllipsoid",  # deprecated
  38]
  39
  40
  41####################################################
  42def merge(*meshs, flag=False):
  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 use `flag`.
  49    In this case a point array of IDs is added to the output.
  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    acts = [a for a in utils.flatten(meshs) if a]
  60
  61    if not acts:
  62        return None
  63
  64    idarr = []
  65    polyapp = vtk.vtkAppendPolyData()
  66    for i, a in enumerate(acts):
  67        try:
  68            poly = a.polydata()
  69        except AttributeError:
  70            # so a vtkPolydata can also be passed
  71            poly = a
  72        polyapp.AddInputData(poly)
  73        if flag:
  74            idarr += [i] * poly.GetNumberOfPoints()
  75    polyapp.Update()
  76    mpoly = polyapp.GetOutput()
  77
  78    if flag:
  79        varr = utils.numpy2vtk(idarr, dtype=np.uint16, name="OriginalMeshID")
  80        mpoly.GetPointData().AddArray(varr)
  81
  82    if isinstance(acts[0], vedo.Mesh):
  83        msh = vedo.Mesh(mpoly)
  84    else:
  85        msh = Points(mpoly)
  86
  87    if isinstance(acts[0], vtk.vtkActor):
  88        cprp = vtk.vtkProperty()
  89        cprp.DeepCopy(acts[0].GetProperty())
  90        msh.SetProperty(cprp)
  91        msh.property = cprp
  92
  93    msh.pipeline = utils.OperationNode(
  94        "merge", parents=acts,
  95        comment=f"#pts {msh._data.GetNumberOfPoints()}",
  96    )
  97    return msh
  98
  99
 100####################################################
 101def visible_points(mesh, area=(), tol=None, invert=False):
 102    """
 103    Extract points based on whether they are visible or not.
 104    Visibility is determined by accessing the z-buffer of a rendering window.
 105    The position of each input point is converted into display coordinates,
 106    and then the z-value at that point is obtained.
 107    If within the user-specified tolerance, the point is considered visible.
 108    Associated data attributes are passed to the output as well.
 109
 110    This filter also allows you to specify a rectangular window in display (pixel)
 111    coordinates in which the visible points must lie.
 112
 113    Arguments:
 114        area : (list)
 115            specify a rectangular region as (xmin,xmax,ymin,ymax)
 116        tol : (float)
 117            a tolerance in normalized display coordinate system
 118        invert : (bool)
 119            select invisible points instead.
 120
 121    Example:
 122        ```python
 123        from vedo import Ellipsoid, show, visible_points
 124        s = Ellipsoid().rotate_y(30)
 125
 126        #Camera options: pos, focal_point, viewup, distance,
 127        camopts = dict(pos=(0,0,25), focal_point=(0,0,0))
 128        show(s, camera=camopts, offscreen=True)
 129
 130        m = visible_points(s)
 131        #print('visible pts:', m.points()) # numpy array
 132        show(m, new=True, axes=1) # optionally draw result on a new window
 133        ```
 134        ![](https://vedo.embl.es/images/feats/visible_points.png)
 135    """
 136    # specify a rectangular region
 137    svp = vtk.vtkSelectVisiblePoints()
 138    svp.SetInputData(mesh.polydata())
 139    svp.SetRenderer(vedo.plotter_instance.renderer)
 140
 141    if len(area) == 4:
 142        svp.SetSelection(area[0], area[1], area[2], area[3])
 143    if tol is not None:
 144        svp.SetTolerance(tol)
 145    if invert:
 146        svp.SelectInvisibleOn()
 147    svp.Update()
 148
 149    m = Points(svp.GetOutput()).pointSize(5)
 150    m.name = "VisiblePoints"
 151    return m
 152
 153
 154def delaunay2d(plist, mode="scipy", boundaries=(), tol=None, alpha=0, offset=0, transform=None):
 155    """
 156    Create a mesh from points in the XY plane.
 157    If `mode='fit'` then the filter computes a best fitting
 158    plane and projects the points onto it.
 159
 160    Arguments:
 161        tol : (float)
 162            specify a tolerance to control discarding of closely spaced points.
 163            This tolerance is specified as a fraction of the diagonal length of the bounding box of the points.
 164        alpha : (float)
 165            for a non-zero alpha value, only edges or triangles contained
 166            within a sphere centered at mesh vertices will be output.
 167            Otherwise, only triangles will be output.
 168        offset : (float)
 169            multiplier to control the size of the initial, bounding Delaunay triangulation.
 170        transform: vtkTransform
 171            a VTK transformation (eg. a thinplate spline)
 172            which is applied to points to generate a 2D problem.
 173            This maps a 3D dataset into a 2D dataset where triangulation can be done on the XY plane.
 174            The points are transformed and triangulated.
 175            The topology of triangulated points is used as the output topology.
 176
 177    Examples:
 178        - [delaunay2d.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/delaunay2d.py)
 179
 180            ![](https://vedo.embl.es/images/basic/delaunay2d.png)
 181    """
 182    if isinstance(plist, Points):
 183        parents = [plist]
 184        plist = plist.points()
 185    else:
 186        parents = []
 187        plist = np.ascontiguousarray(plist)
 188        plist = utils.make3d(plist)
 189
 190    #############################################
 191    if mode == "scipy":
 192        from scipy.spatial import Delaunay as scipy_delaunay
 193
 194        tri = scipy_delaunay(plist[:, 0:2])
 195        return vedo.mesh.Mesh([plist, tri.simplices])
 196    #############################################
 197
 198    pd = vtk.vtkPolyData()
 199    vpts = vtk.vtkPoints()
 200    vpts.SetData(utils.numpy2vtk(plist, dtype=np.float32))
 201    pd.SetPoints(vpts)
 202
 203    delny = vtk.vtkDelaunay2D()
 204    delny.SetInputData(pd)
 205    if tol:
 206        delny.SetTolerance(tol)
 207    delny.SetAlpha(alpha)
 208    delny.SetOffset(offset)
 209    if transform:
 210        if hasattr(transform, "transform"):
 211            transform = transform.transform
 212        delny.SetTransform(transform)
 213
 214    if mode == "xy" and boundaries:
 215        boundary = vtk.vtkPolyData()
 216        boundary.SetPoints(vpts)
 217        cell_array = vtk.vtkCellArray()
 218        for b in boundaries:
 219            cpolygon = vtk.vtkPolygon()
 220            for idd in b:
 221                cpolygon.GetPointIds().InsertNextId(idd)
 222            cell_array.InsertNextCell(cpolygon)
 223        boundary.SetPolys(cell_array)
 224        delny.SetSourceData(boundary)
 225
 226    if mode == "fit":
 227        delny.SetProjectionPlaneMode(vtk.VTK_BEST_FITTING_PLANE)
 228    delny.Update()
 229    msh = vedo.mesh.Mesh(delny.GetOutput()).clean().lighting("off")
 230
 231    msh.pipeline = utils.OperationNode(
 232        "delaunay2d", parents=parents, 
 233        comment=f"#cells {msh._data.GetNumberOfCells()}",
 234    )
 235    return msh
 236
 237
 238def voronoi(pts, padding=0, fit=False, method="vtk"):
 239    """
 240    Generate the 2D Voronoi convex tiling of the input points (z is ignored).
 241    The points are assumed to lie in a plane. The output is a Mesh. Each output cell is a convex polygon.
 242
 243    The 2D Voronoi tessellation is a tiling of space, where each Voronoi tile represents the region nearest
 244    to one of the input points. Voronoi tessellations are important in computational geometry
 245    (and many other fields), and are the dual of Delaunay triangulations.
 246
 247    Thus the triangulation is constructed in the x-y plane, and the z coordinate is ignored
 248    (although carried through to the output).
 249    If you desire to triangulate in a different plane, you can use fit=True.
 250
 251    A brief summary is as follows. Each (generating) input point is associated with
 252    an initial Voronoi tile, which is simply the bounding box of the point set.
 253    A locator is then used to identify nearby points: each neighbor in turn generates a
 254    clipping line positioned halfway between the generating point and the neighboring point,
 255    and orthogonal to the line connecting them. Clips are readily performed by evaluationg the
 256    vertices of the convex Voronoi tile as being on either side (inside,outside) of the clip line.
 257    If two intersections of the Voronoi tile are found, the portion of the tile "outside" the clip
 258    line is discarded, resulting in a new convex, Voronoi tile. As each clip occurs,
 259    the Voronoi "Flower" error metric (the union of error spheres) is compared to the extent of the region
 260    containing the neighboring clip points. The clip region (along with the points contained in it) is grown
 261    by careful expansion (e.g., outward spiraling iterator over all candidate clip points).
 262    When the Voronoi Flower is contained within the clip region, the algorithm terminates and the Voronoi
 263    tile is output. Once complete, it is possible to construct the Delaunay triangulation from the Voronoi
 264    tessellation. Note that topological and geometric information is used to generate a valid triangulation
 265    (e.g., merging points and validating topology).
 266
 267    Arguments:
 268        pts : (list)
 269            list of input points.
 270        padding : (float)
 271            padding distance. The default is 0.
 272        fit : (bool)
 273            detect automatically the best fitting plane. The default is False.
 274
 275    Examples:
 276        - [voronoi1.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/voronoi1.py)
 277
 278            ![](https://vedo.embl.es/images/basic/voronoi1.png)
 279
 280        - [voronoi2.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/voronoi2.py)
 281
 282            ![](https://vedo.embl.es/images/advanced/voronoi2.png)
 283    """
 284    if method == "scipy":
 285        from scipy.spatial import Voronoi as scipy_voronoi
 286
 287        pts = np.asarray(pts)[:, (0, 1)]
 288        vor = scipy_voronoi(pts)
 289        regs = []  # filter out invalid indices
 290        for r in vor.regions:
 291            flag = True
 292            for x in r:
 293                if x < 0:
 294                    flag = False
 295                    break
 296            if flag and len(r):
 297                regs.append(r)
 298
 299        m = vedo.Mesh([vor.vertices, regs], c="orange5")
 300        m.celldata["VoronoiID"] = np.array(list(range(len(regs)))).astype(int)
 301        m.locator = None
 302
 303    elif method == "vtk":
 304        vor = vtk.vtkVoronoi2D()
 305        if isinstance(pts, Points):
 306            vor.SetInputData(pts.polydata())
 307        else:
 308            pts = np.asarray(pts)
 309            if pts.shape[1] == 2:
 310                pts = np.c_[pts, np.zeros(len(pts))]
 311            pd = vtk.vtkPolyData()
 312            vpts = vtk.vtkPoints()
 313            vpts.SetData(utils.numpy2vtk(pts, dtype=np.float32))
 314            pd.SetPoints(vpts)
 315            vor.SetInputData(pd)
 316        vor.SetPadding(padding)
 317        vor.SetGenerateScalarsToPointIds()
 318        if fit:
 319            vor.SetProjectionPlaneModeToBestFittingPlane()
 320        else:
 321            vor.SetProjectionPlaneModeToXYPlane()
 322        vor.Update()
 323        poly = vor.GetOutput()
 324        arr = poly.GetCellData().GetArray(0)
 325        if arr:
 326            arr.SetName("VoronoiID")
 327        m = vedo.Mesh(poly, c="orange5")
 328        m.locator = vor.GetLocator()
 329
 330    else:
 331        vedo.logger.error(f"Unknown method {method} in voronoi()")
 332        raise RuntimeError
 333
 334    m.lw(2).lighting("off").wireframe()
 335    m.name = "Voronoi"
 336    return m
 337
 338
 339def _rotate_points(points, n0=None, n1=(0, 0, 1)):
 340    # Rotate a set of 3D points from direction n0 to direction n1.
 341    # Return the rotated points and the normal to the fitting plane (if n0 is None).
 342    # The pointing direction of the normal in this case is arbitrary.
 343    points = np.asarray(points)
 344
 345    if points.ndim == 1:
 346        points = points[np.newaxis, :]
 347
 348    if len(points[0]) == 2:
 349        return points, (0, 0, 1)
 350
 351    if n0 is None:  # fit plane
 352        datamean = points.mean(axis=0)
 353        vv = np.linalg.svd(points - datamean)[2]
 354        n0 = np.cross(vv[0], vv[1])
 355
 356    n0 = n0 / np.linalg.norm(n0)
 357    n1 = n1 / np.linalg.norm(n1)
 358    k = np.cross(n0, n1)
 359    l = np.linalg.norm(k)
 360    if not l:
 361        k = n0
 362    k /= np.linalg.norm(k)
 363
 364    ct = np.dot(n0, n1)
 365    theta = np.arccos(ct)
 366    st = np.sin(theta)
 367    v = k * (1 - ct)
 368
 369    rpoints = []
 370    for p in points:
 371        a = p * ct
 372        b = np.cross(k, p) * st
 373        c = v * np.dot(k, p)
 374        rpoints.append(a + b + c)
 375
 376    return np.array(rpoints), n0
 377
 378
 379def fit_line(points):
 380    """
 381    Fits a line through points.
 382
 383    Extra info is stored in `Line.slope`, `Line.center`, `Line.variances`.
 384
 385    Examples:
 386        - [fitline.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/fitline.py)
 387
 388            ![](https://vedo.embl.es/images/advanced/fitline.png)
 389    """
 390    if isinstance(points, Points):
 391        points = points.points()
 392    data = np.array(points)
 393    datamean = data.mean(axis=0)
 394    _, dd, vv = np.linalg.svd(data - datamean)
 395    vv = vv[0] / np.linalg.norm(vv[0])
 396    # vv contains the first principal component, i.e. the direction
 397    # vector of the best fit line in the least squares sense.
 398    xyz_min = points.min(axis=0)
 399    xyz_max = points.max(axis=0)
 400    a = np.linalg.norm(xyz_min - datamean)
 401    b = np.linalg.norm(xyz_max - datamean)
 402    p1 = datamean - a * vv
 403    p2 = datamean + b * vv
 404    line = vedo.shapes.Line(p1, p2, lw=1)
 405    line.slope = vv
 406    line.center = datamean
 407    line.variances = dd
 408    return line
 409
 410
 411def fit_circle(points):
 412    """
 413    Fits a circle through a set of 3D points, with a very fast non-iterative method.
 414
 415    Returns the tuple `(center, radius, normal_to_circle)`.
 416
 417    .. warning::
 418        trying to fit s-shaped points will inevitably lead to instabilities and
 419        circles of small radius.
 420
 421    References:
 422        *J.F. Crawford, Nucl. Instr. Meth. 211, 1983, 223-225.*
 423    """
 424    data = np.asarray(points)
 425
 426    offs = data.mean(axis=0)
 427    data, n0 = _rotate_points(data - offs)
 428
 429    xi = data[:, 0]
 430    yi = data[:, 1]
 431
 432    x = sum(xi)
 433    xi2 = xi * xi
 434    xx = sum(xi2)
 435    xxx = sum(xi2 * xi)
 436
 437    y = sum(yi)
 438    yi2 = yi * yi
 439    yy = sum(yi2)
 440    yyy = sum(yi2 * yi)
 441
 442    xiyi = xi * yi
 443    xy = sum(xiyi)
 444    xyy = sum(xiyi * yi)
 445    xxy = sum(xi * xiyi)
 446
 447    N = len(xi)
 448    k = (xx + yy) / N
 449
 450    a1 = xx - x * x / N
 451    b1 = xy - x * y / N
 452    c1 = 0.5 * (xxx + xyy - x * k)
 453
 454    a2 = xy - x * y / N
 455    b2 = yy - y * y / N
 456    c2 = 0.5 * (xxy + yyy - y * k)
 457
 458    d = a2 * b1 - a1 * b2
 459    if not d:
 460        return offs, 0, n0
 461    x0 = (b1 * c2 - b2 * c1) / d
 462    y0 = (c1 - a1 * x0) / b1
 463
 464    R = np.sqrt(x0 * x0 + y0 * y0 - 1 / N * (2 * x0 * x + 2 * y0 * y - xx - yy))
 465
 466    c, _ = _rotate_points([x0, y0, 0], (0, 0, 1), n0)
 467
 468    return c[0] + offs, R, n0
 469
 470
 471def fit_plane(points, signed=False):
 472    """
 473    Fits a plane to a set of points.
 474
 475    Extra info is stored in `Plane.normal`, `Plane.center`, `Plane.variance`.
 476
 477    Arguments:
 478    signed : (bool)
 479        if True flip sign of the normal based on the ordering of the points
 480
 481    Examples:
 482        - [fitline.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/fitline.py)
 483
 484            ![](https://vedo.embl.es/images/advanced/fitline.png)
 485    """
 486    if isinstance(points, Points):
 487        points = points.points()
 488    data = np.asarray(points)
 489    datamean = data.mean(axis=0)
 490    pts = data - datamean
 491    res = np.linalg.svd(pts)
 492    dd, vv = res[1], res[2]
 493    n = np.cross(vv[0], vv[1])
 494    if signed:
 495        v = np.zeros_like(pts)
 496        for i in range(len(pts) - 1):
 497            vi = np.cross(pts[i], pts[i + 1])
 498            v[i] = vi / np.linalg.norm(vi)
 499        ns = np.mean(v, axis=0)  # normal to the points plane
 500        if np.dot(n, ns) < 0:
 501            n = -n
 502    xyz_min = data.min(axis=0)
 503    xyz_max = data.max(axis=0)
 504    s = np.linalg.norm(xyz_max - xyz_min)
 505    pla = vedo.shapes.Plane(datamean, n, s=[s, s])
 506    pla.normal = n
 507    pla.center = datamean
 508    pla.variance = dd[2]
 509    pla.name = "FitPlane"
 510    pla.top = n
 511    return pla
 512
 513
 514def fit_sphere(coords):
 515    """
 516    Fits a sphere to a set of points.
 517
 518    Extra info is stored in `Sphere.radius`, `Sphere.center`, `Sphere.residue`.
 519
 520    Examples:
 521        - [fitspheres1.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/fitspheres1.py)
 522
 523            ![](https://vedo.embl.es/images/advanced/fitspheres1.jpg)
 524    """
 525    if isinstance(coords, Points):
 526        coords = coords.points()
 527    coords = np.array(coords)
 528    n = len(coords)
 529    A = np.zeros((n, 4))
 530    A[:, :-1] = coords * 2
 531    A[:, 3] = 1
 532    f = np.zeros((n, 1))
 533    x = coords[:, 0]
 534    y = coords[:, 1]
 535    z = coords[:, 2]
 536    f[:, 0] = x * x + y * y + z * z
 537    try:
 538        C, residue, rank, _ = np.linalg.lstsq(A, f, rcond=-1)  # solve AC=f
 539    except:
 540        C, residue, rank, _ = np.linalg.lstsq(A, f)  # solve AC=f
 541    if rank < 4:
 542        return None
 543    t = (C[0] * C[0]) + (C[1] * C[1]) + (C[2] * C[2]) + C[3]
 544    radius = np.sqrt(t)[0]
 545    center = np.array([C[0][0], C[1][0], C[2][0]])
 546    if len(residue):
 547        residue = np.sqrt(residue[0]) / n
 548    else:
 549        residue = 0
 550    sph = vedo.shapes.Sphere(center, radius, c=(1, 0, 0)).wireframe(1)
 551    sph.radius = radius
 552    sph.center = center
 553    sph.residue = residue
 554    sph.name = "FitSphere"
 555    return sph
 556
 557
 558def pca_ellipse(points, pvalue=0.673, res=60):
 559    """
 560    Show the oriented PCA 2D ellipse that contains the fraction `pvalue` of points.
 561
 562    Parameter `pvalue` sets the specified fraction of points inside the ellipse.
 563    Normalized directions are stored in `ellipse.axis1`, `ellipse.axis12`
 564    axes sizes are stored in `ellipse.va`, `ellipse.vb`
 565
 566    Arguments:
 567        pvalue : (float)
 568            ellipse will include this fraction of points
 569        res : (int)
 570            resolution of the ellipse
 571
 572    Examples:
 573        - [histo_pca.py](https://github.com/marcomusy/vedo/tree/master/examples/pyplot/histo_pca.py)
 574
 575            ![](https://vedo.embl.es/images/pyplot/histo_pca.png)
 576    """
 577    from scipy.stats import f
 578
 579    if isinstance(points, Points):
 580        coords = points.points()
 581    else:
 582        coords = points
 583    if len(coords) < 4:
 584        vedo.logger.warning("in pca_ellipse(), there are not enough points!")
 585        return None
 586
 587    P = np.array(coords, dtype=float)[:,(0,1)]
 588    cov = np.cov(P, rowvar=0)     # covariance matrix
 589    _, s, R = np.linalg.svd(cov)  # singular value decomposition
 590    p, n = s.size, P.shape[0]
 591    fppf = f.ppf(pvalue, p, n-p)  # f % point function
 592    ua, ub = np.sqrt(s*fppf/2)*2  # semi-axes (largest first)
 593    center = np.mean(P, axis=0)   # centroid of the ellipse
 594
 595    matri = vtk.vtkMatrix4x4()
 596    matri.DeepCopy((
 597        R[0][0] * ua, R[1][0] * ub, 0, center[0],
 598        R[0][1] * ua, R[1][1] * ub, 0, center[1],
 599                   0,            0, 1,         0,
 600        0, 0, 0, 1)
 601    )
 602    vtra = vtk.vtkTransform()
 603    vtra.SetMatrix(matri)
 604
 605    elli = vedo.shapes.Circle(alpha=0.75, res=res)
 606
 607    # assign the transformation
 608    elli.SetScale(vtra.GetScale())
 609    elli.SetOrientation(vtra.GetOrientation())
 610    elli.SetPosition(vtra.GetPosition())
 611
 612    elli.center = np.array(vtra.GetPosition())
 613    elli.nr_of_points = n
 614    elli.va = ua
 615    elli.vb = ub
 616    elli.axis1 = np.array(vtra.TransformPoint([1, 0, 0])) - elli.center
 617    elli.axis2 = np.array(vtra.TransformPoint([0, 1, 0])) - elli.center
 618    elli.axis1 /= np.linalg.norm(elli.axis1)
 619    elli.axis2 /= np.linalg.norm(elli.axis2)
 620    elli.transformation = vtra
 621    elli.name = "PCAEllipse"
 622    return elli
 623
 624@deprecated(reason=vedo.colors.red + "Please use pca_ellipsoid()" + vedo.colors.reset)
 625def pcaEllipsoid(points, pvalue=0.673):
 626    "Deprecated. Please use `pca_ellipsoid()`."
 627    return pca_ellipsoid(points, pvalue)
 628
 629def pca_ellipsoid(points, pvalue=0.673):
 630    """
 631    Show the oriented PCA ellipsoid that contains fraction `pvalue` of points.
 632
 633    Parameter `pvalue` sets the specified fraction of points inside the ellipsoid.
 634
 635    Extra can be calculated with `mesh.asphericity()`, `mesh.asphericity_error()`
 636    (asphericity is equal to 0 for a perfect sphere).
 637
 638    Axes sizes can be accessed in `ellips.va`, `ellips.vb`, `ellips.vc`,
 639    normalized directions are stored in `ellips.axis1`, `ellips.axis12`
 640    and `ellips.axis3`.
 641
 642    .. warning:: the meaning of `ellips.axis1`, has changed wrt `vedo==2022.1.0`
 643        in that it's now the direction wrt the origin (e.i. the center is subtracted)
 644
 645    Examples:
 646        - [pca_ellipsoid.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/pca_ellipsoid.py)
 647
 648            ![](https://vedo.embl.es/images/basic/pca.png)
 649    """
 650    from scipy.stats import f
 651
 652    if isinstance(points, Points):
 653        coords = points.points()
 654    else:
 655        coords = points
 656    if len(coords) < 4:
 657        vedo.logger.warning("in pcaEllipsoid(), there are not enough points!")
 658        return None
 659
 660    P = np.array(coords, ndmin=2, dtype=float)
 661    cov = np.cov(P, rowvar=0)     # covariance matrix
 662    _, s, R = np.linalg.svd(cov)  # singular value decomposition
 663    p, n = s.size, P.shape[0]
 664    fppf = f.ppf(pvalue, p, n-p)*(n-1)*p*(n+1)/n/(n-p)  # f % point function
 665    cfac = 1 + 6/(n-1)            # correction factor for low statistics
 666    ua, ub, uc = np.sqrt(s*fppf)/cfac  # semi-axes (largest first)
 667    center = np.mean(P, axis=0)   # centroid of the hyperellipsoid
 668
 669    elli = vedo.shapes.Ellipsoid((0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1), alpha=0.25)
 670
 671    matri = vtk.vtkMatrix4x4()
 672    matri.DeepCopy((
 673        R[0][0] * ua*2, R[1][0] * ub*2, R[2][0] * uc*2, center[0],
 674        R[0][1] * ua*2, R[1][1] * ub*2, R[2][1] * uc*2, center[1],
 675        R[0][2] * ua*2, R[1][2] * ub*2, R[2][2] * uc*2, center[2],
 676        0, 0, 0, 1)
 677    )
 678    vtra = vtk.vtkTransform()
 679    vtra.SetMatrix(matri)
 680
 681    # assign the transformation
 682    elli.SetScale(vtra.GetScale())
 683    elli.SetOrientation(vtra.GetOrientation())
 684    elli.SetPosition(vtra.GetPosition())
 685
 686    elli.center = np.array(vtra.GetPosition())
 687    elli.nr_of_points = n
 688    elli.va = ua
 689    elli.vb = ub
 690    elli.vc = uc
 691    elli.axis1 = np.array(vtra.TransformPoint([1, 0, 0])) - elli.center
 692    elli.axis2 = np.array(vtra.TransformPoint([0, 1, 0])) - elli.center
 693    elli.axis3 = np.array(vtra.TransformPoint([0, 0, 1])) - elli.center
 694    elli.axis1 /= np.linalg.norm(elli.axis1)
 695    elli.axis2 /= np.linalg.norm(elli.axis2)
 696    elli.axis3 /= np.linalg.norm(elli.axis3)
 697    elli.transformation = vtra
 698    elli.name = "PCAEllipsoid"
 699    return elli
 700
 701
 702###################################################
 703def Point(pos=(0, 0, 0), r=12, c="red", alpha=1):
 704    """
 705    Create a simple point in space.
 706
 707    .. note:: if you are creating many points you should definitely use class `Points` instead!
 708    """
 709    if isinstance(pos, vtk.vtkActor):
 710        pos = pos.GetPosition()
 711    pd = utils.buildPolyData([[0, 0, 0]])
 712    if len(pos) == 2:
 713        pos = (pos[0], pos[1], 0.0)
 714    pt = Points(pd, c, alpha, r)
 715    pt.SetPosition(pos)
 716    pt.name = "Point"
 717    return pt
 718
 719
 720###################################################
 721class Points(BaseActor, vtk.vtkActor):
 722    """Work with pointclouds."""
 723
 724    def __init__(
 725        self, inputobj=None, c=(0.2, 0.2, 0.2), alpha=1, r=4,
 726    ):
 727        """
 728        Build an object made of only vertex points for a list of 2D/3D points.
 729        Both shapes (N, 3) or (3, N) are accepted as input, if N>3.
 730        For very large point clouds a list of colors and alpha can be assigned to each
 731        point in the form c=[(R,G,B,A), ... ] where 0<=R<256, ... 0<=A<256.
 732
 733        Arguments:
 734            inputobj : (list, tuple)
 735            c : (str, list)
 736                Color name or rgb tuple.
 737            alpha : (float)
 738                Transparency in range [0,1].
 739            r : (int)
 740                Point radius in units of pixels.
 741
 742        Example:
 743            ```python
 744            from vedo import *
 745
 746            def fibonacci_sphere(n):
 747                s = np.linspace(0, n, num=n, endpoint=False)
 748                theta = s * 2.399963229728653
 749                y = 1 - s * (2/(n-1))
 750                r = np.sqrt(1 - y * y)
 751                x = np.cos(theta) * r
 752                z = np.sin(theta) * r
 753                return [x,y,z]
 754
 755            Points(fibonacci_sphere(1000)).show(axes=1).close()
 756            ```
 757            ![](https://vedo.embl.es/images/feats/fibonacci.png)
 758        """
 759
 760        vtk.vtkActor.__init__(self)
 761        BaseActor.__init__(self)
 762
 763        self._data = None
 764
 765        self._mapper = vtk.vtkPolyDataMapper()
 766        self.SetMapper(self._mapper)
 767
 768        self._bfprop = None  # backface property holder
 769
 770        self._scals_idx = 0  # index of the active scalar changed from CLI
 771        self._ligthingnr = 0  # index of the lighting mode changed from CLI
 772        self._cmap_name = ""  # remember the name for self._keypress
 773        # self.name = "Points" # better not to give it a name here
 774
 775        self.property = self.GetProperty()
 776        try:
 777            self.property.RenderPointsAsSpheresOn()
 778        except AttributeError:
 779            pass
 780
 781        if inputobj is None:  ####################
 782            self._data = vtk.vtkPolyData()
 783            return
 784        ########################################
 785
 786        self.property.SetRepresentationToPoints()
 787        self.property.SetPointSize(r)
 788        self.property.LightingOff()
 789
 790        if isinstance(inputobj, vedo.BaseActor):
 791            inputobj = inputobj.points()  # numpy
 792
 793        ######
 794        if isinstance(inputobj, vtk.vtkActor):
 795            poly_copy = vtk.vtkPolyData()
 796            pr = vtk.vtkProperty()
 797            pr.DeepCopy(inputobj.GetProperty())
 798            poly_copy.DeepCopy(inputobj.GetMapper().GetInput())
 799            pr.SetRepresentationToPoints()
 800            pr.SetPointSize(r)
 801            self._data = poly_copy
 802            self._mapper.SetInputData(poly_copy)
 803            self._mapper.SetScalarVisibility(inputobj.GetMapper().GetScalarVisibility())
 804            self.SetProperty(pr)
 805            self.property = pr
 806
 807        elif isinstance(inputobj, vtk.vtkPolyData):
 808            if inputobj.GetNumberOfCells() == 0:
 809                carr = vtk.vtkCellArray()
 810                for i in range(inputobj.GetNumberOfPoints()):
 811                    carr.InsertNextCell(1)
 812                    carr.InsertCellPoint(i)
 813                inputobj.SetVerts(carr)
 814            self._data = inputobj  # cache vtkPolyData and mapper for speed
 815
 816        elif utils.is_sequence(inputobj):  # passing point coords
 817            plist = inputobj
 818            n = len(plist)
 819
 820            if n == 3:  # assume plist is in the format [all_x, all_y, all_z]
 821                if utils.is_sequence(plist[0]) and len(plist[0]) > 3:
 822                    plist = np.stack((plist[0], plist[1], plist[2]), axis=1)
 823            elif n == 2:  # assume plist is in the format [all_x, all_y, 0]
 824                if utils.is_sequence(plist[0]) and len(plist[0]) > 3:
 825                    plist = np.stack((plist[0], plist[1], np.zeros(len(plist[0]))), axis=1)
 826
 827            # if n and len(plist[0]) == 2:  # make it 3d
 828            #     plist = np.c_[np.array(plist), np.zeros(len(plist))]
 829            plist = utils.make3d(plist)
 830
 831            if (
 832                utils.is_sequence(c)
 833                and (len(c) > 3 or (utils.is_sequence(c[0]) and len(c[0]) == 4))
 834            ) or utils.is_sequence(alpha):
 835
 836                cols = c
 837
 838                n = len(plist)
 839                if n != len(cols):
 840                    vedo.logger.error(f"mismatch in Points() colors array lengths {n} and {len(cols)}")
 841                    raise RuntimeError()
 842
 843                src = vtk.vtkPointSource()
 844                src.SetNumberOfPoints(n)
 845                src.Update()
 846
 847                vgf = vtk.vtkVertexGlyphFilter()
 848                vgf.SetInputData(src.GetOutput())
 849                vgf.Update()
 850                pd = vgf.GetOutput()
 851
 852                pd.GetPoints().SetData(utils.numpy2vtk(plist, dtype=np.float32))
 853
 854                ucols = vtk.vtkUnsignedCharArray()
 855                ucols.SetNumberOfComponents(4)
 856                ucols.SetName("Points_RGBA")
 857                if utils.is_sequence(alpha):
 858                    if len(alpha) != n:
 859                        vedo.logger.error(f"mismatch in Points() alpha array lengths {n} and {len(cols)}")
 860                        raise RuntimeError()
 861                    alphas = alpha
 862                    alpha = 1
 863                else:
 864                    alphas = (alpha,) * n
 865
 866                if utils.is_sequence(cols):
 867                    c = None
 868                    if len(cols[0]) == 4:
 869                        for i in range(n):  # FAST
 870                            rc, gc, bc, ac = cols[i]
 871                            ucols.InsertNextTuple4(rc, gc, bc, ac)
 872                    else:
 873                        for i in range(n): # SLOW
 874                            rc,gc,bc = colors.get_color(cols[i])
 875                            ucols.InsertNextTuple4(rc*255, gc*255, bc*255, alphas[i]*255)
 876                else:
 877                    c = cols
 878
 879                pd.GetPointData().AddArray(ucols)
 880                pd.GetPointData().SetActiveScalars("Points_RGBA")
 881                self._mapper.SetInputData(pd)
 882                self._mapper.ScalarVisibilityOn()
 883                self._data = pd
 884
 885            else:
 886
 887                pd = utils.buildPolyData(plist)
 888                self._mapper.SetInputData(pd)
 889                c = colors.get_color(c)
 890                self.property.SetColor(c)
 891                self.property.SetOpacity(alpha)
 892                self._data = pd
 893
 894            ##########
 895            self.pipeline = utils.OperationNode(
 896                self, parents=[], comment=f"#pts {self._data.GetNumberOfPoints()}",
 897            )
 898            return
 899            ##########
 900
 901        elif isinstance(inputobj, str):
 902            verts = vedo.io.load(inputobj)
 903            self.filename = inputobj
 904            self._data = verts.polydata()
 905
 906        else:
 907
 908            # try to extract the points from the VTK input data object
 909            try:
 910                vvpts = inputobj.GetPoints()
 911                pd = vtk.vtkPolyData()
 912                pd.SetPoints(vvpts)
 913                for i in range(inputobj.GetPointData().GetNumberOfArrays()):
 914                    arr = inputobj.GetPointData().GetArray(i)
 915                    pd.GetPointData().AddArray(arr)
 916
 917                self._mapper.SetInputData(pd)
 918                c = colors.get_color(c)
 919                self.property.SetColor(c)
 920                self.property.SetOpacity(alpha)
 921                self._data = pd
 922            except:
 923                vedo.logger.error(f"cannot build Points from type {type(inputobj)}")
 924                raise RuntimeError()
 925
 926        c = colors.get_color(c)
 927        self.property.SetColor(c)
 928        self.property.SetOpacity(alpha)
 929
 930        self._mapper.SetInputData(self._data)
 931
 932        self.pipeline = utils.OperationNode(
 933            self, parents=[], comment=f"#pts {self._data.GetNumberOfPoints()}",
 934        )
 935        return
 936
 937    def _repr_html_(self):
 938        """
 939        HTML representation of the Point cloud object for Jupyter Notebooks.
 940        
 941        Returns:
 942            HTML text with the image and some properties.
 943        """
 944        import io
 945        import base64
 946        from PIL import Image
 947
 948        library_name = "vedo.pointcloud.Points"
 949        help_url = "https://vedo.embl.es/docs/vedo/pointcloud.html"
 950
 951        arr = self.thumbnail()
 952        im = Image.fromarray(arr)
 953        buffered = io.BytesIO()
 954        im.save(buffered, format="PNG", quality=100)
 955        encoded = base64.b64encode(buffered.getvalue()).decode("utf-8")
 956        url = "data:image/png;base64," + encoded
 957        image = f"<img src='{url}'></img>"
 958
 959        # mesh statisitics
 960        bounds = "<br/>".join(
 961            [
 962                utils.precision(min_x,4) + " ... " + utils.precision(max_x,4)
 963                for min_x, max_x in zip(self.bounds()[::2], self.bounds()[1::2])
 964            ]
 965        )
 966        average_size = "{size:.3f}".format(size=self.average_size())
 967
 968        help_text = ""
 969        if self.name:
 970            help_text += f"<b> {self.name}: &nbsp&nbsp</b>"
 971        help_text += '<b><a href="' + help_url + '" target="_blank">' + library_name + "</a></b>" 
 972        if self.filename:
 973            dots = ""
 974            if len(self.filename) > 30:
 975                dots = "..."
 976            help_text += f"<br/><code><i>({dots}{self.filename[-30:]})</i></code>"
 977        
 978        pdata = ""
 979        if self._data.GetPointData().GetScalars():
 980            if self._data.GetPointData().GetScalars().GetName():
 981                name = self._data.GetPointData().GetScalars().GetName()
 982                pdata = "<tr><td><b> point data array </b></td><td>" + name + "</td></tr>"
 983
 984        cdata = ""
 985        if self._data.GetCellData().GetScalars():
 986            if self._data.GetCellData().GetScalars().GetName():
 987                name = self._data.GetCellData().GetScalars().GetName()
 988                cdata = "<tr><td><b> cell data array </b></td><td>" + name + "</td></tr>"
 989
 990        all = [
 991            "<table>",
 992            "<tr>", 
 993            "<td>", image, "</td>",
 994            "<td style='text-align: center; vertical-align: center;'><br/>", help_text,
 995            "<table>",
 996            "<tr><td><b> bounds </b> <br/> (x/y/z) </td><td>" + str(bounds) + "</td></tr>",
 997            "<tr><td><b> center of mass </b></td><td>" + utils.precision(self.center_of_mass(),3) + "</td></tr>",
 998            "<tr><td><b> average size </b></td><td>" + str(average_size) + "</td></tr>",
 999            "<tr><td><b> nr. points </b></td><td>" + str(self.npoints) + "</td></tr>",
1000            pdata,
1001            cdata,
1002            "</table>",
1003            "</table>",
1004        ]
1005        return "\n".join(all)
1006
1007
1008    ##################################################################################
1009    def _update(self, polydata):
1010        # Overwrite the polygonal mesh with a new vtkPolyData
1011        self._data = polydata
1012        self.mapper().SetInputData(polydata)
1013        self.mapper().Modified()
1014        return self
1015
1016    def __add__(self, meshs):
1017        if isinstance(meshs, list):
1018            alist = [self]
1019            for l in meshs:
1020                if isinstance(l, vedo.Assembly):
1021                    alist += l.unpack()
1022                else:
1023                    alist += l
1024            return vedo.assembly.Assembly(alist)
1025
1026        if isinstance(meshs, vedo.Assembly):
1027            return meshs + self  # use Assembly.__add__
1028
1029        return vedo.assembly.Assembly([self, meshs])
1030
1031    def polydata(self, transformed=True):
1032        """
1033        Returns the `vtkPolyData` object associated to a `Mesh`.
1034
1035        .. note::
1036            If `transformed=True` return a copy of polydata that corresponds
1037            to the current mesh position in space.
1038        """
1039        if not self._data:
1040            self._data = self.mapper().GetInput()
1041            return self._data
1042
1043        if transformed:
1044            # if self.GetIsIdentity() or self._data.GetNumberOfPoints()==0: # commmentd out on 15th feb 2020
1045            if self._data.GetNumberOfPoints() == 0:
1046                # no need to do much
1047                return self._data
1048
1049            # otherwise make a copy that corresponds to
1050            # the actual position in space of the mesh
1051            M = self.GetMatrix()
1052            transform = vtk.vtkTransform()
1053            transform.SetMatrix(M)
1054            tp = vtk.vtkTransformPolyDataFilter()
1055            tp.SetTransform(transform)
1056            tp.SetInputData(self._data)
1057            tp.Update()
1058            return tp.GetOutput()
1059
1060        return self._data
1061
1062
1063    def vertices(self, pts=None, transformed=True):
1064        """Alias for `points()`."""
1065        return self.points(pts, transformed)
1066
1067    def clone(self, deep=True, transformed=False):
1068        """
1069        Clone a `PointCloud` or `Mesh` object to make an exact copy of it.
1070
1071        Arguments:
1072            deep : (bool)
1073                if False only build a shallow copy of the object (faster copy).
1074
1075            transformed : (bool)
1076                if True reset the current transformation of the copy to unit.
1077
1078        Examples:
1079            - [mirror.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/mirror.py)
1080
1081               ![](https://vedo.embl.es/images/basic/mirror.png)
1082        """
1083        poly = self.polydata(transformed)
1084        poly_copy = vtk.vtkPolyData()
1085        if deep:
1086            poly_copy.DeepCopy(poly)
1087        else:
1088            poly_copy.ShallowCopy(poly)
1089
1090        if isinstance(self, vedo.Mesh):
1091            cloned = vedo.Mesh(poly_copy)
1092        else:
1093            cloned = Points(poly_copy)
1094
1095        pr = vtk.vtkProperty()
1096        pr.DeepCopy(self.GetProperty())
1097        cloned.SetProperty(pr)
1098        cloned.property = pr
1099
1100        if self.GetBackfaceProperty():
1101            bfpr = vtk.vtkProperty()
1102            bfpr.DeepCopy(self.GetBackfaceProperty())
1103            cloned.SetBackfaceProperty(bfpr)
1104
1105        if not transformed:
1106            if self.transform:
1107                # already has a transform which can be non linear, so use that
1108                cloned.SetUserTransform(self.transform)
1109            else:
1110                # assign the same transformation to the copy
1111                cloned.SetOrigin(self.GetOrigin())
1112                cloned.SetScale(self.GetScale())
1113                cloned.SetOrientation(self.GetOrientation())
1114                cloned.SetPosition(self.GetPosition())
1115
1116        mp = cloned.mapper()
1117        sm = self.mapper()
1118        mp.SetScalarVisibility(sm.GetScalarVisibility())
1119        mp.SetScalarRange(sm.GetScalarRange())
1120        mp.SetColorMode(sm.GetColorMode())
1121        lsr = sm.GetUseLookupTableScalarRange()
1122        mp.SetUseLookupTableScalarRange(lsr)
1123        mp.SetScalarMode(sm.GetScalarMode())
1124        lut = sm.GetLookupTable()
1125        if lut:
1126            mp.SetLookupTable(lut)
1127
1128        cloned.SetPickable(self.GetPickable())
1129
1130        cloned.base = np.array(self.base)
1131        cloned.top = np.array(self.top)
1132        cloned.name = str(self.name)
1133        cloned.filename = str(self.filename)
1134        cloned.info = dict(self.info)
1135
1136        # better not to share the same locators with original obj
1137        cloned.point_locator = None
1138        cloned.cell_locator = None
1139
1140        cloned.pipeline = utils.OperationNode(
1141            "clone", parents=[self], shape='diamond', c='#edede9',
1142        )
1143        return cloned
1144
1145    def clone2d(
1146        self,
1147        pos=(0, 0),
1148        coordsys=4,
1149        scale=None,
1150        c=None,
1151        alpha=None,
1152        ps=2,
1153        lw=1,
1154        sendback=False,
1155        layer=0,
1156    ):
1157        """
1158        Copy a 3D Mesh into a static 2D image. Returns a `vtkActor2D`.
1159
1160        Arguments:
1161            coordsys : (int)
1162                the coordinate system, options are
1163                - 0 = Displays
1164                - 1 = Normalized Display
1165                - 2 = Viewport (origin is the bottom-left corner of the window)
1166                - 3 = Normalized Viewport
1167                - 4 = View (origin is the center of the window)
1168                - 5 = World (anchor the 2d image to mesh)
1169
1170            ps : (int)
1171                point size in pixel units
1172
1173            lw : (int)
1174                line width in pixel units
1175
1176            sendback : (bool)
1177                put it behind any other 3D object
1178
1179        Examples:
1180            - [clone2d.py](https://github.com/marcomusy/vedo/tree/master/examples/other/clone2d.py)
1181
1182                ![](https://vedo.embl.es/images/other/clone2d.png)
1183        """
1184        if scale is None:
1185            msiz = self.diagonal_size()
1186            if vedo.plotter_instance and vedo.plotter_instance.window:
1187                sz = vedo.plotter_instance.window.GetSize()
1188                dsiz = utils.mag(sz)
1189                scale = dsiz / msiz / 10
1190            else:
1191                scale = 350 / msiz
1192
1193        cmsh = self.clone()
1194        poly = cmsh.pos(0, 0, 0).scale(scale).polydata()
1195
1196        mapper3d = self.mapper()
1197        cm = mapper3d.GetColorMode()
1198        lut = mapper3d.GetLookupTable()
1199        sv = mapper3d.GetScalarVisibility()
1200        use_lut = mapper3d.GetUseLookupTableScalarRange()
1201        vrange = mapper3d.GetScalarRange()
1202        sm = mapper3d.GetScalarMode()
1203
1204        mapper2d = vtk.vtkPolyDataMapper2D()
1205        mapper2d.ShallowCopy(mapper3d)
1206        mapper2d.SetInputData(poly)
1207        mapper2d.SetColorMode(cm)
1208        mapper2d.SetLookupTable(lut)
1209        mapper2d.SetScalarVisibility(sv)
1210        mapper2d.SetUseLookupTableScalarRange(use_lut)
1211        mapper2d.SetScalarRange(vrange)
1212        mapper2d.SetScalarMode(sm)
1213
1214        act2d = vtk.vtkActor2D()
1215        act2d.SetMapper(mapper2d)
1216        act2d.SetLayerNumber(layer)
1217        csys = act2d.GetPositionCoordinate()
1218        csys.SetCoordinateSystem(coordsys)
1219        act2d.SetPosition(pos)
1220        if c is not None:
1221            c = colors.get_color(c)
1222            act2d.GetProperty().SetColor(c)
1223            mapper2d.SetScalarVisibility(False)
1224        else:
1225            act2d.GetProperty().SetColor(cmsh.color())
1226        if alpha is not None:
1227            act2d.GetProperty().SetOpacity(alpha)
1228        else:
1229            act2d.GetProperty().SetOpacity(cmsh.alpha())
1230        act2d.GetProperty().SetPointSize(ps)
1231        act2d.GetProperty().SetLineWidth(lw)
1232        act2d.GetProperty().SetDisplayLocationToForeground()
1233        if sendback:
1234            act2d.GetProperty().SetDisplayLocationToBackground()
1235
1236        # print(csys.GetCoordinateSystemAsString())
1237        # print(act2d.GetHeight(), act2d.GetWidth(), act2d.GetLayerNumber())
1238        return act2d
1239
1240    def add_trail(self, offset=(0, 0, 0), n=50, c=None, alpha=1, lw=2):
1241        """
1242        Add a trailing line to mesh.
1243        This new mesh is accessible through `mesh.trail`.
1244
1245        Arguments:
1246            offset : (float)
1247                set an offset vector from the object center.
1248            n : (int)
1249                number of segments
1250            lw : (float)
1251                line width of the trail
1252
1253        Examples:
1254            - [trail.py](https://github.com/marcomusy/vedo/tree/master/examples/simulations/trail.py)
1255
1256                ![](https://vedo.embl.es/images/simulations/trail.gif)
1257
1258            - [airplane1.py](https://github.com/marcomusy/vedo/tree/master/examples/simulations/airplane1.py)
1259            - [airplane2.py](https://github.com/marcomusy/vedo/tree/master/examples/simulations/airplane2.py)
1260        """
1261        if self.trail is None:
1262            pos = self.GetPosition()
1263            self.trail_offset = np.asarray(offset)
1264            self.trail_points = [pos] * n
1265
1266            if c is None:
1267                col = self.GetProperty().GetColor()
1268            else:
1269                col = colors.get_color(c)
1270
1271            tline = vedo.shapes.Line(pos, pos, res=n, c=col, alpha=alpha, lw=lw)
1272            self.trail = tline  # holds the Line
1273        return self
1274
1275    def _update_trail(self):
1276        if isinstance(self, vedo.shapes.Arrow):
1277            currentpos = self.tipPoint()  # the tip of Arrow
1278        else:
1279            currentpos = np.array(self.GetPosition())
1280
1281        self.trail_points.append(currentpos)  # cycle
1282        self.trail_points.pop(0)
1283
1284        data = np.array(self.trail_points) - currentpos + self.trail_offset
1285        tpoly = self.trail.polydata(False)
1286        tpoly.GetPoints().SetData(utils.numpy2vtk(data, dtype=np.float32))
1287        self.trail.SetPosition(currentpos)
1288
1289    def _update_shadows(self):
1290        shadows = self.shadows
1291        self.shadows = []
1292        for sha in shadows:
1293            color = sha.GetProperty().GetColor()
1294            opacity = sha.GetProperty().GetOpacity()
1295            self.add_shadow(**sha.info, c=color, alpha=opacity)
1296
1297    def delete_cells_by_point_index(self, indices):
1298        """
1299        Delete a list of vertices identified by any of their vertex index.
1300
1301        See also `delete_cells()`.
1302
1303        Examples:
1304            - [elete_mesh_pts.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/elete_mesh_pts.py)
1305
1306                ![](https://vedo.embl.es/images/basic/deleteMeshPoints.png)
1307        """
1308        cell_ids = vtk.vtkIdList()
1309        data = self.inputdata()
1310        data.BuildLinks()
1311        n = 0
1312        for i in np.unique(indices):
1313            data.GetPointCells(i, cell_ids)
1314            for j in range(cell_ids.GetNumberOfIds()):
1315                data.DeleteCell(cell_ids.GetId(j))  # flag cell
1316                n += 1
1317
1318        data.RemoveDeletedCells()
1319        self.mapper().Modified()
1320        self.pipeline = utils.OperationNode(
1321            f"delete {n} cells\nby point index", parents=[self])
1322        return self
1323
1324    def compute_normals_with_pca(self, n=20, orientation_point=None, invert=False):
1325        """
1326        Generate point normals using PCA (principal component analysis).
1327        Basically this estimates a local tangent plane around each sample point p
1328        by considering a small neighborhood of points around p, and fitting a plane
1329        to the neighborhood (via PCA).
1330
1331        Arguments:
1332            n : (int)
1333                neighborhood size to calculate the normal
1334            orientation_point : (list)
1335                adjust the +/- sign of the normals so that
1336                the normals all point towards a specified point. If None, perform a traversal
1337                of the point cloud and flip neighboring normals so that they are mutually consistent.
1338            invert : (bool)
1339                flip all normals
1340        """
1341        poly = self.polydata()
1342        pcan = vtk.vtkPCANormalEstimation()
1343        pcan.SetInputData(poly)
1344        pcan.SetSampleSize(n)
1345
1346        if orientation_point is not None:
1347            pcan.SetNormalOrientationToPoint()
1348            pcan.SetOrientationPoint(orientation_point)
1349        else:
1350            pcan.SetNormalOrientationToGraphTraversal()
1351
1352        if invert:
1353            pcan.FlipNormalsOn()
1354        pcan.Update()
1355
1356        varr = pcan.GetOutput().GetPointData().GetNormals()
1357        varr.SetName("Normals")
1358        self.inputdata().GetPointData().SetNormals(varr)
1359        self.inputdata().GetPointData().Modified()
1360        return self
1361    
1362    def compute_acoplanarity(self, n=25, radius=None, on='points'):
1363        """
1364        Compute acoplanarity which is a measure of how much a local region of the mesh 
1365        differs from a plane.
1366        The information is stored in a `pointdata` or `celldata` array with name 'Acoplanarity'.
1367        Either `n` (number of neighbour points) or `radius` (radius of local search) can be specified.
1368        If a radius value is given and not enough points fall inside it, then a -1 is stored.
1369
1370        Example:
1371            ```python
1372            from vedo import *
1373            msh = ParametricShape('RandomHills')
1374            msh.compute_acoplanarity(radius=0.1, on='cells')
1375            msh.cmap("coolwarm", on='cells').add_scalarbar()
1376            msh.show(axes=1).close()
1377            ```
1378            ![](https://vedo.embl.es/images/feats/acoplanarity.jpg)
1379        """
1380        acoplanarities = []
1381        if 'point'  in on:
1382            pts = self.points()
1383        elif 'cell' in on:
1384            pts = self.cell_centers()
1385        else:
1386            raise ValueError(f"In compute_acoplanarity() set on to either 'cells' or 'points', not {on}")
1387
1388        for p in utils.progressbar(pts, delay=5, width=15, title=f'{on} acoplanarity'):
1389            if n:
1390                data = self.closest_point(p, n=n)
1391                npts = n
1392            elif radius:
1393                data = self.closest_point(p, radius=radius)
1394                npts = len(data)
1395
1396            try:
1397                center = data.mean(axis=0)
1398                res = np.linalg.svd(data - center)
1399                acoplanarities.append(res[1][2]/npts)
1400            except:
1401                acoplanarities.append(-1.0)
1402
1403        if 'point'  in on:
1404            self.pointdata["Acoplanarity"] = np.array(acoplanarities, dtype=float)
1405        else:
1406            self.celldata["Acoplanarity"] = np.array(acoplanarities, dtype=float)
1407        return self
1408
1409    @deprecated(reason=vedo.colors.red + "Please use distance_to()" + vedo.colors.reset)
1410    def distanceTo(self, pcloud, signed=False, invert=False, name="Distance"):
1411        "Please use `distance_to()`."
1412        return self.distance_to(pcloud, signed, invert, name)
1413
1414    def distance_to(self, pcloud, signed=False, invert=False, name="Distance"):
1415        """
1416        Computes the distance from one point cloud or mesh to another point cloud or mesh.
1417        This new `pointdata` array is saved with default name "Distance".
1418
1419        Keywords `signed` and `invert` are used to compute signed distance,
1420        but the mesh in that case must have polygonal faces (not a simple point cloud),
1421        and normals must also be computed.
1422
1423        Examples:
1424            - [distance2mesh.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/distance2mesh.py)
1425
1426                ![](https://vedo.embl.es/images/basic/distance2mesh.png)
1427        """
1428        if pcloud.inputdata().GetNumberOfPolys():
1429
1430            poly1 = self.polydata()
1431            poly2 = pcloud.polydata()
1432            df = vtk.vtkDistancePolyDataFilter()
1433            df.ComputeSecondDistanceOff()
1434            df.SetInputData(0, poly1)
1435            df.SetInputData(1, poly2)
1436            df.SetSignedDistance(signed)
1437            df.SetNegateDistance(invert)
1438            df.Update()
1439            scals = df.GetOutput().GetPointData().GetScalars()
1440            dists = utils.vtk2numpy(scals)
1441
1442        else:  # has no polygons and vtkDistancePolyDataFilter wants them (dont know why)
1443
1444            if signed:
1445                vedo.logger.warning("distanceTo() called with signed=True but input object has no polygons")
1446
1447            if not pcloud.point_locator:
1448                pcloud.point_locator = vtk.vtkPointLocator()
1449                pcloud.point_locator.SetDataSet(pcloud.polydata())
1450                pcloud.point_locator.BuildLocator()
1451
1452            ids = []
1453            ps1 = self.points()
1454            ps2 = pcloud.points()
1455            for p in ps1:
1456                pid = pcloud.point_locator.FindClosestPoint(p)
1457                ids.append(pid)
1458
1459            deltas = ps2[ids] - ps1
1460            dists = np.linalg.norm(deltas, axis=1).astype(np.float32)
1461            scals = utils.numpy2vtk(dists)
1462
1463        scals.SetName(name)
1464        self.inputdata().GetPointData().AddArray(scals)  # must be self.inputdata() !
1465        self.inputdata().GetPointData().SetActiveScalars(scals.GetName())
1466        rng = scals.GetRange()
1467        self.mapper().SetScalarRange(rng[0], rng[1])
1468        self.mapper().ScalarVisibilityOn()
1469        
1470        self.pipeline = utils.OperationNode(
1471            "distance_to", 
1472            parents=[self, pcloud], 
1473            shape="cylinder",
1474            comment=f"#pts {self._data.GetNumberOfPoints()}",
1475        )
1476        return dists
1477
1478    def alpha(self, opacity=None):
1479        """Set/get mesh's transparency. Same as `mesh.opacity()`."""
1480        if opacity is None:
1481            return self.GetProperty().GetOpacity()
1482
1483        self.GetProperty().SetOpacity(opacity)
1484        bfp = self.GetBackfaceProperty()
1485        if bfp:
1486            if opacity < 1:
1487                self._bfprop = bfp
1488                self.SetBackfaceProperty(None)
1489            else:
1490                self.SetBackfaceProperty(self._bfprop)
1491        return self
1492
1493    def opacity(self, alpha=None):
1494        """Set/get mesh's transparency. Same as `mesh.alpha()`."""
1495        return self.alpha(alpha)
1496
1497    def force_opaque(self, value=True):
1498        """ Force the Mesh, Line or point cloud to be treated as opaque"""
1499        ## force the opaque pass, fixes picking in vtk9
1500        # but causes other bad troubles with lines..
1501        self.SetForceOpaque(value)
1502        return self
1503
1504    def force_translucent(self, value=True):
1505        """ Force the Mesh, Line or point cloud to be treated as translucent"""
1506        self.SetForceTranslucent(value)
1507        return self
1508
1509    def occlusion(self, value=None):
1510        """Occlusion strength in range [0,1]."""
1511        if value is None:
1512            return self.GetProperty().GetOcclusionStrength()
1513        self.GetProperty().SetOcclusionStrength(value)
1514        return self
1515
1516
1517    @deprecated(reason=vedo.colors.red + "Please use point_size()" + vedo.colors.reset)
1518    def pointSize(self, value):
1519        "Deprecated. Please use `point_size()`."
1520        return self.point_size(value)
1521
1522    def point_size(self, value):
1523        """Set/get mesh's point size of vertices. Same as `mesh.ps()`"""
1524        if not value:
1525            self.GetProperty().SetRepresentationToSurface()
1526        else:
1527            self.GetProperty().SetRepresentationToPoints()
1528            self.GetProperty().SetPointSize(value)
1529        return self
1530
1531    def ps(self, pointsize=None):
1532        """Set/get mesh's point size of vertices. Same as `mesh.point_size()`"""
1533        return self.point_size(pointsize)
1534
1535    def render_points_as_spheres(self, value=True):
1536        """Make points look spheric or make them look as squares."""
1537        self.GetProperty().SetRenderPointsAsSpheres(value)
1538        return self
1539
1540    def color(self, c=False, alpha=None):
1541        """
1542        Set/get mesh's color.
1543        If None is passed as input, will use colors from active scalars.
1544        Same as `mesh.c()`.
1545        """
1546        # overrides base.color()
1547        if c is False:
1548            return np.array(self.GetProperty().GetColor())
1549        if c is None:
1550            self.mapper().ScalarVisibilityOn()
1551            return self
1552        self.mapper().ScalarVisibilityOff()
1553        cc = colors.get_color(c)
1554        self.GetProperty().SetColor(cc)
1555        if self.trail:
1556            self.trail.GetProperty().SetColor(cc)
1557        if alpha is not None:
1558            self.alpha(alpha)
1559        return self
1560
1561    def clean(self):
1562        """
1563        Clean pointcloud or mesh by removing coincident points.
1564        """
1565        cpd = vtk.vtkCleanPolyData()
1566        cpd.PointMergingOn()
1567        cpd.ConvertLinesToPointsOn()
1568        cpd.ConvertPolysToLinesOn()
1569        cpd.ConvertStripsToPolysOn()
1570        cpd.SetInputData(self.inputdata())
1571        cpd.Update()
1572        out = self._update(cpd.GetOutput())
1573
1574        out.pipeline = utils.OperationNode(
1575            "clean", parents=[self], 
1576            comment=f"#pts {out._data.GetNumberOfPoints()}",
1577        )
1578        return out
1579
1580    def subsample(self, fraction, absolute=False):
1581        """
1582        Subsample a point cloud by requiring that the points
1583        or vertices are far apart at least by the specified fraction of the object size.
1584        If a Mesh is passed the polygonal faces are not removed
1585        but holes can appear as vertices are removed.
1586
1587        Examples:
1588            - [moving_least_squares1D.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/moving_least_squares1D.py)
1589
1590                ![](https://vedo.embl.es/images/advanced/moving_least_squares1D.png)
1591
1592            - [recosurface.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/recosurface.py)
1593
1594                ![](https://vedo.embl.es/images/advanced/recosurface.png)
1595        """
1596        if not absolute:
1597            if fraction > 1:
1598                vedo.logger.warning(f"subsample(fraction=...), fraction must be < 1, but is {fraction}")
1599            if fraction <= 0:
1600                return self
1601
1602        cpd = vtk.vtkCleanPolyData()
1603        cpd.PointMergingOn()
1604        cpd.ConvertLinesToPointsOn()
1605        cpd.ConvertPolysToLinesOn()
1606        cpd.ConvertStripsToPolysOn()
1607        cpd.SetInputData(self.inputdata())
1608        if absolute:
1609            cpd.SetTolerance(fraction/self.diagonal_size())
1610            # cpd.SetToleranceIsAbsolute(absolute)
1611        else:
1612            cpd.SetTolerance(fraction)
1613        cpd.Update()
1614
1615        ps = 2
1616        if self.GetProperty().GetRepresentation() == 0:
1617            ps = self.GetProperty().GetPointSize()
1618        
1619        out = self._update(cpd.GetOutput()).ps(ps)
1620
1621        out.pipeline = utils.OperationNode(
1622            "subsample", parents=[self],
1623            comment=f"#pts {out._data.GetNumberOfPoints()}",
1624        )
1625        return out
1626
1627    def threshold(self, scalars, above=None, below=None, on="points"):
1628        """
1629        Extracts cells where scalar value satisfies threshold criterion.
1630
1631        Arguments:
1632            scalars : (str)
1633                name of the scalars array.
1634            above : (float)
1635                minimum value of the scalar
1636            below : (float)
1637                maximum value of the scalar
1638            on : (str)
1639                if 'cells' assume array of scalars refers to cell data.
1640
1641        Examples:
1642            - [mesh_threshold.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/mesh_threshold.py)
1643        """
1644        thres = vtk.vtkThreshold()
1645        thres.SetInputData(self.inputdata())
1646
1647        if on.startswith("c"):
1648            asso = vtk.vtkDataObject.FIELD_ASSOCIATION_CELLS
1649        else:
1650            asso = vtk.vtkDataObject.FIELD_ASSOCIATION_POINTS
1651
1652        thres.SetInputArrayToProcess(0, 0, 0, asso, scalars)
1653
1654        if above is None and below is not None:
1655            thres.ThresholdByLower(below)
1656        elif below is None and above is not None:
1657            thres.ThresholdByUpper(above)
1658        else:
1659            thres.ThresholdBetween(above, below)
1660        thres.Update()
1661
1662        gf = vtk.vtkGeometryFilter()
1663        gf.SetInputData(thres.GetOutput())
1664        gf.Update()
1665        return self._update(gf.GetOutput())
1666
1667    def quantize(self, value):
1668        """
1669        The user should input a value and all {x,y,z} coordinates
1670        will be quantized to that absolute grain size.
1671        """
1672        poly = self.inputdata()
1673        qp = vtk.vtkQuantizePolyDataPoints()
1674        qp.SetInputData(poly)
1675        qp.SetQFactor(value)
1676        qp.Update()
1677        out = self._update(qp.GetOutput()).flat()
1678        out.pipeline = utils.OperationNode("quantize", parents=[self])
1679        return out
1680
1681    def average_size(self):
1682        """
1683        Calculate the average size of a mesh.
1684        This is the mean of the vertex distances from the center of mass.
1685        """
1686        coords = self.points()
1687        cm = np.mean(coords, axis=0)
1688        if coords.shape[0] == 0:
1689            return 0.0
1690        cc = coords - cm
1691        return np.mean(np.linalg.norm(cc, axis=1))
1692
1693
1694    @deprecated(reason=vedo.colors.red + "Please use center_of_mass()" + vedo.colors.reset)
1695    def centerOfMass(self):
1696        "Deprecated. Please use `center_of_mass()`"
1697        return self.center_of_mass()
1698
1699    def center_of_mass(self):
1700        """Get the center of mass of mesh."""
1701        cmf = vtk.vtkCenterOfMass()
1702        cmf.SetInputData(self.polydata())
1703        cmf.Update()
1704        c = cmf.GetCenter()
1705        return np.array(c)
1706
1707    def normal_at(self, i):
1708        """Return the normal vector at vertex point `i`."""
1709        normals = self.polydata().GetPointData().GetNormals()
1710        return np.array(normals.GetTuple(i))
1711
1712    def normals(self, cells=False, recompute=True):
1713        """Retrieve vertex normals as a numpy array.
1714
1715        Arguments:
1716            cells : (bool)
1717                if `True` return cell normals.
1718
1719            recompute : (bool)
1720                if `True` normals are recalculated if not already present.
1721                Note that this might modify the number of mesh points.
1722        """
1723        if cells:
1724            vtknormals = self.polydata().GetCellData().GetNormals()
1725        else:
1726            vtknormals = self.polydata().GetPointData().GetNormals()
1727        if not vtknormals and recompute:
1728            try:
1729                self.compute_normals(cells=cells)
1730                if cells:
1731                    vtknormals = self.polydata().GetCellData().GetNormals()
1732                else:
1733                    vtknormals = self.polydata().GetPointData().GetNormals()
1734            except AttributeError:
1735                # can be that 'Points' object has no attribute 'compute_normals'
1736                pass
1737
1738        if not vtknormals:
1739            return np.array([])
1740        return utils.vtk2numpy(vtknormals)
1741
1742    def labels(
1743        self,
1744        content=None,
1745        on="points",
1746        scale=None,
1747        xrot=0,
1748        yrot=0,
1749        zrot=0,
1750        ratio=1,
1751        precision=None,
1752        italic=False,
1753        font="",
1754        justify="bottom-left",
1755        c="black",
1756        alpha=1,
1757        cells=None,
1758    ):
1759        """
1760        Generate value or ID labels for mesh cells or points.
1761        For large nr. of labels use `font="VTK"` which is much faster.
1762
1763        See also:
1764            `labels2d()`, `flagpole()`, `caption()` and `legend()`.
1765
1766        Arguments:
1767            content : (list,int,str)
1768                either 'id', 'cellid', array name or array number.
1769                A array can also be passed (must match the nr. of points or cells).
1770            on : (str)
1771                generate labels for "cells" instead of "points"
1772            scale : (float)
1773                absolute size of labels, if left as None it is automatic
1774            zrot : (float)
1775                local rotation angle of label in degrees
1776            ratio : (int)
1777                skipping ratio, to reduce nr of labels for large meshes
1778            precision : (int)
1779                numeric precision of labels
1780
1781        ```python
1782        from vedo import *
1783        s = Sphere(res=10).linewidth(1).c("orange").compute_normals()
1784        point_ids = s.labels('id', on="points").c('green')
1785        cell_ids  = s.labels('id', on="cells" ).c('black')
1786        show(s, point_ids, cell_ids)
1787        ```
1788        ![](https://vedo.embl.es/images/feats/labels.png)
1789
1790        Examples:
1791            - [boundaries.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/boundaries.py)
1792
1793                ![](https://vedo.embl.es/images/basic/boundaries.png)
1794        """
1795        if cells is not None:  # deprecation message
1796            vedo.logger.warning("In labels(cells=...) please use labels(on='cells') instead")
1797
1798        if "cell" in on or "face" in on:
1799            cells = True
1800
1801        if isinstance(content, str):
1802            if content in ('cellid', 'cellsid'):
1803                cells = True
1804                content = "id"
1805
1806        if cells:
1807            elems = self.cell_centers()
1808            norms = self.normals(cells=True, recompute=False)
1809            ns = np.sqrt(self.ncells)
1810        else:
1811            elems = self.points()
1812            norms = self.normals(cells=False, recompute=False)
1813            ns = np.sqrt(self.npoints)
1814
1815        hasnorms = False
1816        if len(norms):
1817            hasnorms = True
1818
1819        if scale is None:
1820            if not ns:
1821                ns = 100
1822            scale = self.diagonal_size() / ns / 10
1823
1824        arr = None
1825        mode = 0
1826        if content is None:
1827            mode = 0
1828            if cells:
1829                if self.inputdata().GetCellData().GetScalars():
1830                    name = self.inputdata().GetCellData().GetScalars().GetName()
1831                    arr = self.celldata[name]
1832            else:
1833                if self.inputdata().GetPointData().GetScalars():
1834                    name = self.inputdata().GetPointData().GetScalars().GetName()
1835                    arr = self.pointdata[name]
1836        elif isinstance(content, (str, int)):
1837            if content == "id":
1838                mode = 1
1839            elif cells:
1840                mode = 0
1841                arr = self.celldata[content]
1842            else:
1843                mode = 0
1844                arr = self.pointdata[content]
1845        elif utils.is_sequence(content):
1846            mode = 0
1847            arr = content
1848            # print('WEIRD labels() test', content)
1849            # exit()
1850
1851        if arr is None and mode == 0:
1852            vedo.logger.error("in labels(), array not found for points or cells")
1853            return None
1854
1855        tapp = vtk.vtkAppendPolyData()
1856        ninputs = 0
1857
1858        for i, e in enumerate(elems):
1859            if i % ratio:
1860                continue
1861
1862            if mode == 1:
1863                txt_lab = str(i)
1864            else:
1865                if precision:
1866                    txt_lab = utils.precision(arr[i], precision)
1867                else:
1868                    txt_lab = str(arr[i])
1869
1870            if not txt_lab:
1871                continue
1872
1873            if font == "VTK":
1874                tx = vtk.vtkVectorText()
1875                tx.SetText(txt_lab)
1876                tx.Update()
1877                tx_poly = tx.GetOutput()
1878            else:
1879                tx_poly = vedo.shapes.Text3D(txt_lab, font=font, justify=justify)
1880                tx_poly = tx_poly.inputdata()
1881
1882            if tx_poly.GetNumberOfPoints() == 0:
1883                continue  #######################
1884            ninputs += 1
1885
1886            T = vtk.vtkTransform()
1887            T.PostMultiply()
1888            if italic:
1889                T.Concatenate([1,0.2,0,0,
1890                               0,1,0,0,
1891                               0,0,1,0,
1892                               0,0,0,1])
1893            if hasnorms:
1894                ni = norms[i]
1895                if cells:  # center-justify
1896                    bb = tx_poly.GetBounds()
1897                    dx, dy = (bb[1] - bb[0]) / 2, (bb[3] - bb[2]) / 2
1898                    T.Translate(-dx, -dy, 0)
1899                if xrot:
1900                    T.RotateX(xrot)
1901                if yrot:
1902                    T.RotateY(yrot)
1903                if zrot:
1904                    T.RotateZ(zrot)
1905                crossvec = np.cross([0, 0, 1], ni)
1906                angle = np.arccos(np.dot([0, 0, 1], ni)) * 57.3
1907                T.RotateWXYZ(angle, crossvec)
1908                if cells:  # small offset along normal only for cells
1909                    T.Translate(ni * scale / 2)
1910            else:
1911                if xrot:
1912                    T.RotateX(xrot)
1913                if yrot:
1914                    T.RotateY(yrot)
1915                if zrot:
1916                    T.RotateZ(zrot)
1917            T.Scale(scale, scale, scale)
1918            T.Translate(e)
1919            tf = vtk.vtkTransformPolyDataFilter()
1920            tf.SetInputData(tx_poly)
1921            tf.SetTransform(T)
1922            tf.Update()
1923            tapp.AddInputData(tf.GetOutput())
1924
1925        if ninputs:
1926            tapp.Update()
1927            lpoly = tapp.GetOutput()
1928        else:  # return an empty obj
1929            lpoly = vtk.vtkPolyData()
1930
1931        ids = vedo.mesh.Mesh(lpoly, c=c, alpha=alpha)
1932        ids.GetProperty().LightingOff()
1933        ids.SetUseBounds(False)
1934        return ids
1935
1936    def labels2d(
1937            self,
1938            content="id",
1939            on="points",
1940            scale=1,
1941            precision=4,
1942            font="Calco",
1943            justify="bottom-left",
1944            angle=0,
1945            frame=False,
1946            c="black",
1947            bc=None,
1948            alpha=1,
1949        ):
1950        """
1951        Generate value or ID bi-dimensional labels for mesh cells or points.
1952
1953        See also: `labels()`, `flagpole()`, `caption()` and `legend()`.
1954
1955        Arguments:
1956            content : (str)
1957                either 'id', 'cellid', or array name
1958            on : (str)
1959                generate labels for "cells" instead of "points" (the default)
1960            scale : (float)
1961                size scaling of labels
1962            precision : (int)
1963                precision of numeric labels
1964            angle : (float)
1965                local rotation angle of label in degrees
1966            frame : (bool)
1967                draw a frame around the label
1968            bc : (str)
1969                background color of the label
1970
1971        ```python
1972        from vedo import Sphere, show
1973        sph = Sphere(quads=True, res=4).compute_normals().wireframe()
1974        sph.celldata["zvals"] = sph.cell_centers()[:,2]
1975        l2d = sph.labels("zvals", on="cells", precision=2).backcolor('orange9')
1976        show(sph, l2d, axes=1).close()
1977        ```
1978        ![](https://vedo.embl.es/images/feats/labels2d.png)
1979        """
1980        cells = False
1981        if isinstance(content, str):
1982            if content in ('cellid', 'cellsid'):
1983                cells = True
1984                content = "id"
1985        
1986        if "cell" in on:
1987            cells = True
1988        elif "point" in on:
1989            cells = False
1990
1991        if cells:
1992            if content != 'id' and content not in self.celldata.keys():
1993                vedo.logger.error(f"In labels2D: cell array {content} does not exist.")
1994                return None
1995            cellcloud = Points(self.cell_centers())
1996            arr = self.inputdata().GetCellData().GetScalars()
1997            poly = cellcloud.polydata(False)
1998            poly.GetPointData().SetScalars(arr)
1999        else:
2000            poly = self.polydata()
2001            if content != 'id' and content not in self.pointdata.keys():
2002                vedo.logger.error(f"In labels2D: point array {content} does not exist.")
2003                return None
2004            self.pointdata.select(content)
2005
2006        mp = vtk.vtkLabeledDataMapper()
2007
2008        if content == 'id':
2009            mp.SetLabelModeToLabelIds()
2010        else:
2011            mp.SetLabelModeToLabelScalars()
2012            if precision is not None:
2013                mp.SetLabelFormat(f'%-#.{precision}g')
2014
2015        pr = mp.GetLabelTextProperty()
2016        c = colors.get_color(c)
2017        pr.SetColor(c)
2018        pr.SetOpacity(alpha)
2019        pr.SetFrame(frame)
2020        pr.SetFrameColor(c)
2021        pr.SetItalic(False)
2022        pr.BoldOff()
2023        pr.ShadowOff()
2024        pr.UseTightBoundingBoxOn()
2025        pr.SetOrientation(angle)
2026        pr.SetFontFamily(vtk.VTK_FONT_FILE)
2027        fl = utils.get_font_path(font)
2028        pr.SetFontFile(fl)
2029        pr.SetFontSize(int(20*scale))
2030
2031        if "cent" in justify or "mid" in justify:
2032            pr.SetJustificationToCentered()
2033        elif "rig" in justify:
2034            pr.SetJustificationToRight()
2035        elif "left" in justify:
2036            pr.SetJustificationToLeft()
2037        # ------
2038        if "top" in justify:
2039            pr.SetVerticalJustificationToTop()
2040        else:
2041            pr.SetVerticalJustificationToBottom()
2042
2043        if bc is not None:
2044            bc = colors.get_color(bc)
2045            pr.SetBackgroundColor(bc)
2046            pr.SetBackgroundOpacity(alpha)
2047
2048        mp.SetInputData(poly)
2049        a2d = vtk.vtkActor2D()
2050        a2d.SetMapper(mp)
2051        return a2d
2052
2053    def legend(self, txt):
2054        """Book a legend text."""
2055        self.info["legend"] = txt
2056        return self
2057
2058    @deprecated(reason=vedo.colors.red + "Please use flagpole() instead" + vedo.colors.reset)
2059    def vignette(self, *args, **kwargs):
2060        """Deprecated. Use `flagpole()`."""
2061        return self.flagpole(*args, **kwargs)
2062
2063    def flagpole(
2064        self,
2065        txt=None,
2066        point=None,
2067        offset=None,
2068        s=None,
2069        font="",
2070        rounded=True,
2071        c=None,
2072        alpha=1,
2073        lw=2,
2074        italic=0,
2075        padding=0.1,
2076    ):
2077        """
2078        Generate a flag pole style element to describe an object.
2079        Returns a `Mesh` object.
2080
2081        Use flagpole.follow_camera() to make it face the camera in the scene.
2082
2083        See also `flagpost()`.
2084
2085        Arguments:
2086            txt : (str)
2087                Text to display. The default is the filename or the object name.
2088            point : (list)
2089                position of the flagpole pointer. The default is None.
2090            offset : (list)
2091                text offset wrt the application point. The default is None.
2092            s : (float)
2093                size of the flagpole. The default is None.
2094            font : (str)
2095                font face. Check [available fonts here](https://vedo.embl.es/fonts).
2096            rounded : (bool)
2097                draw a rounded or squared box around the text. The default is True.
2098            c : (list)
2099                text and box color. The default is None.
2100            alpha : (float)
2101                opacity of text and box. The default is 1.
2102            lw : (float)
2103                line with of box frame. The default is 2.
2104            italic : (float)
2105                italicness of text. The default is 0.
2106
2107        Examples:
2108            - [intersect2d.py](https://github.com/marcomusy/vedo/tree/master/examples/pyplot/intersect2d.py)
2109
2110                ![](https://vedo.embl.es/images/pyplot/intersect2d.png)
2111
2112            - [goniometer.py](https://github.com/marcomusy/vedo/tree/master/examples/pyplot/goniometer.py)
2113            - [flag_labels1.py](https://github.com/marcomusy/vedo/tree/master/examples/other/flag_labels1.py)
2114            - [flag_labels2.py](https://github.com/marcomusy/vedo/tree/master/examples/other/flag_labels2.py)
2115        """
2116        acts = []
2117
2118        if txt is None:
2119            if self.filename:
2120                txt = self.filename.split("/")[-1]
2121            elif self.name:
2122                txt = self.name
2123            else:
2124                return None
2125
2126        x0, x1, y0, y1, z0, z1 = self.bounds()
2127        d = self.diagonal_size()
2128        if point is None:
2129            if d:
2130                point = self.closest_point([(x0 + x1) / 2, (y0 + y1) / 2, z1])
2131            else:  # it's a Point
2132                point = self.GetPosition()
2133
2134        pt = utils.make3d(point)
2135
2136        if offset is None:
2137            offset = [(x1 - x0) / 2, (y1 - y0) / 6, 0]
2138        offset = utils.make3d(offset)
2139
2140        if s is None:
2141            s = d / 20
2142
2143        sph = None
2144        if d and (z1 - z0) / d > 0.1:
2145            sph = vedo.shapes.Sphere(pt, r=s * 0.4, res=6)
2146
2147        if c is None:
2148            c = np.array(self.color()) / 1.4
2149
2150        lb = vedo.shapes.Text3D(
2151            txt, pos=pt+offset, s=s, font=font,
2152            italic=italic, justify="center-left",
2153        )
2154        acts.append(lb)
2155
2156        if d and not sph:
2157            sph = vedo.shapes.Circle(pt, r=s / 3, res=15)
2158        acts.append(sph)
2159
2160        x0, x1, y0, y1, z0, z1 = lb.GetBounds()
2161        if rounded:
2162            box = vedo.shapes.KSpline(
2163                [(x0, y0, z0), (x1, y0, z0), (x1, y1, z0), (x0, y1, z0)], closed=True
2164            )
2165        else:
2166            box = vedo.shapes.Line(
2167                [(x0,y0,z0), (x1,y0,z0), (x1,y1,z0), (x0,y1,z0), (x0,y0,z0)]
2168            )
2169
2170        cnt = [(x0+x1) / 2, (y0+y1) / 2, (z0+z1) / 2]
2171
2172        box.SetOrigin(cnt)
2173        box.scale([1 + padding, 1 + 2*padding, 1])
2174        acts.append(box)
2175
2176        # pts = box.points()
2177        # bfaces = []
2178        # for i, pt in enumerate(pts):
2179        #     if i:
2180        #         face = [i-1, i, 0]
2181        #         bfaces.append(face)
2182        # bpts = [cnt] + pts.tolist()
2183        # box2 = vedo.Mesh([bpts, bfaces]).z(-cnt[0]/10)#.c('w').alpha(0.1)
2184        # #should be made assembly otherwise later merge() nullifies it
2185        # box2.SetOrigin(cnt)
2186        # acts.append(box2)
2187
2188        x0, x1, y0, y1, z0, z1 = box.bounds()
2189        if x0 < pt[0] < x1:
2190            c0 = box.closest_point(pt)
2191            c1 = [c0[0], c0[1] + (pt[1] - y0) / 4, pt[2]]
2192        elif (pt[0] - x0) < (x1 - pt[0]):
2193            c0 = [x0, (y0 + y1) / 2, pt[2]]
2194            c1 = [x0 + (pt[0] - x0) / 4, (y0 + y1) / 2, pt[2]]
2195        else:
2196            c0 = [x1, (y0 + y1) / 2, pt[2]]
2197            c1 = [x1 + (pt[0] - x1) / 4, (y0 + y1) / 2, pt[2]]
2198
2199        con = vedo.shapes.Line([c0, c1, pt])
2200        acts.append(con)
2201
2202        macts = vedo.merge(acts).c(c).alpha(alpha)
2203        macts.SetOrigin(pt)
2204        macts.bc("t").pickable(False).GetProperty().LightingOff()
2205        macts.GetProperty().SetLineWidth(lw)
2206        macts.UseBoundsOff()
2207        macts.name = "FlagPole"
2208        return macts
2209
2210    def flagpost(
2211        self,
2212        txt=None,
2213        point=None,
2214        offset=None,
2215        s=1,
2216        c="k9",
2217        bc="k1",
2218        alpha=1,
2219        lw=0,
2220        font="Calco",
2221        justify="center-left",
2222        vspacing=1,
2223    ):
2224        """
2225        Generate a flag post style element to describe an object.
2226
2227        Arguments:
2228            txt : (str)
2229                Text to display. The default is the filename or the object name.
2230            point : (list)
2231                position of the flag anchor point. The default is None.
2232            offset : (list)
2233                a 3D displacement or offset. The default is None.
2234            s : (float)
2235                size of the text to be shown
2236            c : (list)
2237                color of text and line
2238            bc : (list)
2239                color of the flag background
2240            alpha : (float)
2241                opacity of text and box.
2242            lw : (int)
2243                line with of box frame. The default is 0.
2244            font : (str)
2245                font name. Use a monospace font for better rendering. The default is "Calco".
2246                Type `vedo -r fonts` for a font demo.
2247                Check [available fonts here](https://vedo.embl.es/fonts).
2248            justify : (str)
2249                internal text justification. The default is "center-left".
2250            vspacing : (float)
2251                vertical spacing between lines.
2252        
2253        Examples:
2254            - [flag_labels2.py](https://github.com/marcomusy/vedo/tree/master/examples/examples/other/flag_labels2.py)
2255
2256            ![](https://vedo.embl.es/images/other/flag_labels2.png)
2257        """
2258        if txt is None:
2259            if self.filename:
2260                txt = self.filename.split("/")[-1]
2261            elif self.name:
2262                txt = self.name
2263            else:
2264                return None
2265
2266        x0, x1, y0, y1, z0, z1 = self.bounds()
2267        d = self.diagonal_size()
2268        if point is None:
2269            if d:
2270                point = self.closest_point([(x0 + x1) / 2, (y0 + y1) / 2, z1])
2271            else:  # it's a Point
2272                point = self.GetPosition()
2273
2274        point = utils.make3d(point)
2275
2276        if offset is None:
2277            offset = [0,0, (z1-z0)/2]
2278        offset = utils.make3d(offset)
2279
2280        fpost = vedo.addons.Flagpost(
2281            txt, 
2282            point, 
2283            point + offset,
2284            s,
2285            c,
2286            bc,
2287            alpha,
2288            lw,
2289            font,
2290            justify,
2291            vspacing,
2292        )
2293        self._caption = fpost
2294        return fpost
2295
2296    def caption(
2297        self,
2298        txt=None,
2299        point=None,
2300        size=(0.30, 0.15),
2301        padding=5,
2302        font="Calco",
2303        justify="center-right",
2304        vspacing=1,
2305        c=None,
2306        alpha=1,
2307        lw=1,
2308        ontop=True,
2309    ):
2310        """
2311        Add a 2D caption to an object which follows the camera movements.
2312        Latex is not supported. Returns the same input object for concatenation.
2313
2314        See also `flagpole()`, `flagpost()`, `labels()` and `legend()`
2315        with similar functionality.
2316
2317        Arguments:
2318            txt : (str)
2319                text to be rendered. The default is the file name.
2320            point : (list)
2321                anchoring point. The default is None.
2322            size : (list)
2323                (width, height) of the caption box. The default is (0.30, 0.15).
2324            padding : (float)
2325                padding space of the caption box in pixels. The default is 5.
2326            font : (str)
2327                font name. Use a monospace font for better rendering. The default is "VictorMono".
2328                Type `vedo -r fonts` for a font demo.
2329                Check [available fonts here](https://vedo.embl.es/fonts).
2330            justify : (str)
2331                internal text justification. The default is "center-right".
2332            vspacing : (float)
2333                vertical spacing between lines. The default is 1.
2334            c : (str)
2335                text and box color. The default is 'lb'.
2336            alpha : (float)
2337                text and box transparency. The default is 1.
2338            lw : (int)
2339                line width in pixels. The default is 1.
2340            ontop : (bool)
2341                keep the 2d caption always on top. The default is True.
2342
2343        Examples:
2344            - [caption.py](https://github.com/marcomusy/vedo/tree/master/examples/pyplot/caption.py)
2345
2346                ![](https://vedo.embl.es/images/pyplot/caption.png)
2347
2348            - [flag_labels1.py](https://github.com/marcomusy/vedo/tree/master/examples/other/flag_labels1.py)
2349            - [flag_labels2.py](https://github.com/marcomusy/vedo/tree/master/examples/other/flag_labels2.py)
2350        """
2351        if txt is None:
2352            if self.filename:
2353                txt = self.filename.split("/")[-1]
2354            elif self.name:
2355                txt = self.name
2356
2357        if not txt:  # disable it
2358            self._caption = None
2359            return self
2360
2361        for r in vedo.shapes._reps:
2362            txt = txt.replace(r[0], r[1])
2363
2364        if c is None:
2365            c = np.array(self.GetProperty().GetColor()) / 2
2366        else:
2367            c = colors.get_color(c)
2368
2369        if point is None:
2370            x0, x1, y0, y1, _, z1 = self.GetBounds()
2371            pt = [(x0 + x1) / 2, (y0 + y1) / 2, z1]
2372            point = self.closest_point(pt)
2373
2374        capt = vtk.vtkCaptionActor2D()
2375        capt.SetAttachmentPoint(point)
2376        capt.SetBorder(True)
2377        capt.SetLeader(True)
2378        sph = vtk.vtkSphereSource()
2379        sph.Update()
2380        capt.SetLeaderGlyphData(sph.GetOutput())
2381        capt.SetMaximumLeaderGlyphSize(5)
2382        capt.SetPadding(padding)
2383        capt.SetCaption(txt)
2384        capt.SetWidth(size[0])
2385        capt.SetHeight(size[1])
2386        capt.SetThreeDimensionalLeader(not ontop)
2387
2388        pra = capt.GetProperty()
2389        pra.SetColor(c)
2390        pra.SetOpacity(alpha)
2391        pra.SetLineWidth(lw)
2392
2393        pr = capt.GetCaptionTextProperty()
2394        pr.SetFontFamily(vtk.VTK_FONT_FILE)
2395        fl = utils.get_font_path(font)
2396        pr.SetFontFile(fl)
2397        pr.ShadowOff()
2398        pr.BoldOff()
2399        pr.FrameOff()
2400        pr.SetColor(c)
2401        pr.SetOpacity(alpha)
2402        pr.SetJustificationToLeft()
2403        if "top" in justify:
2404            pr.SetVerticalJustificationToTop()
2405        if "bottom" in justify:
2406            pr.SetVerticalJustificationToBottom()
2407        if "cent" in justify:
2408            pr.SetVerticalJustificationToCentered()
2409            pr.SetJustificationToCentered()
2410        if "left" in justify:
2411            pr.SetJustificationToLeft()
2412        if "right" in justify:
2413            pr.SetJustificationToRight()
2414        pr.SetLineSpacing(vspacing)
2415        self._caption = capt
2416        return self
2417
2418    @deprecated(reason=vedo.colors.red + "Please use align_to()" + vedo.colors.reset)
2419    def alignTo(self, *a, **b):
2420        "Deprecated, Please use `align_to()`."
2421        return self.align_to(*a, **b)
2422
2423    def align_to(
2424            self,
2425            target,
2426            iters=100,
2427            rigid=False,
2428            invert=False,
2429            use_centroids=False,
2430        ):
2431        """
2432        Aligned to target mesh through the `Iterative Closest Point` algorithm.
2433
2434        The core of the algorithm is to match each vertex in one surface with
2435        the closest surface point on the other, then apply the transformation
2436        that modify one surface to best match the other (in the least-square sense).
2437
2438        Arguments:
2439            rigid : (bool)
2440                if True do not allow scaling
2441            invert : (bool)
2442                if True start by aligning the target to the source but
2443                invert the transformation finally. Useful when the target is smaller
2444                than the source.
2445            use_centroids : (bool)
2446                start by matching the centroids of the two objects.
2447
2448        Examples:
2449            - [align1.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/align1.py)
2450
2451                ![](https://vedo.embl.es/images/basic/align1.png)
2452
2453            - [align2.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/align2.py)
2454
2455                ![](https://vedo.embl.es/images/basic/align2.png)
2456        """
2457        icp = vtk.vtkIterativeClosestPointTransform()
2458        icp.SetSource(self.polydata())
2459        icp.SetTarget(target.polydata())
2460        if invert:
2461            icp.Inverse()
2462        icp.SetMaximumNumberOfIterations(iters)
2463        if rigid:
2464            icp.GetLandmarkTransform().SetModeToRigidBody()
2465        icp.SetStartByMatchingCentroids(use_centroids)
2466        icp.Update()
2467
2468        M = icp.GetMatrix()
2469        if invert:
2470            M.Invert()  # icp.GetInverse() doesnt work!
2471        # self.apply_transform(M)
2472        self.SetUserMatrix(M)
2473
2474        self.transform = self.GetUserTransform()
2475        self.point_locator = None
2476        self.cell_locator = None
2477        
2478        self.pipeline = utils.OperationNode(
2479            "align_to", parents=[self, target], comment=f"rigid = {rigid}",
2480        )
2481        return self
2482
2483    def transform_with_landmarks(
2484            self,
2485            source_landmarks,
2486            target_landmarks,
2487            rigid=False,
2488            affine=False,
2489            least_squares=False,
2490        ):
2491        """
2492        Transform mesh orientation and position based on a set of landmarks points.
2493        The algorithm finds the best matching of source points to target points
2494        in the mean least square sense, in one single step.
2495
2496        If affine is True the x, y and z axes can scale independently but stay collinear.
2497        With least_squares they can vary orientation.
2498
2499        Examples:
2500            - [align5.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/align5.py)
2501
2502                ![](https://vedo.embl.es/images/basic/align5.png)
2503        """
2504
2505        if utils.is_sequence(source_landmarks):
2506            ss = vtk.vtkPoints()
2507            for p in source_landmarks:
2508                ss.InsertNextPoint(p)
2509        else:
2510            ss = source_landmarks.polydata().GetPoints()
2511            if least_squares:
2512                source_landmarks = source_landmarks.points()
2513
2514        if utils.is_sequence(target_landmarks):
2515            st = vtk.vtkPoints()
2516            for p in target_landmarks:
2517                st.InsertNextPoint(p)
2518        else:
2519            st = target_landmarks.polydata().GetPoints()
2520            if least_squares:
2521                target_landmarks = target_landmarks.points()
2522
2523        if ss.GetNumberOfPoints() != st.GetNumberOfPoints():
2524            n1 = ss.GetNumberOfPoints()
2525            n2 = st.GetNumberOfPoints()
2526            vedo.logger.error(f"source and target have different nr of points {n1} vs {n2}")
2527            raise RuntimeError()
2528
2529        lmt = vtk.vtkLandmarkTransform()
2530        lmt.SetSourceLandmarks(ss)
2531        lmt.SetTargetLandmarks(st)
2532        lmt.SetModeToSimilarity()
2533        if rigid:
2534            lmt.SetModeToRigidBody()
2535            lmt.Update()
2536            self.SetUserTransform(lmt)
2537
2538        elif affine:
2539            lmt.SetModeToAffine()
2540            lmt.Update()
2541            self.SetUserTransform(lmt)
2542
2543        elif least_squares:
2544            cms = source_landmarks.mean(axis=0)
2545            cmt = target_landmarks.mean(axis=0)
2546            m = np.linalg.lstsq(source_landmarks-cms, target_landmarks-cmt, rcond=None)[0]
2547            M = vtk.vtkMatrix4x4()
2548            for i in range(3):
2549                for j in range(3):
2550                    M.SetElement(j, i, m[i][j])
2551            lmt = vtk.vtkTransform()
2552            lmt.Translate( cmt)
2553            lmt.Concatenate(M)
2554            lmt.Translate(-cms)
2555            self.apply_transform(lmt, concatenate=True)
2556        else:
2557            self.SetUserTransform(lmt)
2558
2559        self.transform = lmt
2560        self.point_locator = None
2561        self.cell_locator = None
2562        self.pipeline = utils.OperationNode("transform_with_landmarks", parents=[self])
2563        return self
2564
2565    @deprecated(reason=vedo.colors.red + "Please use apply_transform()" + vedo.colors.reset)
2566    def applyTransform(self, *a, **b):
2567        "Deprecated. Please use `apply_transform()`."
2568        return self.apply_transform(*a, **b)
2569
2570    def apply_transform(self, T, reset=False, concatenate=False):
2571        """
2572        Apply a linear or non-linear transformation to the mesh polygonal data.
2573
2574        Arguments:
2575            T : (matrix)
2576                `vtkTransform`, `vtkMatrix4x4` or a 4x4 or 3x3 python or numpy matrix.
2577            reset : (bool)
2578                if True reset the current transformation matrix
2579                to identity after having moved the object, otherwise the internal
2580                matrix will stay the same (to only affect visualization).
2581                It the input transformation has no internal defined matrix (ie. non linear)
2582                then reset will be assumed as True.
2583            concatenate : (bool)
2584                concatenate the transformation with the current existing one
2585
2586        Example:
2587            ```python
2588            from vedo import Cube, show
2589            c1 = Cube().rotate_z(5).x(2).y(1)
2590            print("cube1 position", c1.pos())
2591            T = c1.get_transform()  # rotate by 5 degrees, sum 2 to x and 1 to y
2592            c2 = Cube().c('r4')
2593            c2.apply_transform(T)   # ignore previous movements
2594            c2.apply_transform(T, concatenate=True)
2595            c2.apply_transform(T, concatenate=True)
2596            print("cube2 position", c2.pos())
2597            show(c1, c2, axes=1).close()
2598            ```
2599            ![](https://vedo.embl.es/images/feats/apply_transform.png)
2600        """
2601        self.point_locator = None
2602        self.cell_locator = None
2603
2604        if isinstance(T, vtk.vtkMatrix4x4):
2605            tr = vtk.vtkTransform()
2606            tr.SetMatrix(T)
2607            T = tr
2608
2609        elif utils.is_sequence(T):
2610            M = vtk.vtkMatrix4x4()
2611            n = len(T[0])
2612            for i in range(n):
2613                for j in range(n):
2614                    M.SetElement(i, j, T[i][j])
2615            tr = vtk.vtkTransform()
2616            tr.SetMatrix(M)
2617            T = tr
2618
2619        if reset or not hasattr(T, 'GetScale'): # might be non-linear
2620
2621            tf = vtk.vtkTransformPolyDataFilter()
2622            tf.SetTransform(T)
2623            tf.SetInputData(self.polydata())
2624            tf.Update()
2625
2626            I = vtk.vtkMatrix4x4()
2627            self.PokeMatrix(I)  # reset to identity
2628            self.SetUserTransform(None)
2629
2630            self._update(tf.GetOutput()) ### UPDATE
2631            self.transform = T
2632
2633        else:
2634
2635            if concatenate:
2636
2637                M = vtk.vtkTransform()
2638                M.PostMultiply()
2639                M.SetMatrix(self.GetMatrix())
2640
2641                M.Concatenate(T)
2642
2643                self.SetScale(M.GetScale())
2644                self.SetOrientation(M.GetOrientation())
2645                self.SetPosition(M.GetPosition())
2646                self.transform = M
2647                self.SetUserTransform(None)
2648
2649            else:
2650
2651                self.SetScale(T.GetScale())
2652                self.SetOrientation(T.GetOrientation())
2653                self.SetPosition(T.GetPosition())
2654                self.SetUserTransform(None)
2655
2656                self.transform = T
2657
2658        return self
2659
2660    def normalize(self):
2661        """Scale Mesh average size to unit."""
2662        coords = self.points()
2663        if not coords.shape[0]:
2664            return self
2665        cm = np.mean(coords, axis=0)
2666        pts = coords - cm
2667        xyz2 = np.sum(pts * pts, axis=0)
2668        scale = 1 / np.sqrt(np.sum(xyz2) / len(pts))
2669        t = vtk.vtkTransform()
2670        t.PostMultiply()
2671        # t.Translate(-cm)
2672        t.Scale(scale, scale, scale)
2673        # t.Translate(cm)
2674        tf = vtk.vtkTransformPolyDataFilter()
2675        tf.SetInputData(self.inputdata())
2676        tf.SetTransform(t)
2677        tf.Update()
2678        self.point_locator = None
2679        self.cell_locator = None
2680        return self._update(tf.GetOutput())
2681
2682    def mirror(self, axis="x", origin=(0, 0, 0), reset=False):
2683        """
2684        Mirror the mesh  along one of the cartesian axes
2685
2686        Arguments:
2687            axis : (str)
2688                axis to use for mirroring, must be set to x, y, z or n.
2689                Or any combination of those. Adding 'n' reverses mesh faces (hence normals).
2690            origin : (list)
2691                use this point as the origin of the mirroring transformation.
2692            reset : (bool)
2693                if True keep into account the current position of the object,
2694                and then reset its internal transformation matrix to Identity.
2695
2696        Examples:
2697            - [mirror.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/mirror.py)
2698
2699                ![](https://vedo.embl.es/images/basic/mirror.png)
2700        """
2701        sx, sy, sz = 1, 1, 1
2702        if "x" in axis.lower(): sx = -1
2703        if "y" in axis.lower(): sy = -1
2704        if "z" in axis.lower(): sz = -1
2705        origin = np.array(origin)
2706        tr = vtk.vtkTransform()
2707        tr.PostMultiply()
2708        tr.Translate(-origin)
2709        tr.Scale(sx, sy, sz)
2710        tr.Translate(origin)
2711        tf = vtk.vtkTransformPolyDataFilter()
2712        tf.SetInputData(self.polydata(reset))
2713        tf.SetTransform(tr)
2714        tf.Update()
2715        outpoly = tf.GetOutput()
2716        if reset:
2717            self.PokeMatrix(vtk.vtkMatrix4x4())  # reset to identity
2718        if sx * sy * sz < 0 or "n" in axis:
2719            rs = vtk.vtkReverseSense()
2720            rs.SetInputData(outpoly)
2721            rs.ReverseNormalsOff()
2722            rs.Update()
2723            outpoly = rs.GetOutput()
2724
2725        self.point_locator = None
2726        self.cell_locator = None
2727
2728        out = self._update(outpoly)
2729
2730        out.pipeline = utils.OperationNode(
2731            f"mirror\naxis = {axis}", parents=[self])
2732        return out
2733
2734    def shear(self, x=0, y=0, z=0):
2735        """Apply a shear deformation along one of the main axes"""
2736        t = vtk.vtkTransform()
2737        sx, sy, sz = self.GetScale()
2738        t.SetMatrix([sx, x, 0, 0,
2739                      y,sy, z, 0,
2740                      0, 0,sz, 0,
2741                      0, 0, 0, 1])
2742        self.apply_transform(t, reset=True)
2743        return self
2744
2745    def flip_normals(self):
2746        """Flip all mesh normals. Same as `mesh.mirror('n')`."""
2747        rs = vtk.vtkReverseSense()
2748        rs.SetInputData(self.inputdata())
2749        rs.ReverseCellsOff()
2750        rs.ReverseNormalsOn()
2751        rs.Update()
2752        out = self._update(rs.GetOutput())
2753        self.pipeline = utils.OperationNode("flip_normals", parents=[self])
2754        return out
2755
2756    #####################################################################################
2757    def cmap(
2758        self,
2759        input_cmap,
2760        input_array=None,
2761        on="points",
2762        name="Scalars",
2763        vmin=None,
2764        vmax=None,
2765        n_colors=256,
2766        alpha=1,
2767        logscale=False
2768    ):
2769        """
2770        Set individual point/cell colors by providing a list of scalar values and a color map.
2771
2772        Arguments:
2773            input_cmap : (str, list, vtkLookupTable, matplotlib.colors.LinearSegmentedColormap)
2774                color map scheme to transform a real number into a color.
2775            input_array : (str, list, vtkArray)
2776                can be the string name of an existing array, a numpy array or a `vtkArray`.
2777            on : (str)
2778                either 'points' or 'cells'.
2779                Apply the color map to data which is defined on either points or cells.
2780            name : (str)
2781                give a name to the provided numpy array (if input_array is a numpy array)
2782            vmin : (float)
2783                clip scalars to this minimum value
2784            vmax : (float)
2785                clip scalars to this maximum value
2786            n_colors : (int)
2787                number of distinct colors to be used in colormap table.
2788            alpha : (float, list)
2789                Mesh transparency. Can be a `list` of values one for each vertex.
2790            logscale : (bool)
2791                Use logscale
2792
2793        Examples:
2794            - [mesh_coloring.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/mesh_coloring.py)
2795            - [mesh_alphas.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/mesh_alphas.py)
2796            - [mesh_custom.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/mesh_custom.py)
2797            (and many others)
2798
2799                ![](https://vedo.embl.es/images/basic/mesh_custom.png)
2800        """
2801        self._cmap_name = input_cmap
2802        poly = self.inputdata()
2803
2804        if input_array is None:
2805            if not self.pointdata.keys() and self.celldata.keys():
2806                on = "cells"
2807                if not poly.GetCellData().GetScalars():
2808                    input_array = 0  # pick the first at hand
2809
2810        if on.startswith("point"):
2811            data = poly.GetPointData()
2812            n = poly.GetNumberOfPoints()
2813        elif on.startswith("cell"):
2814            data = poly.GetCellData()
2815            n = poly.GetNumberOfCells()
2816        else:
2817            vedo.logger.error("Must specify in cmap(on=...) to either 'cells' or 'points'")
2818            raise RuntimeError()
2819        
2820        if input_array is None:  # if None try to fetch the active scalars
2821            arr = data.GetScalars()
2822            if not arr:
2823                vedo.logger.error(f"in cmap(), cannot find any {on} active array ...skip coloring.")
2824                return self
2825
2826            if not arr.GetName():  # sometimes arrays dont have a name..
2827                arr.SetName(name)
2828            
2829        elif isinstance(input_array, str):  # if a string is passed
2830            arr = data.GetArray(input_array)
2831            if not arr:
2832                vedo.logger.error(f"in cmap(), cannot find {on} array {input_array} ...skip coloring.")
2833                return self
2834
2835        elif isinstance(input_array, int):  # if an int is passed
2836            if input_array < data.GetNumberOfArrays():
2837                arr = data.GetArray(input_array)
2838            else:
2839                vedo.logger.error(f"in cmap(), cannot find {on} array at {input_array} ...skip coloring.")
2840                return self
2841
2842        elif utils.is_sequence(input_array):  # if a numpy array is passed
2843            npts = len(input_array)
2844            if npts != n:
2845                vedo.logger.error(f"in cmap(), nr. of input {on} scalars {npts} != {n} ...skip coloring.")
2846                return self
2847            arr = utils.numpy2vtk(input_array, name=name, dtype=float)
2848            data.AddArray(arr)
2849            data.Modified()
2850
2851        elif isinstance(input_array, vtk.vtkArray):  # if a vtkArray is passed
2852            arr = input_array
2853            data.AddArray(arr)
2854            data.Modified()
2855
2856        else:
2857            vedo.logger.error(f"in cmap(), cannot understand input type {type(input_array)}")
2858            raise RuntimeError()
2859
2860        # Now we have array "arr"
2861        array_name = arr.GetName()
2862
2863        if arr.GetNumberOfComponents() == 1:
2864            if vmin is None:
2865                vmin = arr.GetRange()[0]
2866            if vmax is None:
2867                vmax = arr.GetRange()[1]
2868        else:
2869            if vmin is None or vmax is None:
2870                vn = utils.mag(utils.vtk2numpy(arr))
2871            if vmin is None:
2872                vmin = vn.min()
2873            if vmax is None:
2874                vmax = vn.max()
2875
2876        # interpolate alphas if they are not constant
2877        if not utils.is_sequence(alpha):
2878            alpha = [alpha] * n_colors
2879        else:
2880            v = np.linspace(0,1, n_colors, endpoint=True)
2881            xp = np.linspace(0,1, len(alpha), endpoint=True)
2882            alpha = np.interp(v, xp, alpha)
2883
2884        ########################### build the look-up table
2885        if isinstance(input_cmap, vtk.vtkLookupTable):  # vtkLookupTable
2886            lut = input_cmap
2887
2888        elif utils.is_sequence(input_cmap):  # manual sequence of colors
2889            lut = vtk.vtkLookupTable()
2890            if logscale:
2891                lut.SetScaleToLog10()
2892            lut.SetRange(vmin, vmax)
2893            ncols = len(input_cmap)
2894            lut.SetNumberOfTableValues(ncols)
2895
2896            for i, c in enumerate(input_cmap):
2897                r, g, b = colors.get_color(c)
2898                lut.SetTableValue(i, r, g, b, alpha[i])
2899            lut.Build()
2900
2901        else:  # assume string cmap name OR matplotlib.colors.LinearSegmentedColormap
2902            lut = vtk.vtkLookupTable()
2903            if logscale:
2904                lut.SetScaleToLog10()
2905            lut.SetVectorModeToMagnitude()
2906            lut.SetRange(vmin, vmax)
2907            lut.SetNumberOfTableValues(n_colors)
2908            mycols = colors.color_map(range(n_colors), input_cmap, 0, n_colors)
2909            for i, c in enumerate(mycols):
2910                r, g, b = c
2911                lut.SetTableValue(i, r, g, b, alpha[i])
2912            lut.Build()
2913
2914        arr.SetLookupTable(lut)
2915
2916        data.SetActiveScalars(array_name)
2917        # data.SetScalars(arr)  # wrong! it deletes array in position 0, never use SetScalars
2918        # data.SetActiveAttribute(array_name, 0) # boh!
2919
2920        if data.GetScalars():
2921            data.GetScalars().SetLookupTable(lut)
2922            data.GetScalars().Modified()
2923
2924        self._mapper.SetLookupTable(lut)
2925        self._mapper.SetColorModeToMapScalars() # so we dont need to convert uint8 scalars
2926
2927        self._mapper.ScalarVisibilityOn()
2928        self._mapper.SetScalarRange(lut.GetRange())
2929        if on.startswith("point"):
2930            self._mapper.SetScalarModeToUsePointData()
2931        else:
2932            self._mapper.SetScalarModeToUseCellData()
2933        if hasattr(self._mapper, "SetArrayName"):
2934            self._mapper.SetArrayName(array_name)
2935
2936        return self
2937
2938    @deprecated(reason=vedo.colors.red + "Please use property mesh.cellcolors" + vedo.colors.reset)
2939    def cellIndividualColors(self, colorlist):
2940        self.cellcolors = colorlist
2941        return self
2942
2943    @deprecated(reason=vedo.colors.red + "Please use property mesh.cellcolors" + vedo.colors.reset)
2944    def cell_individual_colors(self, colorlist):
2945        self.cellcolors = colorlist
2946        return self
2947
2948    @property
2949    def cellcolors(self):
2950        """
2951        Colorize each cell (face) of a mesh by passing
2952        a 1-to-1 list of colors in format [R,G,B] or [R,G,B,A].
2953        Colors levels and opacities must be in the range [0,255].
2954
2955        A single constant color can also be passed as string or RGBA.
2956
2957        A cell array named "CellsRGBA" is automatically created.
2958
2959        Examples:
2960            - [color_mesh_cells1.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/color_mesh_cells1.py)
2961            - [color_mesh_cells2.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/color_mesh_cells2.py)
2962
2963            ![](https://vedo.embl.es/images/basic/colorMeshCells.png)
2964        """
2965        if "CellsRGBA" not in self.celldata.keys():
2966            lut = self.mapper().GetLookupTable()
2967            vscalars = self._data.GetCellData().GetScalars()
2968            if vscalars is None or lut is None:
2969                arr = np.zeros([self.ncells, 4], dtype=np.uint8)
2970                col = np.array(self.property.GetColor())
2971                col = np.round(col*255).astype(np.uint8)
2972                alf = self.property.GetOpacity()
2973                alf = np.round(alf*255).astype(np.uint8)
2974                arr[:,(0,1,2)] = col
2975                arr[:,3] = alf
2976            else:
2977                cols = lut.MapScalars(vscalars, 0, 0)
2978                arr = utils.vtk2numpy(cols)
2979            self.celldata["CellsRGBA"] = arr
2980        self.celldata.select("CellsRGBA")
2981        return self.celldata["CellsRGBA"]
2982
2983    @cellcolors.setter
2984    def cellcolors(self, value):
2985        if isinstance(value, str):
2986            c = colors.get_color(value)
2987            value = np.array([*c, 1]) * 255
2988            value = np.round(value)
2989
2990        value = np.asarray(value)
2991        n = self.ncells
2992
2993        if value.ndim == 1:
2994            value = np.repeat([value], n, axis=0)
2995
2996        if value.shape[1] == 3:
2997            z = np.zeros((n,1), dtype=np.uint8)
2998            value = np.append(value, z+255, axis=1)
2999
3000        assert n == value.shape[0]
3001
3002        self.celldata["CellsRGBA"] = value.astype(np.uint8)
3003        self.celldata.select("CellsRGBA")
3004
3005
3006    @property
3007    def pointcolors(self):
3008        """
3009        Colorize each point (or vertex of a mesh) by passing
3010        a 1-to-1 list of colors in format [R,G,B] or [R,G,B,A].
3011        Colors levels and opacities must be in the range [0,255].
3012
3013        A single constant color can also be passed as string or RGBA.
3014
3015        A point array named "PointsRGBA" is automatically created.
3016        """
3017        if "PointsRGBA" not in self.pointdata.keys():
3018            lut = self.mapper().GetLookupTable()
3019            vscalars = self._data.GetPointData().GetScalars()
3020            if vscalars is None or lut is None:
3021                arr = np.zeros([self.npoints, 4], dtype=np.uint8)
3022                col = np.array(self.property.GetColor())
3023                col = np.round(col*255).astype(np.uint8)
3024                alf = self.property.GetOpacity()
3025                alf = np.round(alf*255).astype(np.uint8)
3026                arr[:,(0,1,2)] = col
3027                arr[:,3] = alf
3028            else:
3029                cols = lut.MapScalars(vscalars, 0, 0)
3030                arr = utils.vtk2numpy(cols)
3031            self.pointdata["PointsRGBA"] = arr
3032        self.pointdata.select("PointsRGBA")
3033        return self.pointdata["PointsRGBA"]
3034
3035    @pointcolors.setter
3036    def pointcolors(self, value):
3037        if isinstance(value, str):
3038            c = colors.get_color(value)
3039            value = np.array([*c, 1]) * 255
3040            value = np.round(value)
3041
3042        value = np.asarray(value)
3043        n = self.npoints
3044
3045        if value.ndim == 1:
3046            value = np.repeat([value], n, axis=0)
3047
3048        if value.shape[1] == 3:
3049            z = np.zeros((n,1), dtype=np.uint8)
3050            value = np.append(value, z+255, axis=1)
3051
3052        assert n == value.shape[0]
3053
3054        self.pointdata["PointsRGBA"] = value.astype(np.uint8)
3055        self.pointdata.select("PointsRGBA")
3056
3057
3058    @deprecated(reason=vedo.colors.red + "Please use interpolate_data_from()" + vedo.colors.reset)
3059    def interpolateDataFrom(self,
3060            source,
3061            radius=None,
3062            N=None,
3063            kernel="shepard",
3064            exclude=("Normals",),
3065            on="points",
3066            nullStrategy=1,
3067            nullValue=0,
3068        ):
3069        "Deprecated. Please use `interpolate_data_from()`."
3070        return self.interpolate_data_from(
3071            source, radius, N, kernel, exclude, on, nullStrategy, nullValue,
3072        )
3073
3074    def interpolate_data_from(
3075        self,
3076        source,
3077        radius=None,
3078        n=None,
3079        kernel="shepard",
3080        exclude=("Normals",),
3081        on="points",
3082        null_strategy=1,
3083        null_value=0,
3084    ):
3085        """
3086        Interpolate over source to port its data onto the current object using various kernels.
3087
3088        If n (number of closest points to use) is set then radius value is ignored.
3089
3090        Arguments:
3091            kernel : (str)
3092                available kernels are [shepard, gaussian, linear]
3093            null_strategy : (int)
3094                specify a strategy to use when encountering a "null" point
3095                during the interpolation process. Null points occur when the local neighborhood
3096                (of nearby points to interpolate from) is empty.
3097
3098                - Case 0: an output array is created that marks points
3099                  as being valid (=1) or null (invalid =0), and the null_value is set as well
3100                - Case 1: the output data value(s) are set to the provided null_value
3101                - Case 2: simply use the closest point to perform the interpolation.
3102            null_value : (float)
3103                see above.
3104
3105        Examples:
3106            - [interpolateMeshArray.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/interpolateMeshArray.py)
3107
3108                ![](https://vedo.embl.es/images/advanced/interpolateMeshArray.png)
3109        """
3110        if radius is None and not n:
3111            vedo.logger.error("in interpolate_data_from(): please set either radius or n")
3112            raise RuntimeError
3113
3114        if on == "points":
3115            points = source.polydata()
3116        elif on == "cells":
3117            poly2 = vtk.vtkPolyData()
3118            poly2.ShallowCopy(source.polydata())
3119            c2p = vtk.vtkCellDataToPointData()
3120            c2p.SetInputData(poly2)
3121            c2p.Update()
3122            points = c2p.GetOutput()
3123        else:
3124            vedo.logger.error("in interpolate_data_from(), on must be on points or cells")
3125            raise RuntimeError()
3126
3127        locator = vtk.vtkPointLocator()
3128        locator.SetDataSet(points)
3129        locator.BuildLocator()
3130
3131        if kernel.lower() == "shepard":
3132            kern = vtk.vtkShepardKernel()
3133            kern.SetPowerParameter(2)
3134        elif kernel.lower() == "gaussian":
3135            kern = vtk.vtkGaussianKernel()
3136            kern.SetSharpness(2)
3137        elif kernel.lower() == "linear":
3138            kern = vtk.vtkLinearKernel()
3139        else:
3140            vedo.logger.error("available kernels are: [shepard, gaussian, linear]")
3141            raise RuntimeError()
3142
3143        if n:
3144            kern.SetNumberOfPoints(n)
3145            kern.SetKernelFootprintToNClosest()
3146        else:
3147            kern.SetRadius(radius)
3148
3149        interpolator = vtk.vtkPointInterpolator()
3150        interpolator.SetInputData(self.polydata())
3151        interpolator.SetSourceData(points)
3152        interpolator.SetKernel(kern)
3153        interpolator.SetLocator(locator)
3154        interpolator.PassFieldArraysOff()
3155        interpolator.SetNullPointsStrategy(null_strategy)
3156        interpolator.SetNullValue(null_value)
3157        interpolator.SetValidPointsMaskArrayName("ValidPointMask")
3158        for ex in exclude:
3159            interpolator.AddExcludedArray(ex)
3160        interpolator.Update()
3161
3162        if on == "cells":
3163            p2c = vtk.vtkPointDataToCellData()
3164            p2c.SetInputData(interpolator.GetOutput())
3165            p2c.Update()
3166            cpoly = p2c.GetOutput()
3167        else:
3168            cpoly = interpolator.GetOutput()
3169
3170        if self.GetIsIdentity() or cpoly.GetNumberOfPoints() == 0:
3171            self._update(cpoly)
3172        else:
3173            # bring the underlying polydata to where _data is
3174            M = vtk.vtkMatrix4x4()
3175            M.DeepCopy(self.GetMatrix())
3176            M.Invert()
3177            tr = vtk.vtkTransform()
3178            tr.SetMatrix(M)
3179            tf = vtk.vtkTransformPolyDataFilter()
3180            tf.SetTransform(tr)
3181            tf.SetInputData(cpoly)
3182            tf.Update()
3183            self._update(tf.GetOutput())
3184    
3185        self.pipeline = utils.OperationNode(
3186            "interpolate_data_from", parents=[self, source], 
3187        )
3188        return self
3189
3190    def add_gaussian_noise(self, sigma=1):
3191        """
3192        Add gaussian noise to point positions.
3193        An extra array is added named "GaussianNoise" with the shifts.
3194
3195        Arguments:
3196            sigma : (float)
3197                nr. of standard deviations, expressed in percent of the diagonal size of mesh.
3198                Can also be a list [sigma_x, sigma_y, sigma_z].
3199
3200        Examples:
3201            ```python
3202            from vedo import Sphere
3203            Sphere().add_gaussian_noise(1.0).point_size(8).show().close()
3204            ```
3205        """
3206        sz = self.diagonal_size()
3207        pts = self.points()
3208        n = len(pts)
3209        ns = (np.random.randn(n, 3) * sigma) * (sz / 100)
3210        vpts = vtk.vtkPoints()
3211        vpts.SetNumberOfPoints(n)
3212        vpts.SetData(utils.numpy2vtk(pts + ns, dtype=np.float32))
3213        self.inputdata().SetPoints(vpts)
3214        self.inputdata().GetPoints().Modified()
3215        self.pointdata["GaussianNoise"] = -ns
3216        self.pipeline = utils.OperationNode(
3217            "gaussian_noise", 
3218            parents=[self], 
3219            shape="egg",
3220            comment=f"sigma = {sigma}",
3221        )
3222        return self
3223
3224    @deprecated(reason=vedo.colors.red + "Please use closest_point()" + vedo.colors.reset)
3225    def closestPoint(self, pt, N=1, radius=None, returnPointId=False, returnCellId=False):
3226        "Deprecated. Please use `closest_point()`."
3227        return self.closest_point(pt, N, radius, returnPointId, returnCellId)
3228
3229    def closest_point(self, pt, n=1, radius=None, return_point_id=False, return_cell_id=False):
3230        """
3231        Find the closest point(s) on a mesh given from the input point `pt`.
3232
3233        Arguments:
3234            n : (int)
3235                if greater than 1, return a list of n ordered closest points
3236            radius : (float)
3237                if given, get all points within that radius. Then n is ignored.
3238            return_point_id : (bool)
3239                return point ID instead of coordinates
3240            return_cell_id : (bool)
3241                return cell ID in which the closest point sits
3242
3243        Examples:
3244            - [align1.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/align1.py)
3245            - [fitplanes.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/fitplanes.py)
3246            - [quadratic_morphing.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/quadratic_morphing.py)
3247
3248        .. note::
3249            The appropriate tree search locator is built on the fly and cached for speed.
3250
3251            If you want to reset it use `mymesh.point_locator=None`
3252        """
3253        # NB: every time the mesh moves or is warped the locators are set to None
3254        if (n > 1 or radius) or (n == 1 and return_point_id):
3255
3256            poly = None
3257            if not self.point_locator:
3258                poly = self.polydata()
3259                self.point_locator = vtk.vtkStaticPointLocator()
3260                self.point_locator.SetDataSet(poly)
3261                self.point_locator.BuildLocator()
3262
3263            ##########
3264            if radius:
3265                vtklist = vtk.vtkIdList()
3266                self.point_locator.FindPointsWithinRadius(radius, pt, vtklist)
3267            elif n > 1:
3268                vtklist = vtk.vtkIdList()
3269                self.point_locator.FindClosestNPoints(n, pt, vtklist)
3270            else:  # n==1 hence return_point_id==True
3271                ########
3272                return self.point_locator.FindClosestPoint(pt)
3273                ########
3274
3275            if return_point_id:
3276                ########
3277                return utils.vtk2numpy(vtklist)
3278                ########
3279
3280            if not poly:
3281                poly = self.polydata()
3282            trgp = []
3283            for i in range(vtklist.GetNumberOfIds()):
3284                trgp_ = [0, 0, 0]
3285                vi = vtklist.GetId(i)
3286                poly.GetPoints().GetPoint(vi, trgp_)
3287                trgp.append(trgp_)
3288            ########
3289            return np.array(trgp)
3290            ########
3291
3292        else:
3293
3294            if not self.cell_locator:
3295                poly = self.polydata()
3296
3297                # As per Miquel example with limbs the vtkStaticCellLocator doesnt work !!
3298                # https://discourse.vtk.org/t/vtkstaticcelllocator-problem-vtk9-0-3/7854/4
3299                if vedo.vtk_version[0] >= 9 and vedo.vtk_version[0] > 0:
3300                    self.cell_locator = vtk.vtkStaticCellLocator()
3301                else:
3302                    self.cell_locator = vtk.vtkCellLocator()
3303
3304                self.cell_locator.SetDataSet(poly)
3305                self.cell_locator.BuildLocator()
3306
3307            trgp = [0, 0, 0]
3308            cid = vtk.mutable(0)
3309            dist2 = vtk.mutable(0)
3310            subid = vtk.mutable(0)
3311            self.cell_locator.FindClosestPoint(pt, trgp, cid, subid, dist2)
3312
3313            if return_cell_id:
3314                return int(cid)
3315
3316            return np.array(trgp)
3317
3318
3319    def hausdorff_distance(self, points):
3320        """
3321        Compute the Hausdorff distance to the input point set.
3322        Returns a single `float`.
3323
3324        Example:
3325            ```python
3326            from vedo import *
3327            t = np.linspace(0, 2*np.pi, 100)
3328            x = 4/3 * sin(t)**3
3329            y = cos(t) - cos(2*t)/3 - cos(3*t)/6 - cos(4*t)/12
3330            pol1 = Line(np.c_[x,y], closed=True).triangulate()
3331            pol2 = Polygon(nsides=5).pos(2,2)
3332            d12 = pol1.distance_to(pol2)
3333            d21 = pol2.distance_to(pol1)
3334            pol1.lw(0).cmap("viridis")
3335            pol2.lw(0).cmap("viridis")
3336            print("distance d12, d21 :", min(d12), min(d21))
3337            print("hausdorff distance:", pol1.hausdorff_distance(pol2))
3338            print("chamfer distance  :", pol1.chamfer_distance(pol2))
3339            show(pol1, pol2, axes=1)
3340            ```
3341            ![](https://vedo.embl.es/images/feats/heart.png)
3342        """
3343        hp = vtk.vtkHausdorffDistancePointSetFilter()
3344        hp.SetInputData(0, self.polydata())
3345        hp.SetInputData(1, points.polydata())
3346        hp.SetTargetDistanceMethodToPointToCell()
3347        hp.Update()
3348        return hp.GetHausdorffDistance()
3349
3350    def chamfer_distance(self, pcloud):
3351        """
3352        Compute the Chamfer distance to the input point set.
3353        Returns a single `float`.
3354        """
3355        if not pcloud.point_locator:
3356            pcloud.point_locator = vtk.vtkPointLocator()
3357            pcloud.point_locator.SetDataSet(pcloud.polydata())
3358            pcloud.point_locator.BuildLocator()
3359        if not self.point_locator:
3360            self.point_locator = vtk.vtkPointLocator()
3361            self.point_locator.SetDataSet(self.polydata())
3362            self.point_locator.BuildLocator()
3363
3364        ps1 = self.points()
3365        ps2 = pcloud.points()
3366
3367        ids12 = []
3368        for p in ps1:
3369            pid12 = pcloud.point_locator.FindClosestPoint(p)
3370            ids12.append(pid12)
3371        deltav = ps2[ids12] - ps1
3372        da = np.mean(np.linalg.norm(deltav, axis=1))
3373
3374        ids21 = []
3375        for p in ps2:
3376            pid21 = self.point_locator.FindClosestPoint(p)
3377            ids21.append(pid21)
3378        deltav = ps1[ids21] - ps2
3379        db = np.mean(np.linalg.norm(deltav, axis=1))
3380        return (da + db) / 2
3381
3382    def remove_outliers(self, radius, neighbors=5):
3383        """
3384        Remove outliers from a cloud of points within the specified `radius` search.
3385
3386        Arguments:
3387            radius : (float)
3388                Specify the local search radius.
3389            neighbors : (int)
3390                Specify the number of neighbors that a point must have,
3391                within the specified radius, for the point to not be considered isolated.
3392
3393        Examples:
3394            - [clustering.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/clustering.py)
3395
3396                ![](https://vedo.embl.es/images/basic/clustering.png)
3397        """
3398        removal = vtk.vtkRadiusOutlierRemoval()
3399        removal.SetInputData(self.polydata())
3400        removal.SetRadius(radius)
3401        removal.SetNumberOfNeighbors(neighbors)
3402        removal.GenerateOutliersOff()
3403        removal.Update()
3404        inputobj = removal.GetOutput()
3405        if inputobj.GetNumberOfCells() == 0:
3406            carr = vtk.vtkCellArray()
3407            for i in range(inputobj.GetNumberOfPoints()):
3408                carr.InsertNextCell(1)
3409                carr.InsertCellPoint(i)
3410            inputobj.SetVerts(carr)
3411        self._update(inputobj)
3412        self.mapper().ScalarVisibilityOff()
3413        self.pipeline = utils.OperationNode("remove\outliers", parents=[self])
3414        return self
3415
3416    def smooth_mls_1d(self, f=0.2, radius=None):
3417        """
3418        Smooth mesh or points with a `Moving Least Squares` variant.
3419        The point data array "Variances" will contain the residue calculated for each point.
3420        Input mesh's polydata is modified.
3421
3422        Arguments:
3423            f : (float)
3424                smoothing factor - typical range is [0,2].
3425            radius : (float)
3426                radius search in absolute units. If set then `f` is ignored.
3427
3428        Examples:
3429            - [moving_least_squares1D.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/moving_least_squares1D.py)
3430            - [skeletonize.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/skeletonize.py)
3431
3432            ![](https://vedo.embl.es/images/advanced/moving_least_squares1D.png)
3433        """
3434        coords = self.points()
3435        ncoords = len(coords)
3436
3437        if radius:
3438            Ncp = 0
3439        else:
3440            Ncp = int(ncoords * f / 10)
3441            if Ncp < 5:
3442                vedo.logger.warning(f"Please choose a fraction higher than {f}")
3443                Ncp = 5
3444
3445        variances, newline = [], []
3446        for p in coords:
3447            points = self.closest_point(p, n=Ncp, radius=radius)
3448            if len(points) < 4:
3449                continue
3450
3451            points = np.array(points)
3452            pointsmean = points.mean(axis=0)  # plane center
3453            _, dd, vv = np.linalg.svd(points - pointsmean)
3454            newp = np.dot(p - pointsmean, vv[0]) * vv[0] + pointsmean
3455            variances.append(dd[1] + dd[2])
3456            newline.append(newp)
3457
3458        vdata = utils.numpy2vtk(np.array(variances))
3459        vdata.SetName("Variances")
3460        self.inputdata().GetPointData().AddArray(vdata)
3461        self.inputdata().GetPointData().Modified()
3462        self.points(newline)
3463        self.pipeline = utils.OperationNode("smooth_mls_1d", parents=[self])
3464        return self
3465
3466    def smooth_mls_2d(self, f=0.2, radius=None):
3467        """
3468        Smooth mesh or points with a `Moving Least Squares` algorithm variant.
3469        The list `mesh.info['variances']` contains the residue calculated for each point.
3470        When a radius is specified points that are isolated will not be moved and will get
3471        a False entry in array `mesh.info['is_valid']`.
3472
3473        Arguments:
3474            f : (float)
3475                smoothing factor - typical range is [0,2].
3476            radius : (float)
3477                radius search in absolute units. If set then `f` is ignored.
3478
3479        Examples:
3480            - [moving_least_squares2D.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/moving_least_squares2D.py)
3481            - [recosurface.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/recosurface.py)
3482
3483                ![](https://vedo.embl.es/images/advanced/recosurface.png)
3484        """
3485        coords = self.points()
3486        ncoords = len(coords)
3487
3488        if radius:
3489            Ncp = 1
3490        else:
3491            Ncp = int(ncoords * f / 100)
3492            if Ncp < 4:
3493                vedo.logger.error(f"MLS2D: Please choose a fraction higher than {f}")
3494                Ncp = 4
3495
3496        variances, newpts, valid = [], [], []
3497        pb = None
3498        if ncoords > 10000:
3499            pb = utils.ProgressBar(0, ncoords)
3500        for p in coords:
3501            if pb:
3502                pb.print("smoothMLS2D working ...")
3503            pts = self.closest_point(p, n=Ncp, radius=radius)
3504            if len(pts) > 3:
3505                ptsmean = pts.mean(axis=0)  # plane center
3506                _, dd, vv = np.linalg.svd(pts - ptsmean)
3507                cv = np.cross(vv[0], vv[1])
3508                t = (np.dot(cv, ptsmean) - np.dot(cv, p)) / np.dot(cv, cv)
3509                newp = p + cv * t
3510                newpts.append(newp)
3511                variances.append(dd[2])
3512                if radius:
3513                    valid.append(True)
3514            else:
3515                newpts.append(p)
3516                variances.append(0)
3517                if radius:
3518                    valid.append(False)
3519
3520        self.info["variances"] = np.array(variances)
3521        self.info["is_valid"] = np.array(valid)
3522        self.points(newpts)
3523
3524        self.pipeline = utils.OperationNode("smooth_mls_2d", parents=[self])
3525        return self 
3526    
3527    def smooth_lloyd_2d(self, interations=2, bounds=None, options="Qbb Qc Qx"):
3528        """Lloyd relaxation of a 2D pointcloud."""
3529        # Credits: https://hatarilabs.com/ih-en/
3530        # tutorial-to-create-a-geospatial-voronoi-sh-mesh-with-python-scipy-and-geopandas
3531        from scipy.spatial import Voronoi as scipy_voronoi
3532
3533        def _constrain_points(points):
3534            # Update any points that have drifted beyond the boundaries of this space
3535            if bounds is not None:
3536                for point in points:
3537                    if point[0] < bounds[0]: point[0] = bounds[0]
3538                    if point[0] > bounds[1]: point[0] = bounds[1]
3539                    if point[1] < bounds[2]: point[1] = bounds[2]
3540                    if point[1] > bounds[3]: point[1] = bounds[3]
3541            return points
3542
3543        def _find_centroid(vertices):
3544            # The equation for the method used here to find the centroid of a
3545            # 2D polygon is given here: https://en.wikipedia.org/wiki/Centroid#Of_a_polygon
3546            area = 0
3547            centroid_x = 0
3548            centroid_y = 0
3549            for i in range(len(vertices)-1):
3550                step = (vertices[i  , 0] * vertices[i+1, 1]) - \
3551                    (vertices[i+1, 0] * vertices[i  , 1])
3552                centroid_x += (vertices[i, 0] + vertices[i+1, 0]) * step
3553                centroid_y += (vertices[i, 1] + vertices[i+1, 1]) * step
3554                area += step
3555            if area:
3556                centroid_x = (1.0 / (3.0 * area)) * centroid_x
3557                centroid_y = (1.0 / (3.0 * area)) * centroid_y
3558            # prevent centroids from escaping bounding box
3559            return _constrain_points([[centroid_x, centroid_y]])[0]
3560
3561        def _relax(voron):
3562            # Moves each point to the centroid of its cell in the voronoi
3563            # map to "relax" the points (i.e. jitter the points so as
3564            # to spread them out within the space).
3565            centroids = []
3566            for idx in voron.point_region:
3567                # the region is a series of indices into voronoi.vertices
3568                # remove point at infinity, designated by index -1
3569                region = [i for i in voron.regions[idx] if i != -1]
3570                # enclose the polygon
3571                region = region + [region[0]]
3572                verts = voron.vertices[region]
3573                # find the centroid of those vertices
3574                centroids.append(_find_centroid(verts))
3575            return _constrain_points(centroids)
3576
3577        if bounds is None:
3578            bounds = self.bounds()
3579
3580        pts = self.points()[:, (0, 1)]
3581        for i in range(interations):
3582            vor = scipy_voronoi(pts, qhull_options=options)
3583            _constrain_points(vor.vertices)
3584            pts = _relax(vor)
3585        # m = vedo.Mesh([pts, self.faces()]) # not yet working properly
3586        out = Points(pts, c="k")
3587        out.pipeline = utils.OperationNode("smooth_lloyd", parents=[self])
3588        return out
3589
3590    def project_on_plane(self, plane="z", point=None, direction=None):
3591        """
3592        Project the mesh on one of the Cartesian planes.
3593
3594        Arguments:
3595            plane : (str, Plane)
3596                if plane is `str`, plane can be one of ['x', 'y', 'z'],
3597                represents x-plane, y-plane and z-plane, respectively.
3598                Otherwise, plane should be an instance of `vedo.shapes.Plane`.
3599            point : (float, array)
3600                if plane is `str`, point should be a float represents the intercept.
3601                Otherwise, point is the camera point of perspective projection
3602            direction : (array)
3603                direction of oblique projection
3604
3605        Note:
3606            Parameters `point` and `direction` are only used if the given plane
3607            is an instance of `vedo.shapes.Plane`. And one of these two params
3608            should be left as `None` to specify the projection type.
3609
3610        Example:
3611            ```python
3612            s.project_on_plane(plane='z') # project to z-plane
3613            plane = Plane(pos=(4, 8, -4), normal=(-1, 0, 1), s=(5,5))
3614            s.project_on_plane(plane=plane)                       # orthogonal projection
3615            s.project_on_plane(plane=plane, point=(6, 6, 6))      # perspective projection
3616            s.project_on_plane(plane=plane, direction=(1, 2, -1)) # oblique projection
3617            ```
3618
3619        Examples:
3620            - [silhouette2.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/silhouette2.py)
3621
3622                ![](https://vedo.embl.es/images/basic/silhouette2.png)
3623        """
3624        coords = self.points()
3625
3626        if plane == "x":
3627            coords[:, 0] = self.GetOrigin()[0]
3628            intercept = self.xbounds()[0] if point is None else point
3629            self.x(intercept)
3630        elif plane == "y":
3631            coords[:, 1] = self.GetOrigin()[1]
3632            intercept = self.ybounds()[0] if point is None else point
3633            self.y(intercept)
3634        elif plane == "z":
3635            coords[:, 2] = self.GetOrigin()[2]
3636            intercept = self.zbounds()[0] if point is None else point
3637            self.z(intercept)
3638
3639        elif isinstance(plane, vedo.shapes.Plane):
3640            normal = plane.normal / np.linalg.norm(plane.normal)
3641            pl = np.hstack((normal, -np.dot(plane.pos(), normal))).reshape(4, 1)
3642            if direction is None and point is None:
3643                # orthogonal projection
3644                pt = np.hstack((normal, [0])).reshape(4, 1)
3645                # proj_mat = pt.T @ pl * np.eye(4) - pt @ pl.T # python3 only
3646                proj_mat = np.matmul(pt.T, pl) * np.eye(4) - np.matmul(pt, pl.T)
3647
3648            elif direction is None:
3649                # perspective projection
3650                pt = np.hstack((np.array(point), [1])).reshape(4, 1)
3651                # proj_mat = pt.T @ pl * np.eye(4) - pt @ pl.T
3652                proj_mat = np.matmul(pt.T, pl) * np.eye(4) - np.matmul(pt, pl.T)
3653
3654            elif point is None:
3655                # oblique projection
3656                pt = np.hstack((np.array(direction), [0])).reshape(4, 1)
3657                # proj_mat = pt.T @ pl * np.eye(4) - pt @ pl.T
3658                proj_mat = np.matmul(pt.T, pl) * np.eye(4) - np.matmul(pt, pl.T)
3659
3660            coords = np.concatenate([coords, np.ones((coords.shape[:-1] + (1,)))], axis=-1)
3661            # coords = coords @ proj_mat.T
3662            coords = np.matmul(coords, proj_mat.T)
3663            coords = coords[:, :3] / coords[:, 3:]
3664
3665        else:
3666            vedo.logger.error(f"unknown plane {plane}")
3667            raise RuntimeError()
3668
3669        self.alpha(0.1)
3670        self.points(coords)
3671        return self
3672
3673    def warp(self, source, target, sigma=1, mode="3d"):
3674        """
3675        `Thin Plate Spline` transformations describe a nonlinear warp transform defined by a set
3676        of source and target landmarks. Any point on the mesh close to a source landmark will
3677        be moved to a place close to the corresponding target landmark.
3678        The points in between are interpolated smoothly using
3679        Bookstein's Thin Plate Spline algorithm.
3680
3681        Transformation object can be accessed with `mesh.transform`.
3682
3683        Arguments:
3684            sigma : (float)
3685                specify the 'stiffness' of the spline.
3686            mode : (str)
3687                set the basis function to either abs(R) (for 3d) or R2LogR (for 2d meshes)
3688
3689        Examples:
3690            - [interpolate_field.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/interpolate_field.py)
3691            - [warp1.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/warp1.py)
3692            - [warp2.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/warp2.py)
3693
3694                ![](https://vedo.embl.es/images/advanced/warp2.png)
3695        """
3696        parents = [self]
3697        if isinstance(source, Points):
3698            parents.append(source)
3699            source = source.points()
3700        if isinstance(target, Points):
3701            parents.append(target)
3702            target = target.points()
3703
3704        ns = len(source)
3705        ptsou = vtk.vtkPoints()
3706        ptsou.SetNumberOfPoints(ns)
3707        for i in range(ns):
3708            ptsou.SetPoint(i, source[i])
3709
3710        nt = len(target)
3711        if ns != nt:
3712            vedo.logger.error(f"#source {ns} != {nt} #target points")
3713            raise RuntimeError()
3714
3715        pttar = vtk.vtkPoints()
3716        pttar.SetNumberOfPoints(nt)
3717        for i in range(ns):
3718            pttar.SetPoint(i, target[i])
3719
3720        T = vtk.vtkThinPlateSplineTransform()
3721        if mode.lower() == "3d":
3722            T.SetBasisToR()
3723        elif mode.lower() == "2d":
3724            T.SetBasisToR2LogR()
3725        else:
3726            vedo.logger.error(f"unknown mode {mode}")
3727            raise RuntimeError()
3728
3729        T.SetSigma(sigma)
3730        T.SetSourceLandmarks(ptsou)
3731        T.SetTargetLandmarks(pttar)
3732        self.transform = T
3733        self.apply_transform(T, reset=True)
3734
3735        self.pipeline = utils.OperationNode(
3736            "warp", parents=parents, 
3737        )
3738        return self
3739
3740    @deprecated(reason=vedo.colors.red + "Please use cut_with_plane()" + vedo.colors.reset)
3741    def cutWithPlane(self, origin=(0, 0, 0), normal=(1, 0, 0)):
3742        "Deprecated. Please use `cut_with_plane()`."
3743        return self.cut_with_plane(origin, normal)
3744
3745    @deprecated(reason=vedo.colors.red + "Please use cut_with_mesh()" + vedo.colors.reset)
3746    def cutWithMesh(self, mesh, invert=False, keep=False):
3747        "Deprecated. Please use `cut_with_mesh()`."
3748        return self.cut_with_mesh(mesh, invert, keep)
3749
3750    def cut_with_plane(self, origin=(0, 0, 0), normal=(1, 0, 0), invert=False):
3751        """
3752        Cut the mesh with the plane defined by a point and a normal.
3753
3754        Arguments:
3755            origin : (array)
3756                the cutting plane goes through this point
3757            normal : (array)
3758                normal of the cutting plane
3759
3760        Example:
3761            ```python
3762            from vedo import Cube
3763            cube = Cube().cut_with_plane(normal=(1,1,1))
3764            cube.back_color('pink').show()
3765            ```
3766            ![](https://vedo.embl.es/images/feats/cut_with_plane_cube.png)
3767
3768        Examples:
3769            - [trail.py](https://github.com/marcomusy/vedo/blob/master/examples/simulations/trail.py)
3770
3771                ![](https://vedo.embl.es/images/simulations/trail.gif)
3772
3773        Check out also:
3774            `cut_with_box()`, `cut_with_cylinder()`, `cut_with_sphere()`.
3775        """
3776        s = str(normal)
3777        if "x" in s:
3778            normal = (1, 0, 0)
3779            if "-" in s:
3780                normal = -np.array(normal)
3781        elif "y" in s:
3782            normal = (0, 1, 0)
3783            if "-" in s:
3784                normal = -np.array(normal)
3785        elif "z" in s:
3786            normal = (0, 0, 1)
3787            if "-" in s:
3788                normal = -np.array(normal)
3789        plane = vtk.vtkPlane()
3790        plane.SetOrigin(origin)
3791        plane.SetNormal(normal)
3792
3793        clipper = vtk.vtkClipPolyData()
3794        clipper.SetInputData(self.polydata(True))  # must be True
3795        clipper.SetClipFunction(plane)
3796        clipper.GenerateClippedOutputOff()
3797        clipper.GenerateClipScalarsOff()
3798        clipper.SetInsideOut(invert)
3799        clipper.SetValue(0)
3800        clipper.Update()
3801
3802        cpoly = clipper.GetOutput()
3803
3804        if self.GetIsIdentity() or cpoly.GetNumberOfPoints() == 0:
3805            self._update(cpoly)
3806        else:
3807            # bring the underlying polydata to where _data is
3808            M = vtk.vtkMatrix4x4()
3809            M.DeepCopy(self.GetMatrix())
3810            M.Invert()
3811            tr = vtk.vtkTransform()
3812            tr.SetMatrix(M)
3813            tf = vtk.vtkTransformPolyDataFilter()
3814            tf.SetTransform(tr)
3815            tf.SetInputData(cpoly)
3816            tf.Update()
3817            self._update(tf.GetOutput())
3818
3819        self.pipeline = utils.OperationNode("cut_with_plane", parents=[self])
3820        return self
3821
3822    def cut_with_planes(self, origins, normals, invert=False):
3823        """
3824        Cut the mesh with a convex set of planes defined by points and normals.
3825
3826        Arguments:
3827            origins : (array)
3828                each cutting plane goes through this point
3829            normals : (array)
3830                normal of each of the cutting planes
3831
3832        Check out also:
3833            `cut_with_box()`, `cut_with_cylinder()`, `cut_with_sphere()`
3834        """
3835
3836        vpoints = vtk.vtkPoints()
3837        for p in utils.make3d(origins):
3838            vpoints.InsertNextPoint(p)
3839        normals = utils.make3d(normals)
3840
3841        planes = vtk.vtkPlanes()
3842        planes.SetPoints(vpoints)
3843        planes.SetNormals(utils.numpy2vtk(normals, dtype=float))
3844
3845        clipper = vtk.vtkClipPolyData()
3846        clipper.SetInputData(self.polydata(True))  # must be True
3847        clipper.SetInsideOut(invert)
3848        clipper.SetClipFunction(planes)
3849        clipper.GenerateClippedOutputOff()
3850        clipper.GenerateClipScalarsOff()
3851        clipper.SetValue(0)
3852        clipper.Update()
3853
3854        cpoly = clipper.GetOutput()
3855
3856        if self.GetIsIdentity() or cpoly.GetNumberOfPoints() == 0:
3857            self._update(cpoly)
3858        else:
3859            # bring the underlying polydata to where _data is
3860            M = vtk.vtkMatrix4x4()
3861            M.DeepCopy(self.GetMatrix())
3862            M.Invert()
3863            tr = vtk.vtkTransform()
3864            tr.SetMatrix(M)
3865            tf = vtk.vtkTransformPolyDataFilter()
3866            tf.SetTransform(tr)
3867            tf.SetInputData(cpoly)
3868            tf.Update()
3869            self._update(tf.GetOutput())
3870
3871        self.pipeline = utils.OperationNode("cut_with_planes", parents=[self])
3872        return self
3873
3874    def cut_with_box(self, bounds, invert=False):
3875        """
3876        Cut the current mesh with a box or a set of boxes.
3877        This is much faster than `cut_with_mesh()`.
3878
3879        Input `bounds` can be either:
3880        - a Mesh or Points object
3881        - a list of 6 number representing a bounding box `[xmin,xmax, ymin,ymax, zmin,zmax]`
3882        - a list of bounding boxes like the above: `[[xmin1,...], [xmin2,...], ...]`
3883
3884        Example:
3885            ```python
3886            from vedo import Sphere, Cube, show
3887            mesh = Sphere(r=1, res=50)
3888            box  = Cube(side=1.5).wireframe()
3889            mesh.cut_with_box(box)
3890            show(mesh, box, axes=1)
3891            ```
3892            ![](https://vedo.embl.es/images/feats/cut_with_box_cube.png)
3893
3894        Check out also:
3895            `cut_with_line()`, `cut_with_plane()`, `cut_with_cylinder()`
3896        """
3897        if isinstance(bounds, Points):
3898            bounds = bounds.bounds()
3899
3900        box = vtk.vtkBox()
3901        if utils.is_sequence(bounds[0]):
3902            for bs in bounds:
3903                box.AddBounds(bs)
3904        else:
3905            box.SetBounds(bounds)
3906
3907        clipper = vtk.vtkClipPolyData()
3908        clipper.SetInputData(self.polydata(True))  # must be True
3909        clipper.SetClipFunction(box)
3910        clipper.SetInsideOut(not invert)
3911        clipper.GenerateClippedOutputOff()
3912        clipper.GenerateClipScalarsOff()
3913        clipper.SetValue(0)
3914        clipper.Update()
3915        cpoly = clipper.GetOutput()
3916
3917        if self.GetIsIdentity() or cpoly.GetNumberOfPoints() == 0:
3918            self._update(cpoly)
3919        else:
3920            # bring the underlying polydata to where _data is
3921            M = vtk.vtkMatrix4x4()
3922            M.DeepCopy(self.GetMatrix())
3923            M.Invert()
3924            tr = vtk.vtkTransform()
3925            tr.SetMatrix(M)
3926            tf = vtk.vtkTransformPolyDataFilter()
3927            tf.SetTransform(tr)
3928            tf.SetInputData(cpoly)
3929            tf.Update()
3930            self._update(tf.GetOutput())
3931
3932        self.pipeline = utils.OperationNode("cut_with_box", parents=[self])
3933        return self
3934
3935    def cut_with_line(self, points, invert=False, closed=True):
3936        """
3937        Cut the current mesh with a line vertically in the z-axis direction like a cookie cutter.
3938        The polyline is defined by a set of points (z-coordinates are ignored).
3939        This is much faster than `cut_with_mesh()`.
3940
3941        Check out also:
3942            `cut_with_box()`, `cut_with_plane()`, `cut_with_sphere()`
3943        """
3944        pplane = vtk.vtkPolyPlane()
3945        if isinstance(points, Points):
3946            points = points.points().tolist()
3947
3948        if closed:
3949            if isinstance(points, np.ndarray):
3950                points = points.tolist()
3951            points.append(points[0])
3952
3953        vpoints = vtk.vtkPoints()
3954        for p in points:
3955            if len(p) == 2:
3956                p = [p[0], p[1], 0.0]
3957            vpoints.InsertNextPoint(p)
3958
3959        n = len(points)
3960        polyline = vtk.vtkPolyLine()
3961        polyline.Initialize(n, vpoints)
3962        polyline.GetPointIds().SetNumberOfIds(n)
3963        for i in range(n):
3964            polyline.GetPointIds().SetId(i, i)
3965        pplane.SetPolyLine(polyline)
3966
3967        clipper = vtk.vtkClipPolyData()
3968        clipper.SetInputData(self.polydata(True))  # must be True
3969        clipper.SetClipFunction(pplane)
3970        clipper.SetInsideOut(invert)
3971        clipper.GenerateClippedOutputOff()
3972        clipper.GenerateClipScalarsOff()
3973        clipper.SetValue(0)
3974        clipper.Update()
3975        cpoly = clipper.GetOutput()
3976
3977        if self.GetIsIdentity() or cpoly.GetNumberOfPoints() == 0:
3978            self._update(cpoly)
3979        else:
3980            # bring the underlying polydata to where _data is
3981            M = vtk.vtkMatrix4x4()
3982            M.DeepCopy(self.GetMatrix())
3983            M.Invert()
3984            tr = vtk.vtkTransform()
3985            tr.SetMatrix(M)
3986            tf = vtk.vtkTransformPolyDataFilter()
3987            tf.SetTransform(tr)
3988            tf.SetInputData(cpoly)
3989            tf.Update()
3990            self._update(tf.GetOutput())
3991
3992        self.pipeline = utils.OperationNode("cut_with_line", parents=[self])
3993        return self
3994
3995    def cut_with_cylinder(self, center=(0, 0, 0), axis=(0, 0, 1), r=1, invert=False):
3996        """
3997        Cut the current mesh with an infinite cylinder.
3998        This is much faster than `cut_with_mesh()`.
3999
4000        Arguments:
4001            center : (array)
4002                the center of the cylinder
4003            normal : (array)
4004                direction of the cylinder axis
4005            r : (float)
4006                radius of the cylinder
4007
4008        Example:
4009            ```python
4010            from vedo import Disc, show
4011            disc = Disc(r1=1, r2=1.2)
4012            mesh = disc.extrude(3, res=50).linewidth(1)
4013            mesh.cut_with_cylinder([0,0,2], r=0.4, axis='y', invert=True)
4014            show(mesh, axes=1)
4015            ```
4016            ![](https://vedo.embl.es/images/feats/cut_with_cylinder.png)
4017
4018        Examples:
4019            - [optics_main1.py](https://github.com/marcomusy/vedo/blob/master/examples/simulations/optics_main1.py)
4020
4021        Check out also:
4022            `cut_with_box()`, `cut_with_plane()`, `cut_with_sphere()`
4023        """
4024        s = str(axis)
4025        if "x" in s:
4026            axis = (1, 0, 0)
4027        elif "y" in s:
4028            axis = (0, 1, 0)
4029        elif "z" in s:
4030            axis = (0, 0, 1)
4031        cyl = vtk.vtkCylinder()
4032        cyl.SetCenter(center)
4033        cyl.SetAxis(axis[0], axis[1], axis[2])
4034        cyl.SetRadius(r)
4035
4036        clipper = vtk.vtkClipPolyData()
4037        clipper.SetInputData(self.polydata(True))  # must be True
4038        clipper.SetClipFunction(cyl)
4039        clipper.SetInsideOut(not invert)
4040        clipper.GenerateClippedOutputOff()
4041        clipper.GenerateClipScalarsOff()
4042        clipper.SetValue(0)
4043        clipper.Update()
4044        cpoly = clipper.GetOutput()
4045
4046        if self.GetIsIdentity() or cpoly.GetNumberOfPoints() == 0:
4047            self._update(cpoly)
4048        else:
4049            # bring the underlying polydata to where _data is
4050            M = vtk.vtkMatrix4x4()
4051            M.DeepCopy(self.GetMatrix())
4052            M.Invert()
4053            tr = vtk.vtkTransform()
4054            tr.SetMatrix(M)
4055            tf = vtk.vtkTransformPolyDataFilter()
4056            tf.SetTransform(tr)
4057            tf.SetInputData(cpoly)
4058            tf.Update()
4059            self._update(tf.GetOutput())
4060
4061        self.pipeline = utils.OperationNode("cut_with_cylinder", parents=[self])
4062        return self
4063
4064    def cut_with_sphere(self, center=(0, 0, 0), r=1, invert=False):
4065        """
4066        Cut the current mesh with an sphere.
4067        This is much faster than `cut_with_mesh()`.
4068
4069        Arguments:
4070            center : (array)
4071                the center of the sphere
4072            r : (float)
4073                radius of the sphere
4074
4075        Example:
4076            ```python
4077            from vedo import Disc, show
4078            disc = Disc(r1=1, r2=1.2)
4079            mesh = disc.extrude(3, res=50).linewidth(1)
4080            mesh.cut_with_sphere([1,-0.7,2], r=1.5, invert=True)
4081            show(mesh, axes=1)
4082            ```
4083            ![](https://vedo.embl.es/images/feats/cut_with_sphere.png)
4084
4085        Check out also:
4086            `cut_with_box()`, `cut_with_plane()`, `cut_with_cylinder()`
4087        """
4088        sph = vtk.vtkSphere()
4089        sph.SetCenter(center)
4090        sph.SetRadius(r)
4091
4092        clipper = vtk.vtkClipPolyData()
4093        clipper.SetInputData(self.polydata(True))  # must be True
4094        clipper.SetClipFunction(sph)
4095        clipper.SetInsideOut(not invert)
4096        clipper.GenerateClippedOutputOff()
4097        clipper.GenerateClipScalarsOff()
4098        clipper.SetValue(0)
4099        clipper.Update()
4100        cpoly = clipper.GetOutput()
4101
4102        if self.GetIsIdentity() or cpoly.GetNumberOfPoints() == 0:
4103            self._update(cpoly)
4104        else:
4105            # bring the underlying polydata to where _data is
4106            M = vtk.vtkMatrix4x4()
4107            M.DeepCopy(self.GetMatrix())
4108            M.Invert()
4109            tr = vtk.vtkTransform()
4110            tr.SetMatrix(M)
4111            tf = vtk.vtkTransformPolyDataFilter()
4112            tf.SetTransform(tr)
4113            tf.SetInputData(cpoly)
4114            tf.Update()
4115            self._update(tf.GetOutput())
4116
4117        self.pipeline = utils.OperationNode("cut_with_sphere", parents=[self])
4118        return self
4119
4120    def cut_with_mesh(self, mesh, invert=False, keep=False):
4121        """
4122        Cut an `Mesh` mesh with another `Mesh`.
4123
4124        Use `invert` to invert the selection.
4125
4126        Use `keep` to keep the cutoff part, in this case an `Assembly` is returned:
4127        the "cut" object and the "discarded" part of the original object.
4128        You can access both via `assembly.unpack()` method.
4129
4130        Example:
4131        ```python
4132        from vedo import *
4133        arr = np.random.randn(100000, 3)/2
4134        pts = Points(arr).c('red3').pos(5,0,0)
4135        cube = Cube().pos(4,0.5,0)
4136        assem = pts.cut_with_mesh(cube, keep=True)
4137        show(assem.unpack(), axes=1).close()
4138        ```        
4139        ![](https://vedo.embl.es/images/feats/cut_with_mesh.png)
4140
4141       Check out also:
4142            `cut_with_box()`, `cut_with_plane()`, `cut_with_cylinder()`
4143       """
4144        polymesh = mesh.polydata()
4145        poly = self.polydata()
4146
4147        # Create an array to hold distance information
4148        signed_distances = vtk.vtkFloatArray()
4149        signed_distances.SetNumberOfComponents(1)
4150        signed_distances.SetName("SignedDistances")
4151
4152        # implicit function that will be used to slice the mesh
4153        ippd = vtk.vtkImplicitPolyDataDistance()
4154        ippd.SetInput(polymesh)
4155
4156        # Evaluate the signed distance function at all of the grid points
4157        for pointId in range(poly.GetNumberOfPoints()):
4158            p = poly.GetPoint(pointId)
4159            signed_distance = ippd.EvaluateFunction(p)
4160            signed_distances.InsertNextValue(signed_distance)
4161
4162        currentscals = poly.GetPointData().GetScalars()
4163        if currentscals:
4164            currentscals = currentscals.GetName()
4165
4166        poly.GetPointData().AddArray(signed_distances)
4167        poly.GetPointData().SetActiveScalars("SignedDistances")
4168
4169        clipper = vtk.vtkClipPolyData()
4170        clipper.SetInputData(poly)
4171        clipper.SetInsideOut(not invert)
4172        clipper.SetGenerateClippedOutput(keep)
4173        clipper.SetValue(0.0)
4174        clipper.Update()
4175        cpoly = clipper.GetOutput()
4176        if keep:
4177            kpoly = clipper.GetOutput(1)
4178
4179        vis = False
4180        if currentscals:
4181            cpoly.GetPointData().SetActiveScalars(currentscals)
4182            vis = self.mapper().GetScalarVisibility()
4183
4184        if self.GetIsIdentity() or cpoly.GetNumberOfPoints() == 0:
4185            self._update(cpoly)
4186        else:
4187            # bring the underlying polydata to where _data is
4188            M = vtk.vtkMatrix4x4()
4189            M.DeepCopy(self.GetMatrix())
4190            M.Invert()
4191            tr = vtk.vtkTransform()
4192            tr.SetMatrix(M)
4193            tf = vtk.vtkTransformPolyDataFilter()
4194            tf.SetTransform(tr)
4195            tf.SetInputData(cpoly)
4196            tf.Update()
4197            self._update(tf.GetOutput())
4198
4199        self.pointdata.remove("SignedDistances")
4200        self.mapper().SetScalarVisibility(vis)
4201        if keep:
4202            if isinstance(self, vedo.Mesh):
4203                cutoff = vedo.Mesh(kpoly)
4204            else:
4205                cutoff = vedo.Points(kpoly)
4206            cutoff.property = vtk.vtkProperty()
4207            cutoff.property.DeepCopy(self.property)
4208            cutoff.SetProperty(cutoff.property)
4209            cutoff.c("k5").alpha(0.2)
4210            return vedo.Assembly([self, cutoff])
4211
4212        self.pipeline = utils.OperationNode("cut_with_mesh", parents=[self, mesh])
4213        return self
4214
4215    def cut_with_point_loop(
4216        self,
4217        points,
4218        invert=False,
4219        on="points",
4220        include_boundary=False,
4221    ):
4222        """
4223        Cut an `Mesh` object with a set of points forming a closed loop.
4224
4225        Arguments:
4226            invert : (bool)
4227                invert selection (inside-out)
4228            on : (str)
4229                if 'cells' will extract the whole cells lying inside (or outside) the point loop
4230            include_boundary : (bool)
4231                include cells lying exactly on the boundary line. Only relevant on 'cells' mode
4232
4233        Examples:
4234            - [cut_with_points1.py](https://github.com/marcomusy/vedo/blob/master/examples/advanced/cut_with_points1.py)
4235            
4236                ![](https://vedo.embl.es/images/advanced/cutWithPoints1.png)
4237            
4238            - [cut_with_points2.py](https://github.com/marcomusy/vedo/blob/master/examples/advanced/cut_with_points2.py)
4239        
4240                ![](https://vedo.embl.es/images/advanced/cutWithPoints2.png)
4241        """
4242        if isinstance(points, Points):
4243            parents = [points]
4244            vpts = points.polydata().GetPoints()
4245            points = points.points()
4246        else:
4247            parents = [self]
4248            vpts = vtk.vtkPoints()
4249            points = utils.make3d(points)
4250            for p in points:
4251                vpts.InsertNextPoint(p)
4252
4253        if "cell" in on:
4254            ippd = vtk.vtkImplicitSelectionLoop()
4255            ippd.SetLoop(vpts)
4256            ippd.AutomaticNormalGenerationOn()
4257            clipper = vtk.vtkExtractPolyDataGeometry()
4258            clipper.SetInputData(self.polydata())
4259            clipper.SetImplicitFunction(ippd)
4260            clipper.SetExtractInside(not invert)
4261            clipper.SetExtractBoundaryCells(include_boundary)
4262        else:
4263            spol = vtk.vtkSelectPolyData()
4264            spol.SetLoop(vpts)
4265            spol.GenerateSelectionScalarsOn()
4266            spol.GenerateUnselectedOutputOff()
4267            spol.SetInputData(self.polydata())
4268            spol.Update()
4269            clipper = vtk.vtkClipPolyData()
4270            clipper.SetInputData(spol.GetOutput())
4271            clipper.SetInsideOut(not invert)
4272            clipper.SetValue(0.0)
4273        clipper.Update()
4274        cpoly = clipper.GetOutput()
4275
4276        if self.GetIsIdentity() or cpoly.GetNumberOfPoints() == 0:
4277            self._update(cpoly)
4278        else:
4279            # bring the underlying polydata to where _data is
4280            M = vtk.vtkMatrix4x4()
4281            M.DeepCopy(self.GetMatrix())
4282            M.Invert()
4283            tr = vtk.vtkTransform()
4284            tr.SetMatrix(M)
4285            tf = vtk.vtkTransformPolyDataFilter()
4286            tf.SetTransform(tr)
4287            tf.SetInputData(clipper.GetOutput())
4288            tf.Update()
4289            self._update(tf.GetOutput())
4290        
4291        self.pipeline = utils.OperationNode(
4292            "cut_with_pointloop", parents=parents)
4293        return self
4294
4295    def cut_with_scalars(self, value, name="", invert=False):
4296        """
4297        Cut a mesh or point cloud with some input scalar point-data.
4298
4299        Arguments:
4300            value : (float)
4301                cutting value
4302            name : (str)
4303                array name of the scalars to be used
4304            invert : (bool)
4305                flip selection
4306
4307        Example:
4308            ```python
4309            from vedo import *
4310            s = Sphere().lw(1)
4311            pts = s.points()
4312            scalars = np.sin(3*pts[:,2]) + pts[:,0]
4313            s.pointdata["somevalues"] = scalars
4314            s.cut_with_scalars(0.3)
4315            s.cmap("Spectral", "somevalues").add_scalarbar()
4316            s.show(axes=1).close()
4317            ```
4318            ![](https://vedo.embl.es/images/feats/cut_with_scalars.png)
4319        """
4320        if name:
4321            self.pointdata.select(name)
4322        clipper = vtk.vtkClipPolyData()
4323        clipper.SetInputData(self._data)
4324        clipper.SetValue(value)
4325        clipper.GenerateClippedOutputOff()
4326        clipper.SetInsideOut(not invert)
4327        clipper.Update()
4328        self._update(clipper.GetOutput())
4329
4330        self.pipeline = utils.OperationNode(
4331            "cut_with_scalars", parents=[self],
4332        )
4333        return self
4334
4335    def implicit_modeller(self, distance=0.05, res=(50,50,50), bounds=(), maxdist=None):
4336        """Find the surface which sits at the specified distance from the input one."""
4337        if not bounds:
4338            bounds = self.bounds()
4339
4340        if not maxdist:
4341            maxdist = self.diagonal_size() / 2
4342
4343        imp = vtk.vtkImplicitModeller()
4344        imp.SetInputData(self.polydata())
4345        imp.SetSampleDimensions(res)
4346        imp.SetMaximumDistance(maxdist)
4347        imp.SetModelBounds(bounds)
4348        contour = vtk.vtkContourFilter()
4349        contour.SetInputConnection(imp.GetOutputPort())
4350        contour.SetValue(0, distance)
4351        contour.Update()
4352        poly = contour.GetOutput()
4353        out = vedo.Mesh(poly, c="lb")
4354
4355        out.pipeline = utils.OperationNode(
4356            "implicit_modeller", parents=[self],
4357        )        
4358        return out
4359
4360
4361    @deprecated(reason=vedo.colors.red + "Please use generate_mesh()" + vedo.colors.reset)
4362    def tomesh(self, *args, **kwargs):
4363        """Deprecated, please use `generate_mesh()`."""
4364        return self.generate_mesh(*args, **kwargs)
4365
4366    def generate_mesh(
4367        self,
4368        line_resolution=None,
4369        mesh_resolution=None,
4370        smooth=0,
4371        jitter=0.001,
4372        grid=None,
4373        quads=False,
4374        invert=False,
4375    ):
4376        """
4377        Generate a polygonal Mesh from a closed contour line.
4378        If line is not closed it will be closed with a straight segment.
4379
4380        Arguments:
4381            line_resolution : (int)
4382                resolution of the contour line. The default is None, in this case
4383                the contour is not resampled.
4384            mesh_resolution : (int)
4385                resolution of the internal triangles not touching the boundary.
4386                The default is None.
4387            smooth : (float)
4388                smoothing of the contour before meshing.
4389            jitter : (float)
4390                add a small noise to the internal points.
4391            grid : (Grid)
4392                manually pass a Grid object. The default is True.
4393            quads : (bool)
4394                generate a mesh of quads instead of triangles.
4395            invert : (bool)
4396                flip the line orientation. The default is False.
4397
4398        Examples:
4399            - [line2mesh_tri.py](https://github.com/marcomusy/vedo/blob/master/examples/advanced/line2mesh_tri.py)
4400
4401                ![](https://vedo.embl.es/images/advanced/line2mesh_tri.jpg)
4402
4403            - [line2mesh_quads.py](https://github.com/marcomusy/vedo/blob/master/examples/advanced/line2mesh_quads.py)
4404        
4405                ![](https://vedo.embl.es/images/advanced/line2mesh_quads.png)
4406        """
4407        if line_resolution is None:
4408            contour = vedo.shapes.Line(self.points())
4409        else:
4410            contour = vedo.shapes.Spline(self.points(), smooth=smooth, res=line_resolution)
4411        contour.clean()
4412
4413        length = contour.length()
4414        density = length / contour.npoints
4415        vedo.logger.debug(f"tomesh():\n\tline length = {length}")
4416        vedo.logger.debug(f"\tdensity = {density} length/pt_separation")
4417
4418        x0, x1 = contour.xbounds()
4419        y0, y1 = contour.ybounds()
4420
4421        if grid is None:
4422            if mesh_resolution is None:
4423                resx = int((x1 - x0) / density + 0.5)
4424                resy = int((y1 - y0) / density + 0.5)
4425                vedo.logger.debug(f"tmesh_resolution = {[resx, resy]}")
4426            else:
4427                if utils.is_sequence(mesh_resolution):
4428                    resx, resy = mesh_resolution
4429                else:
4430                    resx, resy = mesh_resolution, mesh_resolution
4431            grid = vedo.shapes.Grid(
4432                [(x0 + x1) / 2, (y0 + y1) / 2, 0],
4433                s=((x1 - x0) * 1.025, (y1 - y0) * 1.025),
4434                res=(resx, resy),
4435            )
4436        else:
4437            grid = grid.clone()
4438
4439        cpts = contour.points()
4440
4441        # make sure it's closed
4442        p0, p1 = cpts[0], cpts[-1]
4443        nj = max(2, int(utils.mag(p1 - p0) / density + 0.5))
4444        joinline = vedo.shapes.Line(p1, p0, res=nj)
4445        contour = vedo.merge(contour, joinline).subsample(0.0001)
4446
4447        ####################################### quads
4448        if quads:
4449            cmesh = grid.clone().cut_with_point_loop(contour, on="cells", invert=invert)
4450            cmesh.wireframe(False).lw(0.5)
4451            cmesh.pipeline = utils.OperationNode(
4452                "generate_mesh", parents=[self, contour], 
4453                comment=f"#quads {cmesh._data.GetNumberOfCells()}",
4454            )
4455            return cmesh
4456        #############################################
4457
4458        grid_tmp = grid.points()
4459
4460        if jitter:
4461            np.random.seed(0)
4462            sigma = 1.0 / np.sqrt(grid.npoints) * grid.diagonal_size() * jitter
4463            vedo.logger.debug(f"\tsigma jittering = {sigma}")
4464            grid_tmp += np.random.rand(grid.npoints, 3) * sigma
4465            grid_tmp[:, 2] = 0.0
4466
4467        todel = []
4468        density /= np.sqrt(3)
4469        vgrid_tmp = Points(grid_tmp)
4470
4471        for p in contour.points():
4472            out = vgrid_tmp.closest_point(p, radius=density, return_point_id=True)
4473            todel += out.tolist()
4474        # cpoints = contour.points()
4475        # for i, p in enumerate(cpoints):
4476        #     if i:
4477        #         den = utils.mag(p-cpoints[i-1])/1.732
4478        #     else:
4479        #         den = density
4480        #     todel += vgrid_tmp.closest_point(p, radius=den, return_point_id=True)
4481