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