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