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