vedo.pointcloud

Submodule to work with point clouds

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