vedo.mesh
Submodule to work with polygonal meshes
1#!/usr/bin/env python3 2# -*- coding: utf-8 -*- 3import os 4import numpy as np 5from deprecated import deprecated 6 7try: 8 import vedo.vtkclasses as vtk 9except ImportError: 10 import vtkmodules.all as vtk 11 12import vedo 13from vedo.colors import color_map 14from vedo.colors import get_color 15from vedo.pointcloud import Points 16from vedo.utils import buildPolyData, is_sequence, mag, mag2, precision 17from vedo.utils import numpy2vtk, vtk2numpy, OperationNode 18 19__docformat__ = "google" 20 21__doc__ = """ 22Submodule to work with polygonal meshes 23 24 25""" 26 27__all__ = ["Mesh"] 28 29 30#################################################### 31class Mesh(Points): 32 """ 33 Build an instance of object `Mesh` derived from `vedo.PointCloud`. 34 """ 35 def __init__( 36 self, 37 inputobj=None, 38 c=None, 39 alpha=1, 40 ): 41 """ 42 Input can be a list of vertices and their connectivity (faces of the polygonal mesh), 43 or directly a `vtkPolydata` object. 44 For point clouds - e.i. no faces - just substitute the `faces` list with `None`. 45 46 Example: 47 `Mesh( [ [[x1,y1,z1],[x2,y2,z2], ...], [[0,1,2], [1,2,3], ...] ] )` 48 49 Arguments: 50 c : (color) 51 color in RGB format, hex, symbol or name 52 alpha : (float) 53 mesh opacity [0,1] 54 55 Examples: 56 - [buildmesh.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/buildmesh.py) 57 (and many others!) 58 59  60 """ 61 Points.__init__(self) 62 63 self.line_locator = None 64 65 self._mapper.SetInterpolateScalarsBeforeMapping( 66 vedo.settings.interpolate_scalars_before_mapping 67 ) 68 69 if vedo.settings.use_polygon_offset: 70 self._mapper.SetResolveCoincidentTopologyToPolygonOffset() 71 pof, pou = ( 72 vedo.settings.polygon_offset_factor, 73 vedo.settings.polygon_offset_units, 74 ) 75 self._mapper.SetResolveCoincidentTopologyPolygonOffsetParameters(pof, pou) 76 77 inputtype = str(type(inputobj)) 78 79 if inputobj is None: 80 pass 81 82 elif isinstance(inputobj, (Mesh, vtk.vtkActor)): 83 polyCopy = vtk.vtkPolyData() 84 polyCopy.DeepCopy(inputobj.GetMapper().GetInput()) 85 self._data = polyCopy 86 self._mapper.SetInputData(polyCopy) 87 self._mapper.SetScalarVisibility(inputobj.GetMapper().GetScalarVisibility()) 88 pr = vtk.vtkProperty() 89 pr.DeepCopy(inputobj.GetProperty()) 90 self.SetProperty(pr) 91 self.property = pr 92 93 elif isinstance(inputobj, vtk.vtkPolyData): 94 if inputobj.GetNumberOfCells() == 0: 95 carr = vtk.vtkCellArray() 96 for i in range(inputobj.GetNumberOfPoints()): 97 carr.InsertNextCell(1) 98 carr.InsertCellPoint(i) 99 inputobj.SetVerts(carr) 100 self._data = inputobj # cache vtkPolyData and mapper for speed 101 102 elif isinstance(inputobj, (vtk.vtkStructuredGrid, vtk.vtkRectilinearGrid)): 103 gf = vtk.vtkGeometryFilter() 104 gf.SetInputData(inputobj) 105 gf.Update() 106 self._data = gf.GetOutput() 107 108 elif "trimesh" in inputtype: 109 tact = vedo.utils.trimesh2vedo(inputobj) 110 self._data = tact.polydata() 111 112 elif "meshio" in inputtype: 113 if len(inputobj.cells): 114 mcells = [] 115 for cellblock in inputobj.cells: 116 if cellblock.type in ("triangle", "quad"): 117 mcells += cellblock.data.tolist() 118 self._data = buildPolyData(inputobj.points, mcells) 119 else: 120 self._data = buildPolyData(inputobj.points, None) 121 # add arrays: 122 try: 123 if len(inputobj.point_data): 124 for k in inputobj.point_data.keys(): 125 vdata = numpy2vtk(inputobj.point_data[k]) 126 vdata.SetName(str(k)) 127 self._data.GetPointData().AddArray(vdata) 128 except AssertionError: 129 print("Could not add meshio point data, skip.") 130 # try: 131 # if len(inputobj.cell_data): 132 # for k in inputobj.cell_data.keys(): 133 # #print(inputobj.cell_data) 134 # exit() 135 # vdata = numpy2vtk(inputobj.cell_data[k]) 136 # vdata.SetName(str(k)) 137 # self._data.GetCellData().AddArray(vdata) 138 # except AssertionError: 139 # print("Could not add meshio cell data, skip.") 140 141 elif "meshlab" in inputtype: 142 self._data = vedo.utils.meshlab2vedo(inputobj) 143 144 elif is_sequence(inputobj): 145 ninp = len(inputobj) 146 if ninp == 0: 147 self._data = vtk.vtkPolyData() 148 elif ninp == 2: # assume [vertices, faces] 149 self._data = buildPolyData(inputobj[0], inputobj[1]) 150 else: # assume [vertices] or vertices 151 self._data = buildPolyData(inputobj, None) 152 153 elif hasattr(inputobj, "GetOutput"): # passing vtk object 154 if hasattr(inputobj, "Update"): 155 inputobj.Update() 156 if isinstance(inputobj.GetOutput(), vtk.vtkPolyData): 157 self._data = inputobj.GetOutput() 158 else: 159 gf = vtk.vtkGeometryFilter() 160 gf.SetInputData(inputobj.GetOutput()) 161 gf.Update() 162 self._data = gf.GetOutput() 163 164 elif isinstance(inputobj, str): 165 dataset = vedo.io.load(inputobj) 166 self.filename = inputobj 167 if "TetMesh" in str(type(dataset)): 168 self._data = dataset.tomesh().polydata(False) 169 else: 170 self._data = dataset.polydata(False) 171 172 else: 173 try: 174 gf = vtk.vtkGeometryFilter() 175 gf.SetInputData(inputobj) 176 gf.Update() 177 self._data = gf.GetOutput() 178 except: 179 vedo.logger.error(f"cannot build mesh from type {inputtype}") 180 raise RuntimeError() 181 182 self._mapper.SetInputData(self._data) 183 184 self.property = self.GetProperty() 185 self.property.SetInterpolationToPhong() 186 187 # set the color by c or by scalar 188 if self._data: 189 190 arrexists = False 191 192 if c is None: 193 ptdata = self._data.GetPointData() 194 cldata = self._data.GetCellData() 195 exclude = ["normals", "tcoord"] 196 197 if cldata.GetNumberOfArrays(): 198 for i in range(cldata.GetNumberOfArrays()): 199 iarr = cldata.GetArray(i) 200 if iarr: 201 icname = iarr.GetName() 202 if icname and all(s not in icname.lower() for s in exclude): 203 cldata.SetActiveScalars(icname) 204 self._mapper.ScalarVisibilityOn() 205 self._mapper.SetScalarModeToUseCellData() 206 self._mapper.SetScalarRange(iarr.GetRange()) 207 arrexists = True 208 break # stop at first good one 209 210 # point come after so it has priority 211 if ptdata.GetNumberOfArrays(): 212 for i in range(ptdata.GetNumberOfArrays()): 213 iarr = ptdata.GetArray(i) 214 if iarr: 215 ipname = iarr.GetName() 216 if ipname and all(s not in ipname.lower() for s in exclude): 217 ptdata.SetActiveScalars(ipname) 218 self._mapper.ScalarVisibilityOn() 219 self._mapper.SetScalarModeToUsePointData() 220 self._mapper.SetScalarRange(iarr.GetRange()) 221 arrexists = True 222 break # stop at first good one 223 224 if not arrexists: 225 if c is None: 226 c = "gold" 227 c = get_color(c) 228 elif isinstance(c, float) and c <= 1: 229 c = color_map(c, "rainbow", 0, 1) 230 else: 231 c = get_color(c) 232 self.property.SetColor(c) 233 self.property.SetAmbient(0.1) 234 self.property.SetDiffuse(1) 235 self.property.SetSpecular(0.05) 236 self.property.SetSpecularPower(5) 237 self._mapper.ScalarVisibilityOff() 238 239 if alpha is not None: 240 self.property.SetOpacity(alpha) 241 242 n = self._data.GetNumberOfPoints() 243 self.pipeline = OperationNode(self, comment=f"#pts {n}") 244 245 def _repr_html_(self): 246 """ 247 HTML representation of the Mesh object for Jupyter Notebooks. 248 249 Returns: 250 HTML text with the image and some properties. 251 """ 252 import io 253 import base64 254 from PIL import Image 255 256 library_name = "vedo.mesh.Mesh" 257 help_url = "https://vedo.embl.es/docs/vedo/mesh.html" 258 259 arr = self.thumbnail() 260 im = Image.fromarray(arr) 261 buffered = io.BytesIO() 262 im.save(buffered, format="PNG", quality=100) 263 encoded = base64.b64encode(buffered.getvalue()).decode("utf-8") 264 url = "data:image/png;base64," + encoded 265 image = f"<img src='{url}'></img>" 266 267 # mesh statisitics 268 bounds = "<br/>".join( 269 [ 270 precision(min_x,4) + " ... " + precision(max_x,4) 271 for min_x, max_x in zip(self.bounds()[::2], self.bounds()[1::2]) 272 ] 273 ) 274 average_size = "{size:.3f}".format(size=self.average_size()) 275 276 help_text = "" 277 if self.name: 278 help_text += f"<b> {self.name}:   </b>" 279 help_text += '<b><a href="' + help_url + '" target="_blank">' + library_name + "</a></b>" 280 if self.filename: 281 dots = "" 282 if len(self.filename) > 30: 283 dots = "..." 284 help_text += f"<br/><code><i>({dots}{self.filename[-30:]})</i></code>" 285 286 pdata = "" 287 if self._data.GetPointData().GetScalars(): 288 if self._data.GetPointData().GetScalars().GetName(): 289 name = self._data.GetPointData().GetScalars().GetName() 290 pdata = "<tr><td><b> point data array </b></td><td>" + name + "</td></tr>" 291 292 cdata = "" 293 if self._data.GetCellData().GetScalars(): 294 if self._data.GetCellData().GetScalars().GetName(): 295 name = self._data.GetCellData().GetScalars().GetName() 296 cdata = "<tr><td><b> cell data array </b></td><td>" + name + "</td></tr>" 297 298 all = [ 299 "<table>", 300 "<tr>", 301 "<td>", image, "</td>", 302 "<td style='text-align: center; vertical-align: center;'><br/>", help_text, 303 "<table>", 304 "<tr><td><b> bounds </b> <br/> (x/y/z) </td><td>" + str(bounds) + "</td></tr>", 305 "<tr><td><b> center of mass </b></td><td>" + precision(self.center_of_mass(),3) + "</td></tr>", 306 "<tr><td><b> average size </b></td><td>" + str(average_size) + "</td></tr>", 307 "<tr><td><b> nr. points / faces </b></td><td>" 308 + str(self.npoints) + " / " + str(self.ncells) + "</td></tr>", 309 pdata, 310 cdata, 311 "</table>", 312 "</table>", 313 ] 314 return "\n".join(all) 315 316 317 def faces(self): 318 """ 319 Get cell polygonal connectivity ids as a python `list`. 320 The output format is: `[[id0 ... idn], [id0 ... idm], etc]`. 321 """ 322 arr1d = vtk2numpy(self._data.GetPolys().GetData()) 323 if arr1d is None: 324 return [] 325 326 # Get cell connettivity ids as a 1D array. vtk format is: 327 # [nids1, id0 ... idn, niids2, id0 ... idm, etc]. 328 if len(arr1d) == 0: 329 arr1d = vtk2numpy(self._data.GetStrips().GetData()) 330 if arr1d is None: 331 return [] 332 333 i = 0 334 conn = [] 335 n = len(arr1d) 336 if n: 337 while True: 338 cell = [arr1d[i + k] for k in range(1, arr1d[i] + 1)] 339 conn.append(cell) 340 i += arr1d[i] + 1 341 if i >= n: 342 break 343 return conn # cannot always make a numpy array of it! 344 345 def cells(self): 346 """Alias for `faces()`.""" 347 return self.faces() 348 349 def lines(self, flat=False): 350 """ 351 Get lines connectivity ids as a numpy array. 352 Default format is `[[id0,id1], [id3,id4], ...]` 353 354 Arguments: 355 flat : (bool) 356 return a 1D numpy array as e.g. [2, 10,20, 3, 10,11,12, 2, 70,80, ...] 357 """ 358 # Get cell connettivity ids as a 1D array. The vtk format is: 359 # [nids1, id0 ... idn, niids2, id0 ... idm, etc]. 360 arr1d = vtk2numpy(self.polydata(False).GetLines().GetData()) 361 362 if arr1d is None: 363 return [] 364 365 if flat: 366 return arr1d 367 368 i = 0 369 conn = [] 370 n = len(arr1d) 371 for _ in range(n): 372 cell = [arr1d[i + k + 1] for k in range(arr1d[i])] 373 conn.append(cell) 374 i += arr1d[i] + 1 375 if i >= n: 376 break 377 378 return conn # cannot always make a numpy array of it! 379 380 def edges(self): 381 """Return an array containing the edges connectivity.""" 382 extractEdges = vtk.vtkExtractEdges() 383 extractEdges.SetInputData(self._data) 384 # eed.UseAllPointsOn() 385 extractEdges.Update() 386 lpoly = extractEdges.GetOutput() 387 388 arr1d = vtk2numpy(lpoly.GetLines().GetData()) 389 # [nids1, id0 ... idn, niids2, id0 ... idm, etc]. 390 391 i = 0 392 conn = [] 393 n = len(arr1d) 394 for _ in range(n): 395 cell = [arr1d[i + k + 1] for k in range(arr1d[i])] 396 conn.append(cell) 397 i += arr1d[i] + 1 398 if i >= n: 399 break 400 return conn # cannot always make a numpy array of it! 401 402 def texture( 403 self, 404 tname, 405 tcoords=None, 406 interpolate=True, 407 repeat=True, 408 edge_clamp=False, 409 scale=None, 410 ushift=None, 411 vshift=None, 412 seam_threshold=None, 413 ): 414 """ 415 Assign a texture to mesh from image file or predefined texture `tname`. 416 If tname is set to `None` texture is disabled. 417 Input tname can also be an array or a `vtkTexture`. 418 419 Arguments: 420 tname : (numpy.array, str, Picture, vtkTexture, None) 421 the input texture to be applied. Can be a numpy array, a path to an image file, 422 a vedo Picture. The None value disables texture. 423 tcoords : (numpy.array, str) 424 this is the (u,v) texture coordinate array. Can also be a string of an existing array 425 in the mesh. 426 interpolate : (bool) 427 turn on/off linear interpolation of the texture map when rendering. 428 repeat : (bool) 429 repeat of the texture when tcoords extend beyond the [0,1] range. 430 edge_clamp : (bool) 431 turn on/off the clamping of the texture map when 432 the texture coords extend beyond the [0,1] range. 433 Only used when repeat is False, and edge clamping is supported by the graphics card. 434 scale : (bool) 435 scale the texture image by this factor 436 ushift : (bool) 437 shift u-coordinates of texture by this amount 438 vshift : (bool) 439 shift v-coordinates of texture by this amount 440 seam_threshold : (float) 441 try to seal seams in texture by collapsing triangles 442 (test values around 1.0, lower values = stronger collapse) 443 444 Examples: 445 - [texturecubes.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/texturecubes.py) 446 447  448 """ 449 pd = self.polydata(False) 450 451 if tname is None: # disable texture 452 pd.GetPointData().SetTCoords(None) 453 pd.GetPointData().Modified() 454 return self ###################################### 455 456 if isinstance(tname, vtk.vtkTexture): 457 tu = tname 458 outimg = self._data 459 460 elif isinstance(tname, vedo.Picture): 461 tu = vtk.vtkTexture() 462 outimg = tname.inputdata() 463 464 elif is_sequence(tname): 465 tu = vtk.vtkTexture() 466 outimg = vedo.picture._get_img(tname) 467 468 elif isinstance(tname, str): 469 tu = vtk.vtkTexture() 470 471 if "https://" in tname: 472 try: 473 tname = vedo.io.download(tname, verbose=False) 474 except: 475 vedo.logger.error(f"texture {tname} could not be downloaded") 476 return self 477 478 fn = tname + ".jpg" 479 if os.path.exists(tname): 480 fn = tname 481 else: 482 vedo.logger.error(f"texture file {tname} does not exist") 483 return self 484 485 fnl = fn.lower() 486 if ".jpg" in fnl or ".jpeg" in fnl: 487 reader = vtk.vtkJPEGReader() 488 elif ".png" in fnl: 489 reader = vtk.vtkPNGReader() 490 elif ".bmp" in fnl: 491 reader = vtk.vtkBMPReader() 492 else: 493 vedo.logger.error( 494 "in texture() supported files are only PNG, BMP or JPG" 495 ) 496 return self 497 reader.SetFileName(fn) 498 reader.Update() 499 outimg = reader.GetOutput() 500 501 else: 502 vedo.logger.error(f"in texture() cannot understand input {type(tname)}") 503 return self 504 505 if tcoords is not None: 506 507 if isinstance(tcoords, str): 508 509 vtarr = pd.GetPointData().GetArray(tcoords) 510 511 else: 512 513 tcoords = np.asarray(tcoords) 514 if tcoords.ndim != 2: 515 vedo.logger.error("tcoords must be a 2-dimensional array") 516 return self 517 if tcoords.shape[0] != pd.GetNumberOfPoints(): 518 vedo.logger.error("nr of texture coords must match nr of points") 519 return self 520 if tcoords.shape[1] != 2: 521 vedo.logger.error("tcoords texture vector must have 2 components") 522 vtarr = numpy2vtk(tcoords) 523 vtarr.SetName("TCoordinates") 524 525 pd.GetPointData().SetTCoords(vtarr) 526 pd.GetPointData().Modified() 527 528 elif not pd.GetPointData().GetTCoords(): 529 530 # TCoords still void.. 531 # check that there are no texture-like arrays: 532 names = self.pointdata.keys() 533 candidate_arr = "" 534 for name in names: 535 vtarr = pd.GetPointData().GetArray(name) 536 if vtarr.GetNumberOfComponents() != 2: 537 continue 538 t0, t1 = vtarr.GetRange() 539 if t0 >= 0 and t1 <= 1: 540 candidate_arr = name 541 542 if candidate_arr: 543 544 vtarr = pd.GetPointData().GetArray(candidate_arr) 545 pd.GetPointData().SetTCoords(vtarr) 546 pd.GetPointData().Modified() 547 548 else: 549 # last resource is automatic mapping 550 tmapper = vtk.vtkTextureMapToPlane() 551 tmapper.AutomaticPlaneGenerationOn() 552 tmapper.SetInputData(pd) 553 tmapper.Update() 554 tc = tmapper.GetOutput().GetPointData().GetTCoords() 555 if scale or ushift or vshift: 556 ntc = vtk2numpy(tc) 557 if scale: 558 ntc *= scale 559 if ushift: 560 ntc[:, 0] += ushift 561 if vshift: 562 ntc[:, 1] += vshift 563 tc = numpy2vtk(tc) 564 pd.GetPointData().SetTCoords(tc) 565 pd.GetPointData().Modified() 566 567 tu.SetInputData(outimg) 568 tu.SetInterpolate(interpolate) 569 tu.SetRepeat(repeat) 570 tu.SetEdgeClamp(edge_clamp) 571 572 self.property.SetColor(1, 1, 1) 573 self._mapper.ScalarVisibilityOff() 574 self.SetTexture(tu) 575 576 if seam_threshold is not None: 577 tname = self._data.GetPointData().GetTCoords().GetName() 578 grad = self.gradient(tname) 579 ugrad, vgrad = np.split(grad, 2, axis=1) 580 ugradm, vgradm = vedo.utils.mag2(ugrad), vedo.utils.mag2(vgrad) 581 gradm = np.log(ugradm + vgradm) 582 largegrad_ids = np.arange(len(grad))[gradm > seam_threshold * 4] 583 uvmap = self.pointdata[tname] 584 # collapse triangles that have large gradient 585 new_points = self.points(transformed=False) 586 for f in self.faces(): 587 if np.isin(f, largegrad_ids).all(): 588 id1, id2, id3 = f 589 uv1, uv2, uv3 = uvmap[f] 590 d12 = vedo.mag2(uv1 - uv2) 591 d23 = vedo.mag2(uv2 - uv3) 592 d31 = vedo.mag2(uv3 - uv1) 593 idm = np.argmin([d12, d23, d31]) 594 if idm == 0: 595 new_points[id1] = new_points[id3] 596 new_points[id2] = new_points[id3] 597 elif idm == 1: 598 new_points[id2] = new_points[id1] 599 new_points[id3] = new_points[id1] 600 self.points(new_points) 601 602 self.Modified() 603 return self 604 605 @deprecated(reason=vedo.colors.red + "Please use compute_normals()" + vedo.colors.reset) 606 def computeNormals(self, points=True, cells=True, featureAngle=None, consistency=True): 607 "Deprecated. Please use compute_normals()" 608 return self.compute_normals(points, cells, featureAngle, consistency) 609 610 def compute_normals(self, points=True, cells=True, feature_angle=None, consistency=True): 611 """ 612 Compute cell and vertex normals for the mesh. 613 614 Arguments: 615 points : (bool) 616 do the computation for the vertices too 617 cells : (bool) 618 do the computation for the cells too 619 feature_angle : (float) 620 specify the angle that defines a sharp edge. 621 If the difference in angle across neighboring polygons is greater than this value, 622 the shared edge is considered "sharp" and it is split. 623 consistency : (bool) 624 turn on/off the enforcement of consistent polygon ordering. 625 626 .. warning:: 627 If feature_angle is set to a float the Mesh can be modified, and it 628 can have a different nr. of vertices from the original. 629 """ 630 poly = self.polydata(False) 631 pdnorm = vtk.vtkPolyDataNormals() 632 pdnorm.SetInputData(poly) 633 pdnorm.SetComputePointNormals(points) 634 pdnorm.SetComputeCellNormals(cells) 635 pdnorm.SetConsistency(consistency) 636 pdnorm.FlipNormalsOff() 637 if feature_angle: 638 pdnorm.SetSplitting(True) 639 pdnorm.SetFeatureAngle(feature_angle) 640 else: 641 pdnorm.SetSplitting(False) 642 # print(pdnorm.GetNonManifoldTraversal()) 643 pdnorm.Update() 644 return self._update(pdnorm.GetOutput()) 645 646 def reverse(self, cells=True, normals=False): 647 """ 648 Reverse the order of polygonal cells 649 and/or reverse the direction of point and cell normals. 650 Two flags are used to control these operations: 651 652 - `cells=True` reverses the order of the indices in the cell connectivity list. 653 If cell is a list of IDs only those cells will be reversed. 654 655 - `normals=True` reverses the normals by multiplying the normal vector by -1 656 (both point and cell normals, if present). 657 """ 658 poly = self.polydata(False) 659 660 if is_sequence(cells): 661 for cell in cells: 662 poly.ReverseCell(cell) 663 poly.GetCellData().Modified() 664 return self ############## 665 666 rev = vtk.vtkReverseSense() 667 if cells: 668 rev.ReverseCellsOn() 669 else: 670 rev.ReverseCellsOff() 671 if normals: 672 rev.ReverseNormalsOn() 673 else: 674 rev.ReverseNormalsOff() 675 rev.SetInputData(poly) 676 rev.Update() 677 out = self._update(rev.GetOutput()) 678 out.pipeline = OperationNode("reverse", parents=[self]) 679 return out 680 681 def wireframe(self, value=True): 682 """Set mesh's representation as wireframe or solid surface.""" 683 if value: 684 self.property.SetRepresentationToWireframe() 685 else: 686 self.property.SetRepresentationToSurface() 687 return self 688 689 def flat(self): 690 """Set surface interpolation to flat. 691 692 <img src="https://upload.wikimedia.org/wikipedia/commons/6/6b/Phong_components_version_4.png" width="700"> 693 """ 694 self.property.SetInterpolationToFlat() 695 return self 696 697 def phong(self): 698 """Set surface interpolation to "phong".""" 699 self.property.SetInterpolationToPhong() 700 return self 701 702 def backface_culling(self, value=True): 703 """Set culling of polygons based on orientation of normal with respect to camera.""" 704 self.property.SetBackfaceCulling(value) 705 return self 706 707 def render_lines_as_tubes(self, value=True): 708 """Wrap a fake tube around a simple line for visualization""" 709 self.property.SetRenderLinesAsTubes(value) 710 return self 711 712 def frontface_culling(self, value=True): 713 """Set culling of polygons based on orientation of normal with respect to camera.""" 714 self.property.SetFrontfaceCulling(value) 715 return self 716 717 @deprecated(reason=vedo.colors.red + "Please use backcolor()" + vedo.colors.reset) 718 def backColor(self, bc=None): 719 "Deprecated. Please use `backcolor()`" 720 return self.backcolor(bc) 721 722 def backcolor(self, bc=None): 723 """ 724 Set/get mesh's backface color. 725 """ 726 backProp = self.GetBackfaceProperty() 727 728 if bc is None: 729 if backProp: 730 return backProp.GetDiffuseColor() 731 return self 732 733 if self.property.GetOpacity() < 1: 734 return self 735 736 if not backProp: 737 backProp = vtk.vtkProperty() 738 739 backProp.SetDiffuseColor(get_color(bc)) 740 backProp.SetOpacity(self.property.GetOpacity()) 741 self.SetBackfaceProperty(backProp) 742 self._mapper.ScalarVisibilityOff() 743 return self 744 745 def bc(self, backcolor=False): 746 """Shortcut for `mesh.backcolor()`.""" 747 return self.backcolor(backcolor) 748 749 @deprecated(reason=vedo.colors.red + "Please use linewidth()" + vedo.colors.reset) 750 def lineWidth(self, lw): 751 "Deprecated: please use `linewidth()`" 752 return self.linewidth(lw) 753 754 def linewidth(self, lw=None): 755 """Set/get width of mesh edges. Same as `lw()`.""" 756 if lw is not None: 757 if lw == 0: 758 self.property.EdgeVisibilityOff() 759 self.property.SetRepresentationToSurface() 760 return self 761 self.property.EdgeVisibilityOn() 762 self.property.SetLineWidth(lw) 763 else: 764 return self.property.GetLineWidth() 765 return self 766 767 def lw(self, linewidth=None): 768 """Set/get width of mesh edges. Same as `linewidth()`.""" 769 return self.linewidth(linewidth) 770 771 def linecolor(self, lc=None): 772 """Set/get color of mesh edges. Same as `lc()`.""" 773 if lc is None: 774 return self.property.GetEdgeColor() 775 self.property.EdgeVisibilityOn() 776 self.property.SetEdgeColor(get_color(lc)) 777 return self 778 779 def lc(self, linecolor=None): 780 """Set/get color of mesh edges. Same as `linecolor()`.""" 781 return self.linecolor(linecolor) 782 783 def volume(self): 784 """Get/set the volume occupied by mesh.""" 785 mass = vtk.vtkMassProperties() 786 mass.SetGlobalWarningDisplay(0) 787 mass.SetInputData(self.polydata()) 788 mass.Update() 789 return mass.GetVolume() 790 791 def area(self): 792 """Get/set the surface area of mesh.""" 793 mass = vtk.vtkMassProperties() 794 mass.SetGlobalWarningDisplay(0) 795 mass.SetInputData(self.polydata()) 796 mass.Update() 797 return mass.GetSurfaceArea() 798 799 def is_closed(self): 800 """Return `True` if the mesh is watertight.""" 801 fe = vtk.vtkFeatureEdges() 802 fe.BoundaryEdgesOn() 803 fe.FeatureEdgesOff() 804 fe.NonManifoldEdgesOn() 805 fe.SetInputData(self.polydata(False)) 806 fe.Update() 807 ne = fe.GetOutput().GetNumberOfCells() 808 return not bool(ne) 809 810 def is_manifold(self): 811 """Return `True` if the mesh is manifold.""" 812 fe = vtk.vtkFeatureEdges() 813 fe.BoundaryEdgesOff() 814 fe.FeatureEdgesOff() 815 fe.NonManifoldEdgesOn() 816 fe.SetInputData(self.polydata(False)) 817 fe.Update() 818 ne = fe.GetOutput().GetNumberOfCells() 819 return not bool(ne) 820 821 def non_manifold_faces(self, remove=True, tol="auto"): 822 """ 823 Detect and (try to) remove non-manifold faces of a triangular mesh. 824 825 Set `remove` to `False` to mark cells without removing them. 826 Set `tol=0` for zero-tolerance, the result will be manifold but with holes. 827 Set `tol>0` to cut off non-manifold faces, and try to recover the good ones. 828 Set `tol="auto"` to make an automatic choice of the tolerance. 829 """ 830 # mark original point and cell ids 831 self.add_ids() 832 toremove = self.boundaries( 833 boundary_edges=False, 834 non_manifold_edges=True, 835 cell_edge=True, 836 return_cell_ids=True, 837 ) 838 if len(toremove) == 0: 839 return self 840 841 points = self.points() 842 faces = self.faces() 843 centers = self.cell_centers() 844 845 copy = self.clone() 846 copy.delete_cells(toremove).clean() 847 copy.compute_normals(cells=False) 848 normals = copy.normals() 849 deltas, deltas_i = [], [] 850 851 for i in vedo.utils.progressbar(toremove, delay=3, title="recover faces"): 852 pids = copy.closest_point(centers[i], n=3, return_point_id=True) 853 norms = normals[pids] 854 n = np.mean(norms, axis=0) 855 dn = np.linalg.norm(n) 856 if not dn: 857 continue 858 n = n / dn 859 860 p0, p1, p2 = points[faces[i]][:3] 861 v = np.cross(p1-p0, p2-p0) 862 lv = np.linalg.norm(v) 863 if not lv: 864 continue 865 v = v / lv 866 867 cosa = 1 - np.dot(n, v) 868 deltas.append(cosa) 869 deltas_i.append(i) 870 871 recover = [] 872 if len(deltas): 873 mean_delta = np.mean(deltas) 874 err_delta = np.std(deltas) 875 txt = "" 876 if tol == "auto": #automatic choice 877 tol = mean_delta / 5 878 txt = f"\n Automatic tol. : {tol: .4f}" 879 for i, cosa in zip(deltas_i, deltas): 880 if cosa < tol: 881 recover.append(i) 882 883 vedo.logger.info( 884 f"\n --------- Non manifold faces ---------" 885 f"\n Average tol. : {mean_delta: .4f} +- {err_delta: .4f}{txt}" 886 f"\n Removed faces : {len(toremove)}" 887 f"\n Recovered faces: {len(recover)}" 888 ) 889 890 toremove = list(set(toremove) - set(recover)) 891 892 if not remove: 893 mark = np.zeros(self.ncells, dtype=np.uint8) 894 mark[recover] = 1 895 mark[toremove] = 2 896 self.celldata["NonManifoldCell"] = mark 897 else: 898 self.delete_cells(toremove) 899 900 self.pipeline = OperationNode( 901 "non_manifold_faces", parents=[self], 902 comment=f"#cells {self._data.GetNumberOfCells()}", 903 ) 904 return self 905 906 def shrink(self, fraction=0.85): 907 """Shrink the triangle polydata in the representation of the input mesh. 908 909 Examples: 910 - [shrink.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/shrink.py) 911 912  913 """ 914 shrink = vtk.vtkShrinkPolyData() 915 shrink.SetInputData(self._data) 916 shrink.SetShrinkFactor(fraction) 917 shrink.Update() 918 self.point_locator = None 919 self.cell_locator = None 920 out = self._update(shrink.GetOutput()) 921 922 out.pipeline = OperationNode("shrink", parents=[self]) 923 return out 924 925 def stretch(self, q1, q2): 926 """ 927 Stretch mesh between points `q1` and `q2`. 928 929 Examples: 930 - [aspring.py](https://github.com/marcomusy/vedo/tree/master/examples/simulations/aspring.py) 931 932  933 934 .. note:: 935 for `Mesh` objects, two vectors `mesh.base`, and `mesh.top` must be defined. 936 """ 937 if self.base is None: 938 vedo.logger.error( 939 "in stretch() must define vectors mesh.base and mesh.top at creation" 940 ) 941 raise RuntimeError() 942 943 p1, p2 = self.base, self.top 944 q1, q2, z = np.asarray(q1), np.asarray(q2), np.array([0, 0, 1]) 945 a = p2 - p1 946 b = q2 - q1 947 plength = np.linalg.norm(a) 948 qlength = np.linalg.norm(b) 949 T = vtk.vtkTransform() 950 T.PostMultiply() 951 T.Translate(-p1) 952 cosa = np.dot(a, z) / plength 953 n = np.cross(a, z) 954 if np.linalg.norm(n): 955 T.RotateWXYZ(np.rad2deg(np.arccos(cosa)), n) 956 T.Scale(1, 1, qlength / plength) 957 958 cosa = np.dot(b, z) / qlength 959 n = np.cross(b, z) 960 if np.linalg.norm(n): 961 T.RotateWXYZ(-np.rad2deg(np.arccos(cosa)), n) 962 else: 963 if np.dot(b, z) < 0: 964 T.RotateWXYZ(180, [1,0,0]) 965 966 T.Translate(q1) 967 968 self.SetUserMatrix(T.GetMatrix()) 969 return self 970 971 def crop( 972 self, 973 top=None, 974 bottom=None, 975 right=None, 976 left=None, 977 front=None, 978 back=None, 979 ): 980 """ 981 Crop an `Mesh` object. 982 Use this method at creation (before moving the object). 983 984 Arguments: 985 top : (float) 986 fraction to crop from the top plane (positive z) 987 bottom : (float) 988 fraction to crop from the bottom plane (negative z) 989 front : (float) 990 fraction to crop from the front plane (positive y) 991 back : (float) 992 fraction to crop from the back plane (negative y) 993 right : (float) 994 fraction to crop from the right plane (positive x) 995 left : (float) 996 fraction to crop from the left plane (negative x) 997 998 Example: 999 ```python 1000 from vedo import Sphere 1001 Sphere().crop(right=0.3, left=0.1).show() 1002 ``` 1003  1004 """ 1005 cu = vtk.vtkBox() 1006 pos = np.array(self.GetPosition()) 1007 x0, x1, y0, y1, z0, z1 = self.bounds() 1008 x0, y0, z0 = [x0, y0, z0] - pos 1009 x1, y1, z1 = [x1, y1, z1] - pos 1010 1011 dx, dy, dz = x1 - x0, y1 - y0, z1 - z0 1012 if top: 1013 z1 = z1 - top * dz 1014 if bottom: 1015 z0 = z0 + bottom * dz 1016 if front: 1017 y1 = y1 - front * dy 1018 if back: 1019 y0 = y0 + back * dy 1020 if right: 1021 x1 = x1 - right * dx 1022 if left: 1023 x0 = x0 + left * dx 1024 bounds = (x0, x1, y0, y1, z0, z1) 1025 1026 cu.SetBounds(bounds) 1027 1028 clipper = vtk.vtkClipPolyData() 1029 clipper.SetInputData(self._data) 1030 clipper.SetClipFunction(cu) 1031 clipper.InsideOutOn() 1032 clipper.GenerateClippedOutputOff() 1033 clipper.GenerateClipScalarsOff() 1034 clipper.SetValue(0) 1035 clipper.Update() 1036 self._update(clipper.GetOutput()) 1037 self.point_locator = None 1038 1039 self.pipeline = OperationNode( 1040 "crop", parents=[self], 1041 comment=f"#pts {self._data.GetNumberOfPoints()}", 1042 ) 1043 return self 1044 1045 1046 def cap(self, return_cap=False): 1047 """ 1048 Generate a "cap" on a clipped mesh, or caps sharp edges. 1049 1050 Examples: 1051 - [cut_and_cap.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/cut_and_cap.py) 1052 1053  1054 1055 See also: `join()`, `join_segments()`, `slice()`. 1056 """ 1057 poly = self._data 1058 1059 fe = vtk.vtkFeatureEdges() 1060 fe.SetInputData(poly) 1061 fe.BoundaryEdgesOn() 1062 fe.FeatureEdgesOff() 1063 fe.NonManifoldEdgesOff() 1064 fe.ManifoldEdgesOff() 1065 fe.Update() 1066 1067 stripper = vtk.vtkStripper() 1068 stripper.SetInputData(fe.GetOutput()) 1069 stripper.JoinContiguousSegmentsOn() 1070 stripper.Update() 1071 1072 boundaryPoly = vtk.vtkPolyData() 1073 boundaryPoly.SetPoints(stripper.GetOutput().GetPoints()) 1074 boundaryPoly.SetPolys(stripper.GetOutput().GetLines()) 1075 1076 rev = vtk.vtkReverseSense() 1077 rev.ReverseCellsOn() 1078 rev.SetInputData(boundaryPoly) 1079 rev.Update() 1080 1081 tf = vtk.vtkTriangleFilter() 1082 tf.SetInputData(rev.GetOutput()) 1083 tf.Update() 1084 1085 if return_cap: 1086 m = Mesh(tf.GetOutput()) 1087 # assign the same transformation to the copy 1088 m.SetOrigin(self.GetOrigin()) 1089 m.SetScale(self.GetScale()) 1090 m.SetOrientation(self.GetOrientation()) 1091 m.SetPosition(self.GetPosition()) 1092 1093 m.pipeline = OperationNode( 1094 "cap", parents=[self], 1095 comment=f"#pts {m._data.GetNumberOfPoints()}", 1096 ) 1097 return m 1098 1099 polyapp = vtk.vtkAppendPolyData() 1100 polyapp.AddInputData(poly) 1101 polyapp.AddInputData(tf.GetOutput()) 1102 polyapp.Update() 1103 out = self._update(polyapp.GetOutput()).clean() 1104 1105 out.pipeline = OperationNode( 1106 "capped", parents=[self], 1107 comment=f"#pts {out._data.GetNumberOfPoints()}", 1108 ) 1109 return out 1110 1111 def join(self, polys=True, reset=False): 1112 """ 1113 Generate triangle strips and/or polylines from 1114 input polygons, triangle strips, and lines. 1115 1116 Input polygons are assembled into triangle strips only if they are triangles; 1117 other types of polygons are passed through to the output and not stripped. 1118 Use mesh.triangulate() to triangulate non-triangular polygons prior to running 1119 this filter if you need to strip all the data. 1120 1121 Also note that if triangle strips or polylines are present in the input 1122 they are passed through and not joined nor extended. 1123 If you wish to strip these use mesh.triangulate() to fragment the input 1124 into triangles and lines prior to applying join(). 1125 1126 Arguments: 1127 polys : (bool) 1128 polygonal segments will be joined if they are contiguous 1129 reset : (bool) 1130 reset points ordering 1131 1132 Warning: 1133 If triangle strips or polylines exist in the input data 1134 they will be passed through to the output data. 1135 This filter will only construct triangle strips if triangle polygons 1136 are available; and will only construct polylines if lines are available. 1137 1138 Example: 1139 ```python 1140 from vedo import * 1141 c1 = Cylinder(pos=(0,0,0), r=2, height=3, axis=(1,.0,0), alpha=.1).triangulate() 1142 c2 = Cylinder(pos=(0,0,2), r=1, height=2, axis=(0,.3,1), alpha=.1).triangulate() 1143 intersect = c1.intersect_with(c2).join(reset=True) 1144 spline = Spline(intersect).c('blue').lw(5) 1145 show(c1, c2, spline, intersect.labels('id'), axes=1).close() 1146 ``` 1147  1148 """ 1149 sf = vtk.vtkStripper() 1150 sf.SetPassThroughCellIds(True) 1151 sf.SetPassThroughPointIds(True) 1152 sf.SetJoinContiguousSegments(polys) 1153 sf.SetInputData(self.polydata(False)) 1154 sf.Update() 1155 if reset: 1156 poly = sf.GetOutput() 1157 cpd = vtk.vtkCleanPolyData() 1158 cpd.PointMergingOn() 1159 cpd.ConvertLinesToPointsOn() 1160 cpd.ConvertPolysToLinesOn() 1161 cpd.ConvertStripsToPolysOn() 1162 cpd.SetInputData(poly) 1163 cpd.Update() 1164 poly = cpd.GetOutput() 1165 vpts = poly.GetCell(0).GetPoints().GetData() 1166 poly.GetPoints().SetData(vpts) 1167 return self._update(poly) 1168 1169 out = self._update(sf.GetOutput()) 1170 out.pipeline = OperationNode("join", parents=[self], 1171 comment=f"#pts {out._data.GetNumberOfPoints()}", 1172 ) 1173 return out 1174 1175 def join_segments(self, closed=True, tol=1e-03): 1176 """ 1177 Join line segments into contiguous lines. 1178 Useful to call with `triangulate()` method. 1179 1180 Returns: 1181 list of `shapes.Lines` 1182 1183 Example: 1184 ```python 1185 from vedo import * 1186 msh = Torus().alpha(0.1).wireframe() 1187 intersection = msh.intersect_with_plane(normal=[1,1,1]).c('purple5') 1188 slices = [s.triangulate() for s in intersection.join_segments()] 1189 show(msh, intersection, merge(slices), axes=1, viewup='z') 1190 ``` 1191  1192 """ 1193 vlines = [] 1194 for ipiece, outline in enumerate(self.split()): 1195 1196 outline.clean() 1197 pts = outline.points() 1198 if len(pts) < 3: 1199 continue 1200 avesize = outline.average_size() 1201 lines = outline.lines() 1202 # print("---lines", lines, "in piece", ipiece) 1203 tol = avesize / pts.shape[0] * tol 1204 1205 k = 0 1206 joinedpts = [pts[k]] 1207 for _ in range(len(pts)): 1208 pk = pts[k] 1209 for j, line in enumerate(lines): 1210 1211 id0, id1 = line[0], line[-1] 1212 p0, p1 = pts[id0], pts[id1] 1213 1214 if np.linalg.norm(p0 - pk) < tol: 1215 n = len(line) 1216 for m in range(1, n): 1217 joinedpts.append(pts[line[m]]) 1218 # joinedpts.append(p1) 1219 k = id1 1220 lines.pop(j) 1221 break 1222 1223 elif np.linalg.norm(p1 - pk) < tol: 1224 n = len(line) 1225 for m in reversed(range(0, n-1)): 1226 joinedpts.append(pts[line[m]]) 1227 # joinedpts.append(p0) 1228 k = id0 1229 lines.pop(j) 1230 break 1231 1232 if len(joinedpts) > 1: 1233 newline = vedo.shapes.Line(joinedpts, closed=closed) 1234 newline.clean() 1235 newline.SetProperty(self.GetProperty()) 1236 newline.property = self.GetProperty() 1237 newline.pipeline = OperationNode("join_segments", parents=[self], 1238 comment=f"#pts {newline._data.GetNumberOfPoints()}", 1239 ) 1240 vlines.append(newline) 1241 1242 return vlines 1243 1244 def slice(self, origin=[0,0,0], normal=[1,0,0]): 1245 """ 1246 Slice a mesh with a plane and fill the contour. 1247 1248 Example: 1249 ```python 1250 from vedo import * 1251 msh = Mesh(dataurl+"bunny.obj").alpha(0.1).wireframe() 1252 mslice = msh.slice(normal=[0,1,0.3], origin=[0,0.16,0]) 1253 mslice.c('purple5') 1254 show(msh, mslice, axes=1) 1255 ``` 1256  1257 1258 See also: `join()`, `join_segments()`, `cap()`, `cut_with_plane()`. 1259 """ 1260 intersection = self.intersect_with_plane(origin=origin, normal=normal) 1261 slices = [s.triangulate() for s in intersection.join_segments()] 1262 mslices = vedo.pointcloud.merge(slices) 1263 if mslices: 1264 mslices.name = "MeshSlice" 1265 mslices.pipeline = OperationNode( 1266 "slice", parents=[self], comment=f"normal = {normal}") 1267 return mslices 1268 1269 1270 def triangulate(self, verts=True, lines=True): 1271 """ 1272 Converts mesh polygons into triangles. 1273 1274 If the input mesh is only made of 2D lines (no faces) the output will be a triangulation 1275 that fills the internal area. The contours may be concave, and may even contain holes, 1276 i.e. a contour may contain an internal contour winding in the opposite 1277 direction to indicate that it is a hole. 1278 1279 Arguments: 1280 verts : (bool) 1281 if True, break input vertex cells into individual vertex cells 1282 (one point per cell). If False, the input vertex cells will be ignored. 1283 lines : (bool) 1284 if True, break input polylines into line segments. 1285 If False, input lines will be ignored and the output will have no lines. 1286 """ 1287 if self._data.GetNumberOfPolys() or self._data.GetNumberOfStrips(): 1288 tf = vtk.vtkTriangleFilter() 1289 tf.SetPassLines(lines) 1290 tf.SetPassVerts(verts) 1291 1292 elif self._data.GetNumberOfLines(): 1293 tf = vtk.vtkContourTriangulator() 1294 1295 else: 1296 vedo.logger.debug("input in triangulate() seems to be void! Skip.") 1297 return self 1298 1299 tf.SetInputData(self._data) 1300 tf.Update() 1301 out = self._update(tf.GetOutput()) 1302 1303 out.pipeline = OperationNode( 1304 "triangulate", parents=[self], comment=f"#cells {out._data.GetNumberOfCells()}", 1305 ) 1306 return out 1307 1308 1309 def compute_cell_area(self, name="Area"): 1310 """Add to this mesh a cell data array containing the areas of the polygonal faces""" 1311 csf = vtk.vtkCellSizeFilter() 1312 csf.SetInputData(self.polydata(False)) 1313 csf.SetComputeArea(True) 1314 csf.SetComputeVolume(False) 1315 csf.SetComputeLength(False) 1316 csf.SetComputeVertexCount(False) 1317 csf.SetAreaArrayName(name) 1318 csf.Update() 1319 return self._update(csf.GetOutput()) 1320 1321 def compute_cell_vertex_count(self, name="VertexCount"): 1322 """Add to this mesh a cell data array containing the nr of vertices 1323 that a polygonal face has.""" 1324 csf = vtk.vtkCellSizeFilter() 1325 csf.SetInputData(self.polydata(False)) 1326 csf.SetComputeArea(False) 1327 csf.SetComputeVolume(False) 1328 csf.SetComputeLength(False) 1329 csf.SetComputeVertexCount(True) 1330 csf.SetVertexCountArrayName(name) 1331 csf.Update() 1332 return self._update(csf.GetOutput()) 1333 1334 def compute_quality(self, metric=6): 1335 """ 1336 Calculate metrics of quality for the elements of a triangular mesh. 1337 This method adds to the mesh a cell array named "Quality". 1338 See class [vtkMeshQuality](https://vtk.org/doc/nightly/html/classvtkMeshQuality.html) 1339 for explanation. 1340 1341 Arguments: 1342 metric : (int) 1343 type of available estimators are: 1344 - EDGE RATIO, 0 1345 - ASPECT RATIO, 1 1346 - RADIUS RATIO, 2 1347 - ASPECT FROBENIUS, 3 1348 - MED ASPECT FROBENIUS, 4 1349 - MAX ASPECT FROBENIUS, 5 1350 - MIN_ANGLE, 6 1351 - COLLAPSE RATIO, 7 1352 - MAX ANGLE, 8 1353 - CONDITION, 9 1354 - SCALED JACOBIAN, 10 1355 - SHEAR, 11 1356 - RELATIVE SIZE SQUARED, 12 1357 - SHAPE, 13 1358 - SHAPE AND SIZE, 14 1359 - DISTORTION, 15 1360 - MAX EDGE RATIO, 16 1361 - SKEW, 17 1362 - TAPER, 18 1363 - VOLUME, 19 1364 - STRETCH, 20 1365 - DIAGONAL, 21 1366 - DIMENSION, 22 1367 - ODDY, 23 1368 - SHEAR AND SIZE, 24 1369 - JACOBIAN, 25 1370 - WARPAGE, 26 1371 - ASPECT GAMMA, 27 1372 - AREA, 28 1373 - ASPECT BETA, 29 1374 1375 Examples: 1376 - [meshquality.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/meshquality.py) 1377 1378  1379 """ 1380 qf = vtk.vtkMeshQuality() 1381 qf.SetInputData(self.polydata(False)) 1382 qf.SetTriangleQualityMeasure(metric) 1383 qf.SaveCellQualityOn() 1384 qf.Update() 1385 pd = qf.GetOutput() 1386 self._update(pd) 1387 self.pipeline = OperationNode("compute_quality", parents=[self]) 1388 return self 1389 1390 def check_validity(self, tol=0): 1391 """ 1392 Return an array of possible problematic faces following this convention: 1393 - Valid = 0 1394 - WrongNumberOfPoints = 1 1395 - IntersectingEdges = 2 1396 - IntersectingFaces = 4 1397 - NoncontiguousEdges = 8 1398 - Nonconvex = 10 1399 - OrientedIncorrectly = 20 1400 1401 Arguments: 1402 tol : (float) 1403 value is used as an epsilon for floating point 1404 equality checks throughout the cell checking process. 1405 """ 1406 vald = vtk.vtkCellValidator() 1407 if tol: 1408 vald.SetTolerance(tol) 1409 vald.SetInputData(self._data) 1410 vald.Update() 1411 varr = vald.GetOutput().GetCellData().GetArray("ValidityState") 1412 return vtk2numpy(varr) 1413 1414 1415 @deprecated(reason=vedo.colors.red + "Please use compute_curvature()" + vedo.colors.reset) 1416 def addCurvatureScalars(self, method=0): 1417 """Deprecated. Please use `compute_curvature()`""" 1418 return self.compute_curvature(method) 1419 1420 def compute_curvature(self, method=0): 1421 """ 1422 Add scalars to `Mesh` that contains the curvature calculated in three different ways. 1423 1424 Variable `method` can be: 1425 - 0 = gaussian 1426 - 1 = mean curvature 1427 - 2 = max curvature 1428 - 3 = min curvature 1429 1430 Example: 1431 ```python 1432 from vedo import Torus 1433 Torus().compute_curvature().add_scalarbar().show(axes=1).close() 1434 ``` 1435  1436 """ 1437 curve = vtk.vtkCurvatures() 1438 curve.SetInputData(self._data) 1439 curve.SetCurvatureType(method) 1440 curve.Update() 1441 self._update(curve.GetOutput()) 1442 self._mapper.ScalarVisibilityOn() 1443 return self 1444 1445 def compute_connectivity(self): 1446 """ 1447 Flag a mesh by connectivity: each disconnected region will receive a different Id. 1448 You can access the array of ids through `mesh.pointdata["RegionId"]`. 1449 """ 1450 cf = vtk.vtkConnectivityFilter() 1451 cf.SetInputData(self.polydata(False)) 1452 cf.SetExtractionModeToAllRegions() 1453 cf.ColorRegionsOn() 1454 cf.Update() 1455 return self._update(cf.GetOutput()) 1456 1457 def compute_elevation(self, low=(0, 0, 0), high=(0, 0, 1), vrange=(0, 1)): 1458 """ 1459 Add to `Mesh` a scalar array that contains distance along a specified direction. 1460 1461 Arguments: 1462 low : (list) 1463 one end of the line (small scalar values) 1464 high : (list) 1465 other end of the line (large scalar values) 1466 vrange : (list) 1467 set the range of the scalar 1468 1469 Example: 1470 ```python 1471 from vedo import Sphere 1472 s = Sphere().compute_elevation(low=(0,0,0), high=(1,1,1)) 1473 s.add_scalarbar().show(axes=1).close() 1474 ``` 1475  1476 """ 1477 ef = vtk.vtkElevationFilter() 1478 ef.SetInputData(self.polydata()) 1479 ef.SetLowPoint(low) 1480 ef.SetHighPoint(high) 1481 ef.SetScalarRange(vrange) 1482 ef.Update() 1483 self._update(ef.GetOutput()) 1484 self._mapper.ScalarVisibilityOn() 1485 return self 1486 1487 @deprecated(reason=vedo.colors.red + "Please use add_shadow()" + vedo.colors.reset) 1488 def addShadow(self, *a, **b): 1489 """Please use `add_shadow()`""" 1490 return self.add_shadow(*a, **b) 1491 1492 def add_shadow( 1493 self, 1494 plane, 1495 point, 1496 direction=None, 1497 c=(0.6, 0.6, 0.6), 1498 alpha=1, 1499 culling=0, 1500 ): 1501 """ 1502 Generate a shadow out of an `Mesh` on one of the three Cartesian planes. 1503 The output is a new `Mesh` representing the shadow. 1504 This new mesh is accessible through `mesh.shadow`. 1505 By default the shadow mesh is placed on the bottom wall of the bounding box. 1506 1507 See also `pointcloud.project_on_plane()`. 1508 1509 Arguments: 1510 plane : (str, Plane) 1511 if plane is `str`, plane can be one of `['x', 'y', 'z']`, 1512 represents x-plane, y-plane and z-plane, respectively. 1513 Otherwise, plane should be an instance of `vedo.shapes.Plane` 1514 point : (float, array) 1515 if plane is `str`, point should be a float represents the intercept. 1516 Otherwise, point is the camera point of perspective projection 1517 direction : (list) 1518 direction of oblique projection 1519 culling : (int) 1520 choose between front [1] or backface [-1] culling or None. 1521 1522 Examples: 1523 - [shadow1.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/shadow1.py) 1524 - [airplane1.py](https://github.com/marcomusy/vedo/tree/master/examples/simulations/airplane1.py) 1525 - [airplane2.py](https://github.com/marcomusy/vedo/tree/master/examples/simulations/airplane2.py) 1526 1527  1528 """ 1529 shad = self.clone() 1530 shad.name = "Shadow" 1531 pts = shad.points() 1532 if "x" == plane: 1533 # shad = shad.project_on_plane('x') 1534 # instead do it manually so in case of alpha<1 we dont see glitches due to coplanar points 1535 # we leave a small tolerance of 0.1% in thickness 1536 x0, x1 = self.xbounds() 1537 pts[:, 0] = (pts[:, 0] - (x0 + x1) / 2) / 1000 + self.GetOrigin()[0] 1538 shad.points(pts) 1539 shad.x(point) 1540 elif "y" == plane: 1541 x0, x1 = self.ybounds() 1542 pts[:, 1] = (pts[:, 1] - (x0 + x1) / 2) / 1000 + self.GetOrigin()[1] 1543 shad.points(pts) 1544 shad.y(point) 1545 elif "z" == plane: 1546 x0, x1 = self.zbounds() 1547 pts[:, 2] = (pts[:, 2] - (x0 + x1) / 2) / 1000 + self.GetOrigin()[2] 1548 shad.points(pts) 1549 shad.z(point) 1550 else: 1551 shad = shad.project_on_plane(plane, point, direction) 1552 1553 shad.c(c).alpha(alpha).flat() 1554 1555 if culling in (1, True): 1556 shad.frontface_culling() 1557 elif culling == -1: 1558 shad.backface_culling() 1559 1560 shad.GetProperty().LightingOff() 1561 shad.SetPickable(False) 1562 shad.SetUseBounds(True) 1563 1564 if shad not in self.shadows: 1565 self.shadows.append(shad) 1566 shad.info = dict(plane=plane, point=point, direction=direction) 1567 return self 1568 1569 def subdivide(self, n=1, method=0, mel=None): 1570 """ 1571 Increase the number of vertices of a surface mesh. 1572 1573 Arguments: 1574 n : (int) 1575 number of subdivisions. 1576 method : (int) 1577 Loop(0), Linear(1), Adaptive(2), Butterfly(3), Centroid(4) 1578 mel : (float) 1579 Maximum Edge Length (applicable to Adaptive method only). 1580 """ 1581 triangles = vtk.vtkTriangleFilter() 1582 triangles.SetInputData(self._data) 1583 triangles.Update() 1584 originalMesh = triangles.GetOutput() 1585 if method == 0: 1586 sdf = vtk.vtkLoopSubdivisionFilter() 1587 elif method == 1: 1588 sdf = vtk.vtkLinearSubdivisionFilter() 1589 elif method == 2: 1590 sdf = vtk.vtkAdaptiveSubdivisionFilter() 1591 if mel is None: 1592 mel = self.diagonal_size() / np.sqrt(self._data.GetNumberOfPoints()) / n 1593 sdf.SetMaximumEdgeLength(mel) 1594 elif method == 3: 1595 sdf = vtk.vtkButterflySubdivisionFilter() 1596 elif method == 4: 1597 sdf = vtk.vtkDensifyPolyData() 1598 else: 1599 vedo.logger.error(f"in subdivide() unknown method {method}") 1600 raise RuntimeError() 1601 1602 if method != 2: 1603 sdf.SetNumberOfSubdivisions(n) 1604 1605 sdf.SetInputData(originalMesh) 1606 sdf.Update() 1607 out = sdf.GetOutput() 1608 1609 self.pipeline = OperationNode( 1610 "subdivide", parents=[self], comment=f"#pts {out.GetNumberOfPoints()}", 1611 ) 1612 return self._update(out) 1613 1614 def decimate(self, fraction=0.5, n=None, method="quadric", boundaries=False): 1615 """ 1616 Downsample the number of vertices in a mesh to `fraction`. 1617 1618 Arguments: 1619 fraction : (float) 1620 the desired target of reduction. 1621 n : (int) 1622 the desired number of final points (`fraction` is recalculated based on it). 1623 method : (str) 1624 can be either 'quadric' or 'pro'. In the first case triagulation 1625 will look like more regular, irrespective of the mesh original curvature. 1626 In the second case triangles are more irregular but mesh is more precise on more 1627 curved regions. 1628 boundaries : (bool) 1629 in "pro" mode decide whether to leave boundaries untouched or not 1630 1631 .. note:: Setting `fraction=0.1` leaves 10% of the original number of vertices 1632 """ 1633 poly = self._data 1634 if n: # N = desired number of points 1635 npt = poly.GetNumberOfPoints() 1636 fraction = n / npt 1637 if fraction >= 1: 1638 return self 1639 1640 if "quad" in method: 1641 decimate = vtk.vtkQuadricDecimation() 1642 # decimate.SetVolumePreservation(True) 1643 1644 else: 1645 decimate = vtk.vtkDecimatePro() 1646 decimate.PreserveTopologyOn() 1647 if boundaries: 1648 decimate.BoundaryVertexDeletionOff() 1649 else: 1650 decimate.BoundaryVertexDeletionOn() 1651 decimate.SetInputData(poly) 1652 decimate.SetTargetReduction(1 - fraction) 1653 decimate.Update() 1654 out = decimate.GetOutput() 1655 1656 self.pipeline = OperationNode( 1657 "decimate", parents=[self], comment=f"#pts {out.GetNumberOfPoints()}", 1658 ) 1659 return self._update(out) 1660 1661 def collapse_edges(self, distance, iterations=1): 1662 """Collapse mesh edges so that are all above distance.""" 1663 self.clean() 1664 x0, x1, y0, y1, z0, z1 = self.bounds() 1665 fs = min(x1 - x0, y1 - y0, z1 - z0) / 10 1666 d2 = distance * distance 1667 if distance > fs: 1668 vedo.logger.error( 1669 f"distance parameter is too large, should be < {fs}, skip!" 1670 ) 1671 return self 1672 for _ in range(iterations): 1673 medges = self.edges() 1674 pts = self.points() 1675 newpts = np.array(pts) 1676 moved = [] 1677 for e in medges: 1678 if len(e) == 2: 1679 id0, id1 = e 1680 p0, p1 = pts[id0], pts[id1] 1681 d = mag2(p1 - p0) 1682 if d < d2 and id0 not in moved and id1 not in moved: 1683 p = (p0 + p1) / 2 1684 newpts[id0] = p 1685 newpts[id1] = p 1686 moved += [id0, id1] 1687 1688 self.points(newpts) 1689 self.clean() 1690 self.compute_normals() 1691 1692 self.pipeline = OperationNode( 1693 "collapse_edges", parents=[self], 1694 comment=f"#pts {self._data.GetNumberOfPoints()}", 1695 ) 1696 return self 1697 1698 1699 def smooth( 1700 self, niter=15, pass_band=0.1, edge_angle=15, feature_angle=60, boundary=False 1701 ): 1702 """ 1703 Adjust mesh point positions using the `Windowed Sinc` function interpolation kernel. 1704 1705 Arguments: 1706 niter : (int) 1707 number of iterations. 1708 pass_band : (float) 1709 set the pass_band value for the windowed sinc filter. 1710 edge_angle : (float) 1711 edge angle to control smoothing along edges (either interior or boundary). 1712 feature_angle : (float) 1713 specifies the feature angle for sharp edge identification. 1714 boundary : (bool) 1715 specify if boundary should also be smoothed or kept unmodified 1716 1717 Examples: 1718 - [mesh_smoother1.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/mesh_smoother1.py) 1719 1720  1721 """ 1722 poly = self._data 1723 cl = vtk.vtkCleanPolyData() 1724 cl.SetInputData(poly) 1725 cl.Update() 1726 smf = vtk.vtkWindowedSincPolyDataFilter() 1727 smf.SetInputData(cl.GetOutput()) 1728 smf.SetNumberOfIterations(niter) 1729 smf.SetEdgeAngle(edge_angle) 1730 smf.SetFeatureAngle(feature_angle) 1731 smf.SetPassBand(pass_band) 1732 smf.NormalizeCoordinatesOn() 1733 smf.NonManifoldSmoothingOn() 1734 smf.FeatureEdgeSmoothingOn() 1735 smf.SetBoundarySmoothing(boundary) 1736 smf.Update() 1737 out = self._update(smf.GetOutput()) 1738 1739 out.pipeline = OperationNode("smooth", parents=[self], 1740 comment=f"#pts {out._data.GetNumberOfPoints()}", 1741 ) 1742 return out 1743 1744 @deprecated(reason=vedo.colors.red + "Please use fill_holes()" + vedo.colors.reset) 1745 def fillHoles(self, size=None): 1746 """Deprecated. Please use `fill_holes()`""" 1747 return self.fill_holes(size) 1748 1749 def fill_holes(self, size=None): 1750 """ 1751 Identifies and fills holes in input mesh. 1752 Holes are identified by locating boundary edges, linking them together into loops, 1753 and then triangulating the resulting loops. 1754 1755 Arguments: 1756 size : (float) 1757 Approximate limit to the size of the hole that can be filled. The default is None. 1758 1759 Examples: 1760 - [fillholes.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/fillholes.py) 1761 """ 1762 fh = vtk.vtkFillHolesFilter() 1763 if not size: 1764 mb = self.diagonal_size() 1765 size = mb / 10 1766 fh.SetHoleSize(size) 1767 fh.SetInputData(self._data) 1768 fh.Update() 1769 out = self._update(fh.GetOutput()) 1770 1771 out.pipeline = OperationNode( 1772 "fill_holes", parents=[self], 1773 comment=f"#pts {out._data.GetNumberOfPoints()}", 1774 ) 1775 return out 1776 1777 def is_inside(self, point, tol=1e-05): 1778 """Return True if point is inside a polydata closed surface.""" 1779 poly = self.polydata() 1780 points = vtk.vtkPoints() 1781 points.InsertNextPoint(point) 1782 pointsPolydata = vtk.vtkPolyData() 1783 pointsPolydata.SetPoints(points) 1784 sep = vtk.vtkSelectEnclosedPoints() 1785 sep.SetTolerance(tol) 1786 sep.CheckSurfaceOff() 1787 sep.SetInputData(pointsPolydata) 1788 sep.SetSurfaceData(poly) 1789 sep.Update() 1790 return sep.IsInside(0) 1791 1792 def inside_points(self, pts, invert=False, tol=1e-05, return_ids=False): 1793 """ 1794 Return the point cloud that is inside mesh surface as a new Points object. 1795 1796 If return_ids is True a list of IDs is returned and in addition input points 1797 are marked by a pointdata array named "IsInside". 1798 1799 Example: 1800 `print(pts.pointdata["IsInside"])` 1801 1802 Examples: 1803 - [pca_ellipsoid.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/pca_ellipsoid.py) 1804 1805  1806 """ 1807 if isinstance(pts, Points): 1808 pointsPolydata = pts.polydata() 1809 ptsa = pts.points() 1810 else: 1811 ptsa = np.asarray(pts) 1812 vpoints = vtk.vtkPoints() 1813 vpoints.SetData(numpy2vtk(ptsa, dtype=np.float32)) 1814 pointsPolydata = vtk.vtkPolyData() 1815 pointsPolydata.SetPoints(vpoints) 1816 1817 sep = vtk.vtkSelectEnclosedPoints() 1818 # sep = vtk.vtkExtractEnclosedPoints() 1819 sep.SetTolerance(tol) 1820 sep.SetInputData(pointsPolydata) 1821 sep.SetSurfaceData(self.polydata()) 1822 sep.SetInsideOut(invert) 1823 sep.Update() 1824 1825 varr = sep.GetOutput().GetPointData().GetArray("SelectedPoints") 1826 mask = vtk2numpy(varr).astype(bool) 1827 ids = np.array(range(len(ptsa)), dtype=int)[mask] 1828 1829 if isinstance(pts, Points): 1830 varr.SetName("IsInside") 1831 pts.inputdata().GetPointData().AddArray(varr) 1832 1833 if return_ids: 1834 return ids 1835 1836 pcl = Points(ptsa[ids]) 1837 pcl.name = "InsidePoints" 1838 1839 pcl.pipeline = OperationNode( 1840 "inside_points", parents=[self, ptsa], 1841 comment=f"#pts {pcl._data.GetNumberOfPoints()}", 1842 ) 1843 return pcl 1844 1845 def boundaries( 1846 self, 1847 boundary_edges=True, 1848 manifold_edges=False, 1849 non_manifold_edges=False, 1850 feature_angle=None, 1851 return_point_ids=False, 1852 return_cell_ids=False, 1853 cell_edge=False, 1854 ): 1855 """ 1856 Return the boundary lines of an input mesh. 1857 Check also `vedo.base.BaseActor.mark_boundaries()` method. 1858 1859 Arguments: 1860 boundary_edges : (bool) 1861 Turn on/off the extraction of boundary edges. 1862 manifold_edges : (bool) 1863 Turn on/off the extraction of manifold edges. 1864 non_manifold_edges : (bool) 1865 Turn on/off the extraction of non-manifold edges. 1866 feature_angle : (bool) 1867 Specify the min angle btw 2 faces for extracting edges. 1868 return_point_ids : (bool) 1869 return a numpy array of point indices 1870 return_cell_ids : (bool) 1871 return a numpy array of cell indices 1872 cell_edge : (bool) 1873 set to `True` if a cell need to share an edge with 1874 the boundary line, or `False` if a single vertex is enough 1875 1876 Examples: 1877 - [boundaries.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/boundaries.py) 1878 1879  1880 """ 1881 fe = vtk.vtkFeatureEdges() 1882 fe.SetBoundaryEdges(boundary_edges) 1883 fe.SetNonManifoldEdges(non_manifold_edges) 1884 fe.SetManifoldEdges(manifold_edges) 1885 # fe.SetPassLines(True) # vtk9.2 1886 fe.ColoringOff() 1887 fe.SetFeatureEdges(False) 1888 if feature_angle is not None: 1889 fe.SetFeatureEdges(True) 1890 fe.SetFeatureAngle(feature_angle) 1891 1892 if return_point_ids or return_cell_ids: 1893 idf = vtk.vtkIdFilter() 1894 idf.SetInputData(self.polydata()) 1895 idf.SetPointIdsArrayName("BoundaryIds") 1896 idf.SetPointIds(True) 1897 idf.Update() 1898 1899 fe.SetInputData(idf.GetOutput()) 1900 fe.Update() 1901 1902 vid = fe.GetOutput().GetPointData().GetArray("BoundaryIds") 1903 npid = vtk2numpy(vid).astype(int) 1904 1905 if return_point_ids: 1906 return npid 1907 1908 if return_cell_ids: 1909 n = 1 if cell_edge else 0 1910 inface = [] 1911 for i, face in enumerate(self.faces()): 1912 # isin = np.any([vtx in npid for vtx in face]) 1913 isin = 0 1914 for vtx in face: 1915 isin += int(vtx in npid) 1916 if isin > n: 1917 break 1918 if isin > n: 1919 inface.append(i) 1920 return np.array(inface).astype(int) 1921 1922 return self 1923 1924 else: 1925 1926 fe.SetInputData(self.polydata()) 1927 fe.Update() 1928 msh = Mesh(fe.GetOutput(), c="p").lw(5).lighting("off") 1929 1930 msh.pipeline = OperationNode( 1931 "boundaries", 1932 parents=[self], 1933 shape="octagon", 1934 comment=f"#pts {msh._data.GetNumberOfPoints()}", 1935 ) 1936 return msh 1937 1938 1939 def imprint(self, loopline, tol=0.01): 1940 """ 1941 Imprint the contact surface of one object onto another surface. 1942 1943 Arguments: 1944 loopline : vedo.shapes.Line 1945 a Line object to be imprinted onto the mesh. 1946 tol : (float), optional 1947 projection tolerance which controls how close the imprint surface must be to the target. 1948 1949 Example: 1950 ```python 1951 from vedo import * 1952 grid = Grid()#.triangulate() 1953 circle = Circle(r=0.3, res=24).pos(0.11,0.12) 1954 line = Line(circle, closed=True, lw=4, c='r4') 1955 grid.imprint(line) 1956 show(grid, line, axes=1).close() 1957 ``` 1958  1959 """ 1960 loop = vtk.vtkContourLoopExtraction() 1961 loop.SetInputData(loopline.polydata()) 1962 loop.Update() 1963 1964 clean_loop = vtk.vtkCleanPolyData() 1965 clean_loop.SetInputData(loop.GetOutput()) 1966 clean_loop.Update() 1967 1968 imp = vtk.vtkImprintFilter() 1969 imp.SetTargetData(self.polydata()) 1970 imp.SetImprintData(clean_loop.GetOutput()) 1971 imp.SetTolerance(tol) 1972 imp.BoundaryEdgeInsertionOn() 1973 imp.TriangulateOutputOn() 1974 imp.Update() 1975 out = self._update(imp.GetOutput()) 1976 1977 out.pipeline = OperationNode("imprint", parents=[self], 1978 comment=f"#pts {out._data.GetNumberOfPoints()}", 1979 ) 1980 return out 1981 1982 def connected_vertices(self, index): 1983 """Find all vertices connected to an input vertex specified by its index. 1984 1985 Examples: 1986 - [connected_vtx.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/connected_vtx.py) 1987 1988  1989 """ 1990 poly = self._data 1991 1992 cell_idlist = vtk.vtkIdList() 1993 poly.GetPointCells(index, cell_idlist) 1994 1995 idxs = [] 1996 for i in range(cell_idlist.GetNumberOfIds()): 1997 point_idlist = vtk.vtkIdList() 1998 poly.GetCellPoints(cell_idlist.GetId(i), point_idlist) 1999 for j in range(point_idlist.GetNumberOfIds()): 2000 idj = point_idlist.GetId(j) 2001 if idj == index: 2002 continue 2003 if idj in idxs: 2004 continue 2005 idxs.append(idj) 2006 2007 return idxs 2008 2009 def connected_cells(self, index, return_ids=False): 2010 """Find all cellls connected to an input vertex specified by its index.""" 2011 2012 # Find all cells connected to point index 2013 dpoly = self._data 2014 idlist = vtk.vtkIdList() 2015 dpoly.GetPointCells(index, idlist) 2016 2017 ids = vtk.vtkIdTypeArray() 2018 ids.SetNumberOfComponents(1) 2019 rids = [] 2020 for k in range(idlist.GetNumberOfIds()): 2021 cid = idlist.GetId(k) 2022 ids.InsertNextValue(cid) 2023 rids.append(int(cid)) 2024 if return_ids: 2025 return rids 2026 2027 selection_node = vtk.vtkSelectionNode() 2028 selection_node.SetFieldType(vtk.vtkSelectionNode.CELL) 2029 selection_node.SetContentType(vtk.vtkSelectionNode.INDICES) 2030 selection_node.SetSelectionList(ids) 2031 selection = vtk.vtkSelection() 2032 selection.AddNode(selection_node) 2033 extractSelection = vtk.vtkExtractSelection() 2034 extractSelection.SetInputData(0, dpoly) 2035 extractSelection.SetInputData(1, selection) 2036 extractSelection.Update() 2037 gf = vtk.vtkGeometryFilter() 2038 gf.SetInputData(extractSelection.GetOutput()) 2039 gf.Update() 2040 return Mesh(gf.GetOutput()).lw(1) 2041 2042 2043 def silhouette(self, direction=None, border_edges=True, feature_angle=False): 2044 """ 2045 Return a new line `Mesh` which corresponds to the outer `silhouette` 2046 of the input as seen along a specified `direction`, this can also be 2047 a `vtkCamera` object. 2048 2049 Arguments: 2050 direction : (list) 2051 viewpoint direction vector. 2052 If *None* this is guessed by looking at the minimum 2053 of the sides of the bounding box. 2054 border_edges : (bool) 2055 enable or disable generation of border edges 2056 feature_angle : (float) 2057 minimal angle for sharp edges detection. 2058 If set to `False` the functionality is disabled. 2059 2060 Examples: 2061 - [silhouette1.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/silhouette1.py) 2062 2063  2064 """ 2065 sil = vtk.vtkPolyDataSilhouette() 2066 sil.SetInputData(self.polydata()) 2067 sil.SetBorderEdges(border_edges) 2068 if feature_angle is False: 2069 sil.SetEnableFeatureAngle(0) 2070 else: 2071 sil.SetEnableFeatureAngle(1) 2072 sil.SetFeatureAngle(feature_angle) 2073 2074 if direction is None and vedo.plotter_instance and vedo.plotter_instance.camera: 2075 sil.SetCamera(vedo.plotter_instance.camera) 2076 m = Mesh() 2077 m.mapper().SetInputConnection(sil.GetOutputPort()) 2078 2079 elif isinstance(direction, vtk.vtkCamera): 2080 sil.SetCamera(direction) 2081 m = Mesh() 2082 m.mapper().SetInputConnection(sil.GetOutputPort()) 2083 2084 elif direction == "2d": 2085 sil.SetVector(3.4, 4.5, 5.6) # random 2086 sil.SetDirectionToSpecifiedVector() 2087 sil.Update() 2088 m = Mesh(sil.GetOutput()) 2089 2090 elif is_sequence(direction): 2091 sil.SetVector(direction) 2092 sil.SetDirectionToSpecifiedVector() 2093 sil.Update() 2094 m = Mesh(sil.GetOutput()) 2095 else: 2096 vedo.logger.error(f"in silhouette() unknown direction type {type(direction)}") 2097 vedo.logger.error("first render the scene with show() or specify camera/direction") 2098 return self 2099 2100 m.lw(2).c((0, 0, 0)).lighting("off") 2101 m.mapper().SetResolveCoincidentTopologyToPolygonOffset() 2102 m.pipeline = OperationNode("silhouette", parents=[self]) 2103 return m 2104 2105 @deprecated(reason=vedo.colors.red + "Please use follow_camera()" + vedo.colors.reset) 2106 def followCamera(self, cam=None): 2107 "Deprecated. Please use `follow_camera()`" 2108 return self.follow_camera(cam) 2109 2110 def follow_camera(self, camera=None): 2111 """ 2112 Return an object that will follow camera movements and stay locked to it. 2113 Use `mesh.follow_camera(False)` to disable it. 2114 2115 A `vtkCamera` object can also be passed. 2116 """ 2117 if camera is False: 2118 try: 2119 self.SetCamera(None) 2120 return self 2121 except AttributeError: 2122 return self 2123 2124 factor = Follower(self, camera) 2125 2126 if isinstance(camera, vtk.vtkCamera): 2127 factor.SetCamera(camera) 2128 else: 2129 plt = vedo.plotter_instance 2130 if plt and plt.renderer and plt.renderer.GetActiveCamera(): 2131 factor.SetCamera(plt.renderer.GetActiveCamera()) 2132 else: 2133 factor._isfollower = True # postpone to show() call 2134 2135 return factor 2136 2137 def isobands(self, n=10, vmin=None, vmax=None): 2138 """ 2139 Return a new `Mesh` representing the isobands of the active scalars. 2140 This is a new mesh where the scalar is now associated to cell faces and 2141 used to colorize the mesh. 2142 2143 Arguments: 2144 n : (int) 2145 number of isobands in the range 2146 vmin : (float) 2147 minimum of the range 2148 vmax : (float) 2149 maximum of the range 2150 2151 Examples: 2152 - [isolines.py](https://github.com/marcomusy/vedo/tree/master/examples/pyplot/isolines.py) 2153 """ 2154 r0, r1 = self._data.GetScalarRange() 2155 if vmin is None: 2156 vmin = r0 2157 if vmax is None: 2158 vmax = r1 2159 2160 # -------------------------------- 2161 bands = [] 2162 dx = (vmax - vmin) / float(n) 2163 b = [vmin, vmin + dx / 2.0, vmin + dx] 2164 i = 0 2165 while i < n: 2166 bands.append(b) 2167 b = [b[0] + dx, b[1] + dx, b[2] + dx] 2168 i += 1 2169 2170 # annotate, use the midpoint of the band as the label 2171 lut = self.mapper().GetLookupTable() 2172 labels = [] 2173 for b in bands: 2174 labels.append("{:4.2f}".format(b[1])) 2175 values = vtk.vtkVariantArray() 2176 for la in labels: 2177 values.InsertNextValue(vtk.vtkVariant(la)) 2178 for i in range(values.GetNumberOfTuples()): 2179 lut.SetAnnotation(i, values.GetValue(i).ToString()) 2180 2181 bcf = vtk.vtkBandedPolyDataContourFilter() 2182 bcf.SetInputData(self.polydata()) 2183 # Use either the minimum or maximum value for each band. 2184 for i, band in enumerate(bands): 2185 bcf.SetValue(i, band[2]) 2186 # We will use an indexed lookup table. 2187 bcf.SetScalarModeToIndex() 2188 bcf.GenerateContourEdgesOff() 2189 bcf.Update() 2190 bcf.GetOutput().GetCellData().GetScalars().SetName("IsoBands") 2191 m1 = Mesh(bcf.GetOutput()).compute_normals(cells=True) 2192 m1.mapper().SetLookupTable(lut) 2193 2194 m1.pipeline = OperationNode("isobands", parents=[self]) 2195 return m1 2196 2197 def isolines(self, n=10, vmin=None, vmax=None): 2198 """ 2199 Return a new `Mesh` representing the isolines of the active scalars. 2200 2201 Arguments: 2202 n : (int) 2203 number of isolines in the range 2204 vmin : (float) 2205 minimum of the range 2206 vmax : (float) 2207 maximum of the range 2208 2209 Examples: 2210 - [isolines.py](https://github.com/marcomusy/vedo/tree/master/examples/pyplot/isolines.py) 2211 2212  2213 """ 2214 bcf = vtk.vtkContourFilter() 2215 bcf.SetInputData(self.polydata()) 2216 r0, r1 = self._data.GetScalarRange() 2217 if vmin is None: 2218 vmin = r0 2219 if vmax is None: 2220 vmax = r1 2221 bcf.GenerateValues(n, vmin, vmax) 2222 bcf.Update() 2223 sf = vtk.vtkStripper() 2224 sf.SetJoinContiguousSegments(True) 2225 sf.SetInputData(bcf.GetOutput()) 2226 sf.Update() 2227 cl = vtk.vtkCleanPolyData() 2228 cl.SetInputData(sf.GetOutput()) 2229 cl.Update() 2230 msh = Mesh(cl.GetOutput(), c="k").lighting("off") 2231 msh.mapper().SetResolveCoincidentTopologyToPolygonOffset() 2232 msh.pipeline = OperationNode("isolines", parents=[self]) 2233 return msh 2234 2235 def extrude(self, zshift=1, rotation=0, dR=0, cap=True, res=1): 2236 """ 2237 Sweep a polygonal data creating a "skirt" from free edges and lines, and lines from vertices. 2238 The input dataset is swept around the z-axis to create new polygonal primitives. 2239 For example, sweeping a line results in a cylindrical shell, and sweeping a circle creates a torus. 2240 2241 You can control whether the sweep of a 2D object (i.e., polygon or triangle strip) 2242 is capped with the generating geometry. 2243 Also, you can control the angle of rotation, and whether translation along the z-axis 2244 is performed along with the rotation. (Translation is useful for creating "springs"). 2245 You also can adjust the radius of the generating geometry using the "dR" keyword. 2246 2247 The skirt is generated by locating certain topological features. 2248 Free edges (edges of polygons or triangle strips only used by one polygon or triangle strips) 2249 generate surfaces. This is true also of lines or polylines. Vertices generate lines. 2250 2251 This filter can be used to model axisymmetric objects like cylinders, bottles, and wine glasses; 2252 or translational/rotational symmetric objects like springs or corkscrews. 2253 2254 Warning: 2255 Some polygonal objects have no free edges (e.g., sphere). When swept, this will result 2256 in two separate surfaces if capping is on, or no surface if capping is off. 2257 2258 Examples: 2259 - [extrude.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/extrude.py) 2260 2261  2262 """ 2263 if is_sequence(zshift): 2264 # ms = [] # todo 2265 # poly0 = self.clone().polydata() 2266 # for i in range(len(zshift)-1): 2267 # rf = vtk.vtkRotationalExtrusionFilter() 2268 # rf.SetInputData(poly0) 2269 # rf.SetResolution(res) 2270 # rf.SetCapping(0) 2271 # rf.SetAngle(rotation) 2272 # rf.SetTranslation(zshift) 2273 # rf.SetDeltaRadius(dR) 2274 # rf.Update() 2275 # poly1 = rf.GetOutput() 2276 return self 2277 2278 rf = vtk.vtkRotationalExtrusionFilter() 2279 # rf = vtk.vtkLinearExtrusionFilter() 2280 rf.SetInputData(self.polydata(False)) # must not be transformed 2281 rf.SetResolution(res) 2282 rf.SetCapping(cap) 2283 rf.SetAngle(rotation) 2284 rf.SetTranslation(zshift) 2285 rf.SetDeltaRadius(dR) 2286 rf.Update() 2287 m = Mesh(rf.GetOutput(), c=self.c(), alpha=self.alpha()) 2288 prop = vtk.vtkProperty() 2289 prop.DeepCopy(self.property) 2290 m.SetProperty(prop) 2291 m.property = prop 2292 # assign the same transformation 2293 m.SetOrigin(self.GetOrigin()) 2294 m.SetScale(self.GetScale()) 2295 m.SetOrientation(self.GetOrientation()) 2296 m.SetPosition(self.GetPosition()) 2297 2298 m.compute_normals(cells=False).flat().lighting("default") 2299 2300 m.pipeline = OperationNode( 2301 "extrude", parents=[self], 2302 comment=f"#pts {m._data.GetNumberOfPoints()}", 2303 ) 2304 return m 2305 2306 def split(self, maxdepth=1000, flag=False): 2307 """ 2308 Split a mesh by connectivity and order the pieces by increasing area. 2309 2310 Arguments: 2311 maxdepth : (int) 2312 only consider this maximum number of mesh parts. 2313 flag : (bool) 2314 if set to True return the same single object, 2315 but add a "RegionId" array to flag the mesh subparts 2316 2317 Examples: 2318 - [splitmesh.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/splitmesh.py) 2319 2320  2321 """ 2322 pd = self.polydata(False) 2323 cf = vtk.vtkConnectivityFilter() 2324 cf.SetInputData(pd) 2325 cf.SetExtractionModeToAllRegions() 2326 cf.SetColorRegions(True) 2327 cf.Update() 2328 2329 if flag: 2330 return self._update(cf.GetOutput()) 2331 2332 if not cf.GetOutput().GetNumberOfPoints(): 2333 return [] 2334 2335 a = Mesh(cf.GetOutput()) 2336 alist = [] 2337 for t in range(max(a.pointdata["RegionId"]) + 1): 2338 if t == maxdepth: 2339 break 2340 suba = a.clone().threshold("RegionId", t - 0.1, t + 0.1) 2341 area = suba.area() 2342 alist.append([suba, area]) 2343 2344 alist.sort(key=lambda x: x[1]) 2345 alist.reverse() 2346 blist = [] 2347 for i, l in enumerate(alist): 2348 l[0].color(i + 1).phong() 2349 l[0].mapper().ScalarVisibilityOff() 2350 blist.append(l[0]) 2351 2352 if i<10: 2353 l[0].pipeline = OperationNode( 2354 f"split mesh {i}", parents=[self], 2355 comment=f"#pts {l[0]._data.GetNumberOfPoints()}", 2356 ) 2357 return blist 2358 2359 @deprecated(reason=vedo.colors.red + "Please use extract_largest_region()" + vedo.colors.reset) 2360 def extractLargestRegion(self): 2361 "Deprecated. Please use `extract_largest_region()`" 2362 return self.extract_largest_region() 2363 2364 def extract_largest_region(self): 2365 """ 2366 Extract the largest connected part of a mesh and discard all the smaller pieces. 2367 2368 Examples: 2369 - [largestregion.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/largestregion.py) 2370 """ 2371 conn = vtk.vtkConnectivityFilter() 2372 conn.SetExtractionModeToLargestRegion() 2373 conn.ScalarConnectivityOff() 2374 conn.SetInputData(self._data) 2375 conn.Update() 2376 m = Mesh(conn.GetOutput()) 2377 pr = vtk.vtkProperty() 2378 pr.DeepCopy(self.property) 2379 m.SetProperty(pr) 2380 m.property = pr 2381 # assign the same transformation 2382 m.SetOrigin(self.GetOrigin()) 2383 m.SetScale(self.GetScale()) 2384 m.SetOrientation(self.GetOrientation()) 2385 m.SetPosition(self.GetPosition()) 2386 vis = self._mapper.GetScalarVisibility() 2387 m.mapper().SetScalarVisibility(vis) 2388 2389 m.pipeline = OperationNode( 2390 "extract_largest_region", parents=[self], 2391 comment=f"#pts {m._data.GetNumberOfPoints()}", 2392 ) 2393 return m 2394 2395 def boolean(self, operation, mesh2, method=0, tol=None): 2396 """Volumetric union, intersection and subtraction of surfaces. 2397 2398 Use `operation` for the allowed operations `['plus', 'intersect', 'minus']`. 2399 2400 Two possible algorithms are available by changing `method`. 2401 2402 Example: 2403 - [boolean.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/boolean.py) 2404 2405  2406 """ 2407 if method == 0: 2408 bf = vtk.vtkBooleanOperationPolyDataFilter() 2409 elif method == 1: 2410 bf = vtk.vtkLoopBooleanPolyDataFilter() 2411 else: 2412 raise ValueError(f"Unknown method={method}") 2413 2414 poly1 = self.compute_normals().polydata() 2415 poly2 = mesh2.compute_normals().polydata() 2416 2417 if operation.lower() in ("plus", "+"): 2418 bf.SetOperationToUnion() 2419 elif operation.lower() == "intersect": 2420 bf.SetOperationToIntersection() 2421 elif operation.lower() in ("minus", "-"): 2422 bf.SetOperationToDifference() 2423 2424 if tol: 2425 bf.SetTolerance(tol) 2426 2427 bf.SetInputData(0, poly1) 2428 bf.SetInputData(1, poly2) 2429 bf.Update() 2430 2431 msh = Mesh(bf.GetOutput(), c=None) 2432 msh.flat() 2433 msh.name = self.name + operation + mesh2.name 2434 2435 msh.pipeline = OperationNode( 2436 "boolean " + operation, 2437 parents=[self, mesh2], 2438 shape="cylinder", 2439 comment=f"#pts {msh._data.GetNumberOfPoints()}", 2440 ) 2441 return msh 2442 2443 @deprecated(reason=vedo.colors.red + "Please use intersect_with()" + vedo.colors.reset) 2444 def intersectWith(self, mesh2, tol=1e-06): 2445 "Deprecated. Please use `intersect_with()`" 2446 return self.intersect_with(mesh2, tol) 2447 2448 @deprecated(reason=vedo.colors.red + "Please use intersect_with_line()" + vedo.colors.reset) 2449 def intersectWithLine(self, p0, p1=None, returnIds=False, tol=0): 2450 "Deprecated. Please use `intersect_with_line()`" 2451 return self.intersect_with_line(p0, p1, returnIds, tol) 2452 2453 def intersect_with(self, mesh2, tol=1e-06): 2454 """ 2455 Intersect this Mesh with the input surface to return a set of lines. 2456 2457 Examples: 2458 - [surf_intersect.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/surf_intersect.py) 2459 2460  2461 """ 2462 bf = vtk.vtkIntersectionPolyDataFilter() 2463 bf.SetGlobalWarningDisplay(0) 2464 poly1 = self.polydata() 2465 poly2 = mesh2.polydata() 2466 bf.SetTolerance(tol) 2467 bf.SetInputData(0, poly1) 2468 bf.SetInputData(1, poly2) 2469 bf.Update() 2470 msh = Mesh(bf.GetOutput(), "k", 1).lighting("off") 2471 msh.GetProperty().SetLineWidth(3) 2472 msh.name = "SurfaceIntersection" 2473 2474 msh.pipeline = OperationNode( 2475 "intersect_with", parents=[self, mesh2], comment=f"#pts {msh.npoints}", 2476 ) 2477 return msh 2478 2479 def intersect_with_line(self, p0, p1=None, return_ids=False, tol=0): 2480 """ 2481 Return the list of points intersecting the mesh 2482 along the segment defined by two points `p0` and `p1`. 2483 2484 Use `return_ids` to return the cell ids along with point coords 2485 2486 Example: 2487 ```python 2488 from vedo import * 2489 s = Spring() 2490 pts = s.intersect_with_line([0,0,0], [1,0.1,0]) 2491 ln = Line([0,0,0], [1,0.1,0], c='blue') 2492 ps = Points(pts, r=10, c='r') 2493 show(s, ln, ps, bg='white').close() 2494 ``` 2495  2496 """ 2497 if isinstance(p0, Points): 2498 p0, p1 = p0.points() 2499 2500 if not self.line_locator: 2501 self.line_locator = vtk.vtkOBBTree() 2502 self.line_locator.SetDataSet(self.polydata()) 2503 if not tol: 2504 tol = mag(np.asarray(p1) - np.asarray(p0)) / 10000 2505 self.line_locator.SetTolerance(tol) 2506 self.line_locator.BuildLocator() 2507 2508 intersectPoints = vtk.vtkPoints() 2509 idlist = vtk.vtkIdList() 2510 self.line_locator.IntersectWithLine(p0, p1, intersectPoints, idlist) 2511 pts = [] 2512 for i in range(intersectPoints.GetNumberOfPoints()): 2513 intersection = [0, 0, 0] 2514 intersectPoints.GetPoint(i, intersection) 2515 pts.append(intersection) 2516 pts = np.array(pts) 2517 2518 if return_ids: 2519 pts_ids = [] 2520 for i in range(idlist.GetNumberOfIds()): 2521 cid = idlist.GetId(i) 2522 pts_ids.append(cid) 2523 return (pts, np.array(pts_ids).astype(np.uint32)) 2524 2525 return pts 2526 2527 def intersect_with_plane(self, origin=[0,0,0], normal=[1,0,0]): 2528 """ 2529 Intersect this Mesh with a plane to return a set of lines. 2530 2531 Example: 2532 ```python 2533 from vedo import * 2534 sph = Sphere() 2535 mi = sph.clone().intersect_with_plane().join() 2536 print(mi.lines()) 2537 show(sph, mi, axes=1).close() 2538 ``` 2539  2540 """ 2541 plane = vtk.vtkPlane() 2542 plane.SetOrigin(origin) 2543 plane.SetNormal(normal) 2544 2545 cutter = vtk.vtkPolyDataPlaneCutter() 2546 cutter.SetInputData(self.polydata()) 2547 cutter.SetPlane(plane) 2548 cutter.InterpolateAttributesOn() 2549 cutter.ComputeNormalsOff() 2550 cutter.Update() 2551 2552 msh = Mesh(cutter.GetOutput(), "k", 1).lighting("off") 2553 msh.GetProperty().SetLineWidth(3) 2554 msh.name = "PlaneIntersection" 2555 2556 msh.pipeline = OperationNode( 2557 "intersect_with_plan", parents=[self], 2558 comment=f"#pts {msh._data.GetNumberOfPoints()}", 2559 ) 2560 return msh 2561 2562 def intersect_with_multiplanes(self, origin=(0,0,0), normal=(1,0,0), dmin=-1, dmax=1, n=10): 2563 """ 2564 Generate a set of lines from cutting a mesh in n intervals 2565 between a minimum and maximum distance from a plane of given origin and normal. 2566 2567 Arguments: 2568 origin : (list) 2569 the point of the cutting plane 2570 normal : (list) 2571 normal vector to the cutting plane 2572 dmin : (float) 2573 negative distance below the plane 2574 dmax : (float) 2575 positive distance above the plane 2576 n : (int) 2577 number of cuts 2578 """ 2579 plane = vtk.vtkPlane() 2580 poly = self.polydata() 2581 plane.SetOrigin(poly.GetCenter()) 2582 plane.SetNormal(normal) 2583 2584 bounds = np.array(self.bounds()) 2585 min_bound = bounds[[0,2,4]] 2586 max_bound = bounds[[1,3,5]] 2587 2588 center = np.array(poly.GetCenter()) 2589 distance_min = np.linalg.norm(min_bound - center) 2590 distance_max = np.linalg.norm(max_bound - center) 2591 2592 cutter = vtk.vtkCutter() 2593 cutter.SetCutFunction(plane) 2594 cutter.SetInputData(poly) 2595 cutter.GenerateValues(n, -distance_min, distance_max) 2596 cutter.Update() 2597 2598 stripper = vtk.vtkStripper() 2599 stripper.SetJoinContiguousSegments(True) 2600 stripper.SetInputData(cutter.GetOutput()) 2601 stripper.Update() 2602 2603 msh = Mesh(stripper.GetOutput()) 2604 msh.property.LightingOff() 2605 msh.property.SetColor(get_color("k2")) 2606 msh.lw(2) 2607 2608 msh.pipeline = OperationNode( 2609 "intersect_with_multiplanes", parents=[self], 2610 comment=f"#pts {msh._data.GetNumberOfPoints()}", 2611 ) 2612 return msh 2613 2614 def collide_with(self, mesh2, tol=0, return_bool=False): 2615 """ 2616 Collide this Mesh with the input surface. 2617 Information is stored in `ContactCells1` and `ContactCells2`. 2618 """ 2619 ipdf = vtk.vtkCollisionDetectionFilter() 2620 # ipdf.SetGlobalWarningDisplay(0) 2621 2622 transform0 = vtk.vtkTransform() 2623 transform1 = vtk.vtkTransform() 2624 2625 # ipdf.SetBoxTolerance(tol) 2626 ipdf.SetCellTolerance(tol) 2627 ipdf.SetInputData(0, self.polydata()) 2628 ipdf.SetInputData(1, mesh2.polydata()) 2629 ipdf.SetTransform(0, transform0) 2630 ipdf.SetTransform(1, transform1) 2631 if return_bool: 2632 ipdf.SetCollisionModeToFirstContact() 2633 else: 2634 ipdf.SetCollisionModeToAllContacts() 2635 ipdf.Update() 2636 2637 if return_bool: 2638 return bool(ipdf.GetNumberOfContacts()) 2639 2640 msh = Mesh(ipdf.GetContactsOutput(), "k", 1).lighting("off") 2641 msh.metadata["ContactCells1"] = vtk2numpy(ipdf.GetOutput(0).GetFieldData().GetArray("ContactCells")) 2642 msh.metadata["ContactCells2"] = vtk2numpy(ipdf.GetOutput(1).GetFieldData().GetArray("ContactCells")) 2643 msh.GetProperty().SetLineWidth(3) 2644 msh.name = "SurfaceCollision" 2645 2646 msh.pipeline = OperationNode( 2647 "collide_with", parents=[self, mesh2], 2648 comment=f"#pts {msh._data.GetNumberOfPoints()}", 2649 ) 2650 return msh 2651 2652 def geodesic(self, start, end): 2653 """ 2654 Dijkstra algorithm to compute the geodesic line. 2655 Takes as input a polygonal mesh and performs a single source shortest path calculation. 2656 2657 Arguments: 2658 start : (int, list) 2659 start vertex index or close point `[x,y,z]` 2660 end : (int, list) 2661 end vertex index or close point `[x,y,z]` 2662 2663 Examples: 2664 - [geodesic.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/geodesic.py) 2665 2666  2667 """ 2668 if is_sequence(start): 2669 cc = self.points() 2670 pa = Points(cc) 2671 start = pa.closest_point(start, return_point_id=True) 2672 end = pa.closest_point(end, return_point_id=True) 2673 2674 dijkstra = vtk.vtkDijkstraGraphGeodesicPath() 2675 dijkstra.SetInputData(self.polydata()) 2676 dijkstra.SetStartVertex(end) # inverted in vtk 2677 dijkstra.SetEndVertex(start) 2678 dijkstra.Update() 2679 2680 weights = vtk.vtkDoubleArray() 2681 dijkstra.GetCumulativeWeights(weights) 2682 2683 idlist = dijkstra.GetIdList() 2684 ids = [idlist.GetId(i) for i in range(idlist.GetNumberOfIds())] 2685 2686 length = weights.GetMaxId() + 1 2687 arr = np.zeros(length) 2688 for i in range(length): 2689 arr[i] = weights.GetTuple(i)[0] 2690 2691 poly = dijkstra.GetOutput() 2692 2693 vdata = numpy2vtk(arr) 2694 vdata.SetName("CumulativeWeights") 2695 poly.GetPointData().AddArray(vdata) 2696 2697 vdata2 = numpy2vtk(ids, dtype=np.uint) 2698 vdata2.SetName("VertexIDs") 2699 poly.GetPointData().AddArray(vdata2) 2700 poly.GetPointData().Modified() 2701 2702 dmesh = Mesh(poly, c="k") 2703 prop = vtk.vtkProperty() 2704 prop.DeepCopy(self.property) 2705 prop.SetLineWidth(3) 2706 prop.SetOpacity(1) 2707 dmesh.SetProperty(prop) 2708 dmesh.property = prop 2709 dmesh.name = "GeodesicLine" 2710 2711 dmesh.pipeline = OperationNode( 2712 "GeodesicLine", parents=[self], 2713 comment=f"#pts {dmesh._data.GetNumberOfPoints()}", 2714 ) 2715 return dmesh 2716 2717 ##################################################################### 2718 ### Stuff returning a Volume 2719 ### 2720 def binarize(self, spacing=(1, 1, 1), invert=False, direction_matrix=None, 2721 image_size=None, origin=None, fg_val=255, bg_val=0): 2722 """ 2723 Convert a `Mesh` into a `Volume` 2724 where the foreground (exterior) voxels value is fg_val (255 by default) 2725 and the background (interior) voxels value is bg_val (0 by default). 2726 2727 Examples: 2728 - [mesh2volume.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/mesh2volume.py) 2729 2730  2731 """ 2732 # https://vtk.org/Wiki/VTK/Examples/Cxx/PolyData/PolyDataToImageData 2733 pd = self.polydata() 2734 2735 whiteImage = vtk.vtkImageData() 2736 if direction_matrix: 2737 whiteImage.SetDirectionMatrix(direction_matrix) 2738 2739 dim = [0, 0, 0] if not image_size else image_size 2740 2741 bounds = self.bounds() 2742 if not image_size: # compute dimensions 2743 dim = [0, 0, 0] 2744 for i in [0, 1, 2]: 2745 dim[i] = int( 2746 np.ceil((bounds[i * 2 + 1] - bounds[i * 2]) / spacing[i]) 2747 ) 2748 2749 whiteImage.SetDimensions(dim) 2750 whiteImage.SetSpacing(spacing) 2751 whiteImage.SetExtent(0, dim[0] - 1, 0, dim[1] - 1, 0, dim[2] - 1) 2752 2753 if not origin: 2754 origin = [0, 0, 0] 2755 origin[0] = bounds[0] + spacing[0] / 2 2756 origin[1] = bounds[2] + spacing[1] / 2 2757 origin[2] = bounds[4] + spacing[2] / 2 2758 whiteImage.SetOrigin(origin) 2759 whiteImage.AllocateScalars(vtk.VTK_UNSIGNED_CHAR, 1) 2760 2761 # fill the image with foreground voxels: 2762 inval = bg_val if invert else fg_val 2763 whiteImage.GetPointData().GetScalars().Fill(inval) 2764 2765 # polygonal data --> image stencil: 2766 pol2stenc = vtk.vtkPolyDataToImageStencil() 2767 pol2stenc.SetInputData(pd) 2768 pol2stenc.SetOutputOrigin(whiteImage.GetOrigin()) 2769 pol2stenc.SetOutputSpacing(whiteImage.GetSpacing()) 2770 pol2stenc.SetOutputWholeExtent(whiteImage.GetExtent()) 2771 pol2stenc.Update() 2772 2773 # cut the corresponding white image and set the background: 2774 outval = fg_val if invert else bg_val 2775 2776 imgstenc = vtk.vtkImageStencil() 2777 imgstenc.SetInputData(whiteImage) 2778 imgstenc.SetStencilConnection(pol2stenc.GetOutputPort()) 2779 imgstenc.SetReverseStencil(invert) 2780 imgstenc.SetBackgroundValue(outval) 2781 imgstenc.Update() 2782 vol = vedo.Volume(imgstenc.GetOutput()) 2783 2784 vol.pipeline = OperationNode( 2785 "binarize", 2786 parents=[self], 2787 comment=f"dim = {tuple(vol.dimensions())}", 2788 c="#e9c46a:#0096c7" 2789 ) 2790 return vol 2791 2792 2793 @deprecated(reason=vedo.colors.red + "Please use signed_distance()" + vedo.colors.reset) 2794 def signedDistance(self, **k): 2795 "Deprecated. Please use `signed_distance()`" 2796 return self.signed_distance(**k) 2797 2798 def signed_distance(self, bounds=None, dims=(20, 20, 20), invert=False, maxradius=None): 2799 """ 2800 Compute the `Volume` object whose voxels contains the signed distance from 2801 the mesh. 2802 2803 Arguments: 2804 bounds : (list) 2805 bounds of the output volume 2806 dims : (list) 2807 dimensions (nr. of voxels) of the output volume 2808 invert : (bool) 2809 flip the sign 2810 2811 Examples: 2812 - [volumeFromMesh.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/volumeFromMesh.py) 2813 """ 2814 if maxradius is not None: 2815 vedo.logger.warning( 2816 "in signedDistance(maxradius=...) is ignored. (Only valid for pointclouds)." 2817 ) 2818 if bounds is None: 2819 bounds = self.bounds() 2820 sx = (bounds[1] - bounds[0]) / dims[0] 2821 sy = (bounds[3] - bounds[2]) / dims[1] 2822 sz = (bounds[5] - bounds[4]) / dims[2] 2823 2824 img = vtk.vtkImageData() 2825 img.SetDimensions(dims) 2826 img.SetSpacing(sx, sy, sz) 2827 img.SetOrigin(bounds[0], bounds[2], bounds[4]) 2828 img.AllocateScalars(vtk.VTK_FLOAT, 1) 2829 2830 imp = vtk.vtkImplicitPolyDataDistance() 2831 imp.SetInput(self.polydata()) 2832 b2 = bounds[2] 2833 b4 = bounds[4] 2834 d0, d1, d2 = dims 2835 2836 for i in range(d0): 2837 x = i * sx + bounds[0] 2838 for j in range(d1): 2839 y = j * sy + b2 2840 for k in range(d2): 2841 v = imp.EvaluateFunction((x, y, k * sz + b4)) 2842 if invert: 2843 v = -v 2844 img.SetScalarComponentFromFloat(i, j, k, 0, v) 2845 2846 vol = vedo.Volume(img) 2847 vol.name = "SignedVolume" 2848 2849 vol.pipeline = OperationNode( 2850 "signed_distance", parents=[self], 2851 comment=f"dim = {tuple(vol.dimensions())}", 2852 c="#e9c46a:#0096c7" 2853 ) 2854 return vol 2855 2856 def tetralize( 2857 self, 2858 side=0.02, 2859 nmax=300_000, 2860 gap=None, 2861 subsample=False, 2862 uniform=True, 2863 seed=0, 2864 debug=False, 2865 ): 2866 """ 2867 Tetralize a closed polygonal mesh. Return a `TetMesh`. 2868 2869 Arguments: 2870 side : (float) 2871 desired side of the single tetras as fraction of the bounding box diagonal. 2872 Typical values are in the range (0.01 - 0.03) 2873 nmax : (int) 2874 maximum random numbers to be sampled in the bounding box 2875 gap : (float) 2876 keep this minimum distance from the surface, 2877 if None an automatic choice is made. 2878 subsample : (bool) 2879 subsample input surface, the geometry might be affected 2880 (the number of original faces reduceed), but higher tet quality might be obtained. 2881 uniform : (bool) 2882 generate tets more uniformly packed in the interior of the mesh 2883 seed : (int) 2884 random number generator seed 2885 debug : (bool) 2886 show an intermediate plot with sampled points 2887 2888 Examples: 2889 - [tetralize_surface.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/tetralize_surface.py) 2890 2891  2892 """ 2893 surf = self.clone().clean().compute_normals() 2894 d = surf.diagonal_size() 2895 if gap is None: 2896 gap = side * d * np.sqrt(2 / 3) 2897 n = int(min((1 / side) ** 3, nmax)) 2898 2899 # fill the space w/ points 2900 x0, x1, y0, y1, z0, z1 = surf.bounds() 2901 2902 if uniform: 2903 pts = vedo.utils.pack_spheres([x0, x1, y0, y1, z0, z1], side * d * 1.42) 2904 pts += np.random.randn(len(pts), 3) * side * d * 1.42 / 100 # some small jitter 2905 else: 2906 disp = np.array([x0 + x1, y0 + y1, z0 + z1]) / 2 2907 np.random.seed(seed) 2908 pts = (np.random.rand(n, 3) - 0.5) * np.array( 2909 [x1 - x0, y1 - y0, z1 - z0] 2910 ) + disp 2911 2912 normals = surf.celldata["Normals"] 2913 cc = surf.cell_centers() 2914 subpts = cc - normals * gap * 1.05 2915 pts = pts.tolist() + subpts.tolist() 2916 2917 if debug: 2918 print(".. tetralize(): subsampling and cleaning") 2919 2920 fillpts = surf.inside_points(pts) 2921 fillpts.subsample(side) 2922 2923 if gap: 2924 fillpts.distance_to(surf) 2925 fillpts.threshold("Distance", above=gap) 2926 2927 if subsample: 2928 surf.subsample(side) 2929 2930 tmesh = vedo.tetmesh.delaunay3d(vedo.merge(fillpts, surf)) 2931 tcenters = tmesh.cell_centers() 2932 2933 ids = surf.inside_points(tcenters, return_ids=True) 2934 ins = np.zeros(tmesh.ncells) 2935 ins[ids] = 1 2936 2937 if debug: 2938 # vedo.pyplot.histogram(fillpts.pointdata["Distance"], xtitle=f"gap={gap}").show().close() 2939 edges = self.edges() 2940 points = self.points() 2941 elen = mag(points[edges][:, 0, :] - points[edges][:, 1, :]) 2942 histo = vedo.pyplot.histogram(elen, xtitle="edge length", xlim=(0, 3 * side * d)) 2943 print(".. edges min, max", elen.min(), elen.max()) 2944 fillpts.cmap("bone") 2945 vedo.show([ 2946 [ f"This is a debug plot.\n\nGenerated points: {n}\ngap: {gap}", 2947 surf.wireframe().alpha(0.2), 2948 vedo.addons.Axes(surf), 2949 fillpts, 2950 Points(subpts).c("r4").ps(3), 2951 ], 2952 [ f"Edges mean length: {np.mean(elen)}\n\nPress q to continue", histo 2953 ], 2954 ], 2955 N=2, 2956 sharecam=False, 2957 new=True, 2958 ).close() 2959 print(".. thresholding") 2960 2961 tmesh.celldata["inside"] = ins.astype(np.uint8) 2962 tmesh.threshold("inside", above=0.9) 2963 tmesh.celldata.remove("inside") 2964 2965 if debug: 2966 print(f".. tetralize() completed, ntets = {tmesh.ncells}") 2967 2968 tmesh.pipeline = OperationNode( 2969 "tetralize", 2970 parents=[self], 2971 comment=f"#tets = {tmesh.ncells}", 2972 c="#e9c46a:#9e2a2b" 2973 ) 2974 return tmesh 2975 2976#################################################### 2977class Follower(vedo.base.BaseActor, vtk.vtkFollower): 2978 2979 def __init__(self, actor, camera=None): 2980 2981 vtk.vtkFollower.__init__(self) 2982 vedo.base.BaseActor.__init__(self) 2983 2984 self.name = actor.name 2985 self._isfollower = False 2986 2987 self.SetMapper(actor.GetMapper()) 2988 2989 self.SetProperty(actor.GetProperty()) 2990 self.SetBackfaceProperty(actor.GetBackfaceProperty()) 2991 self.SetTexture(actor.GetTexture()) 2992 2993 self.SetCamera(camera) 2994 self.SetOrigin(actor.GetOrigin()) 2995 self.SetScale(actor.GetScale()) 2996 self.SetOrientation(actor.GetOrientation()) 2997 self.SetPosition(actor.GetPosition()) 2998 self.SetUseBounds(actor.GetUseBounds()) 2999 3000 self.PickableOff() 3001 3002 self.pipeline = OperationNode( 3003 "Follower", 3004 parents=[actor], 3005 shape='component', 3006 c="#d9ed92", 3007 )
32class Mesh(Points): 33 """ 34 Build an instance of object `Mesh` derived from `vedo.PointCloud`. 35 """ 36 def __init__( 37 self, 38 inputobj=None, 39 c=None, 40 alpha=1, 41 ): 42 """ 43 Input can be a list of vertices and their connectivity (faces of the polygonal mesh), 44 or directly a `vtkPolydata` object. 45 For point clouds - e.i. no faces - just substitute the `faces` list with `None`. 46 47 Example: 48 `Mesh( [ [[x1,y1,z1],[x2,y2,z2], ...], [[0,1,2], [1,2,3], ...] ] )` 49 50 Arguments: 51 c : (color) 52 color in RGB format, hex, symbol or name 53 alpha : (float) 54 mesh opacity [0,1] 55 56 Examples: 57 - [buildmesh.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/buildmesh.py) 58 (and many others!) 59 60  61 """ 62 Points.__init__(self) 63 64 self.line_locator = None 65 66 self._mapper.SetInterpolateScalarsBeforeMapping( 67 vedo.settings.interpolate_scalars_before_mapping 68 ) 69 70 if vedo.settings.use_polygon_offset: 71 self._mapper.SetResolveCoincidentTopologyToPolygonOffset() 72 pof, pou = ( 73 vedo.settings.polygon_offset_factor, 74 vedo.settings.polygon_offset_units, 75 ) 76 self._mapper.SetResolveCoincidentTopologyPolygonOffsetParameters(pof, pou) 77 78 inputtype = str(type(inputobj)) 79 80 if inputobj is None: 81 pass 82 83 elif isinstance(inputobj, (Mesh, vtk.vtkActor)): 84 polyCopy = vtk.vtkPolyData() 85 polyCopy.DeepCopy(inputobj.GetMapper().GetInput()) 86 self._data = polyCopy 87 self._mapper.SetInputData(polyCopy) 88 self._mapper.SetScalarVisibility(inputobj.GetMapper().GetScalarVisibility()) 89 pr = vtk.vtkProperty() 90 pr.DeepCopy(inputobj.GetProperty()) 91 self.SetProperty(pr) 92 self.property = pr 93 94 elif isinstance(inputobj, vtk.vtkPolyData): 95 if inputobj.GetNumberOfCells() == 0: 96 carr = vtk.vtkCellArray() 97 for i in range(inputobj.GetNumberOfPoints()): 98 carr.InsertNextCell(1) 99 carr.InsertCellPoint(i) 100 inputobj.SetVerts(carr) 101 self._data = inputobj # cache vtkPolyData and mapper for speed 102 103 elif isinstance(inputobj, (vtk.vtkStructuredGrid, vtk.vtkRectilinearGrid)): 104 gf = vtk.vtkGeometryFilter() 105 gf.SetInputData(inputobj) 106 gf.Update() 107 self._data = gf.GetOutput() 108 109 elif "trimesh" in inputtype: 110 tact = vedo.utils.trimesh2vedo(inputobj) 111 self._data = tact.polydata() 112 113 elif "meshio" in inputtype: 114 if len(inputobj.cells): 115 mcells = [] 116 for cellblock in inputobj.cells: 117 if cellblock.type in ("triangle", "quad"): 118 mcells += cellblock.data.tolist() 119 self._data = buildPolyData(inputobj.points, mcells) 120 else: 121 self._data = buildPolyData(inputobj.points, None) 122 # add arrays: 123 try: 124 if len(inputobj.point_data): 125 for k in inputobj.point_data.keys(): 126 vdata = numpy2vtk(inputobj.point_data[k]) 127 vdata.SetName(str(k)) 128 self._data.GetPointData().AddArray(vdata) 129 except AssertionError: 130 print("Could not add meshio point data, skip.") 131 # try: 132 # if len(inputobj.cell_data): 133 # for k in inputobj.cell_data.keys(): 134 # #print(inputobj.cell_data) 135 # exit() 136 # vdata = numpy2vtk(inputobj.cell_data[k]) 137 # vdata.SetName(str(k)) 138 # self._data.GetCellData().AddArray(vdata) 139 # except AssertionError: 140 # print("Could not add meshio cell data, skip.") 141 142 elif "meshlab" in inputtype: 143 self._data = vedo.utils.meshlab2vedo(inputobj) 144 145 elif is_sequence(inputobj): 146 ninp = len(inputobj) 147 if ninp == 0: 148 self._data = vtk.vtkPolyData() 149 elif ninp == 2: # assume [vertices, faces] 150 self._data = buildPolyData(inputobj[0], inputobj[1]) 151 else: # assume [vertices] or vertices 152 self._data = buildPolyData(inputobj, None) 153 154 elif hasattr(inputobj, "GetOutput"): # passing vtk object 155 if hasattr(inputobj, "Update"): 156 inputobj.Update() 157 if isinstance(inputobj.GetOutput(), vtk.vtkPolyData): 158 self._data = inputobj.GetOutput() 159 else: 160 gf = vtk.vtkGeometryFilter() 161 gf.SetInputData(inputobj.GetOutput()) 162 gf.Update() 163 self._data = gf.GetOutput() 164 165 elif isinstance(inputobj, str): 166 dataset = vedo.io.load(inputobj) 167 self.filename = inputobj 168 if "TetMesh" in str(type(dataset)): 169 self._data = dataset.tomesh().polydata(False) 170 else: 171 self._data = dataset.polydata(False) 172 173 else: 174 try: 175 gf = vtk.vtkGeometryFilter() 176 gf.SetInputData(inputobj) 177 gf.Update() 178 self._data = gf.GetOutput() 179 except: 180 vedo.logger.error(f"cannot build mesh from type {inputtype}") 181 raise RuntimeError() 182 183 self._mapper.SetInputData(self._data) 184 185 self.property = self.GetProperty() 186 self.property.SetInterpolationToPhong() 187 188 # set the color by c or by scalar 189 if self._data: 190 191 arrexists = False 192 193 if c is None: 194 ptdata = self._data.GetPointData() 195 cldata = self._data.GetCellData() 196 exclude = ["normals", "tcoord"] 197 198 if cldata.GetNumberOfArrays(): 199 for i in range(cldata.GetNumberOfArrays()): 200 iarr = cldata.GetArray(i) 201 if iarr: 202 icname = iarr.GetName() 203 if icname and all(s not in icname.lower() for s in exclude): 204 cldata.SetActiveScalars(icname) 205 self._mapper.ScalarVisibilityOn() 206 self._mapper.SetScalarModeToUseCellData() 207 self._mapper.SetScalarRange(iarr.GetRange()) 208 arrexists = True 209 break # stop at first good one 210 211 # point come after so it has priority 212 if ptdata.GetNumberOfArrays(): 213 for i in range(ptdata.GetNumberOfArrays()): 214 iarr = ptdata.GetArray(i) 215 if iarr: 216 ipname = iarr.GetName() 217 if ipname and all(s not in ipname.lower() for s in exclude): 218 ptdata.SetActiveScalars(ipname) 219 self._mapper.ScalarVisibilityOn() 220 self._mapper.SetScalarModeToUsePointData() 221 self._mapper.SetScalarRange(iarr.GetRange()) 222 arrexists = True 223 break # stop at first good one 224 225 if not arrexists: 226 if c is None: 227 c = "gold" 228 c = get_color(c) 229 elif isinstance(c, float) and c <= 1: 230 c = color_map(c, "rainbow", 0, 1) 231 else: 232 c = get_color(c) 233 self.property.SetColor(c) 234 self.property.SetAmbient(0.1) 235 self.property.SetDiffuse(1) 236 self.property.SetSpecular(0.05) 237 self.property.SetSpecularPower(5) 238 self._mapper.ScalarVisibilityOff() 239 240 if alpha is not None: 241 self.property.SetOpacity(alpha) 242 243 n = self._data.GetNumberOfPoints() 244 self.pipeline = OperationNode(self, comment=f"#pts {n}") 245 246 def _repr_html_(self): 247 """ 248 HTML representation of the Mesh object for Jupyter Notebooks. 249 250 Returns: 251 HTML text with the image and some properties. 252 """ 253 import io 254 import base64 255 from PIL import Image 256 257 library_name = "vedo.mesh.Mesh" 258 help_url = "https://vedo.embl.es/docs/vedo/mesh.html" 259 260 arr = self.thumbnail() 261 im = Image.fromarray(arr) 262 buffered = io.BytesIO() 263 im.save(buffered, format="PNG", quality=100) 264 encoded = base64.b64encode(buffered.getvalue()).decode("utf-8") 265 url = "data:image/png;base64," + encoded 266 image = f"<img src='{url}'></img>" 267 268 # mesh statisitics 269 bounds = "<br/>".join( 270 [ 271 precision(min_x,4) + " ... " + precision(max_x,4) 272 for min_x, max_x in zip(self.bounds()[::2], self.bounds()[1::2]) 273 ] 274 ) 275 average_size = "{size:.3f}".format(size=self.average_size()) 276 277 help_text = "" 278 if self.name: 279 help_text += f"<b> {self.name}:   </b>" 280 help_text += '<b><a href="' + help_url + '" target="_blank">' + library_name + "</a></b>" 281 if self.filename: 282 dots = "" 283 if len(self.filename) > 30: 284 dots = "..." 285 help_text += f"<br/><code><i>({dots}{self.filename[-30:]})</i></code>" 286 287 pdata = "" 288 if self._data.GetPointData().GetScalars(): 289 if self._data.GetPointData().GetScalars().GetName(): 290 name = self._data.GetPointData().GetScalars().GetName() 291 pdata = "<tr><td><b> point data array </b></td><td>" + name + "</td></tr>" 292 293 cdata = "" 294 if self._data.GetCellData().GetScalars(): 295 if self._data.GetCellData().GetScalars().GetName(): 296 name = self._data.GetCellData().GetScalars().GetName() 297 cdata = "<tr><td><b> cell data array </b></td><td>" + name + "</td></tr>" 298 299 all = [ 300 "<table>", 301 "<tr>", 302 "<td>", image, "</td>", 303 "<td style='text-align: center; vertical-align: center;'><br/>", help_text, 304 "<table>", 305 "<tr><td><b> bounds </b> <br/> (x/y/z) </td><td>" + str(bounds) + "</td></tr>", 306 "<tr><td><b> center of mass </b></td><td>" + precision(self.center_of_mass(),3) + "</td></tr>", 307 "<tr><td><b> average size </b></td><td>" + str(average_size) + "</td></tr>", 308 "<tr><td><b> nr. points / faces </b></td><td>" 309 + str(self.npoints) + " / " + str(self.ncells) + "</td></tr>", 310 pdata, 311 cdata, 312 "</table>", 313 "</table>", 314 ] 315 return "\n".join(all) 316 317 318 def faces(self): 319 """ 320 Get cell polygonal connectivity ids as a python `list`. 321 The output format is: `[[id0 ... idn], [id0 ... idm], etc]`. 322 """ 323 arr1d = vtk2numpy(self._data.GetPolys().GetData()) 324 if arr1d is None: 325 return [] 326 327 # Get cell connettivity ids as a 1D array. vtk format is: 328 # [nids1, id0 ... idn, niids2, id0 ... idm, etc]. 329 if len(arr1d) == 0: 330 arr1d = vtk2numpy(self._data.GetStrips().GetData()) 331 if arr1d is None: 332 return [] 333 334 i = 0 335 conn = [] 336 n = len(arr1d) 337 if n: 338 while True: 339 cell = [arr1d[i + k] for k in range(1, arr1d[i] + 1)] 340 conn.append(cell) 341 i += arr1d[i] + 1 342 if i >= n: 343 break 344 return conn # cannot always make a numpy array of it! 345 346 def cells(self): 347 """Alias for `faces()`.""" 348 return self.faces() 349 350 def lines(self, flat=False): 351 """ 352 Get lines connectivity ids as a numpy array. 353 Default format is `[[id0,id1], [id3,id4], ...]` 354 355 Arguments: 356 flat : (bool) 357 return a 1D numpy array as e.g. [2, 10,20, 3, 10,11,12, 2, 70,80, ...] 358 """ 359 # Get cell connettivity ids as a 1D array. The vtk format is: 360 # [nids1, id0 ... idn, niids2, id0 ... idm, etc]. 361 arr1d = vtk2numpy(self.polydata(False).GetLines().GetData()) 362 363 if arr1d is None: 364 return [] 365 366 if flat: 367 return arr1d 368 369 i = 0 370 conn = [] 371 n = len(arr1d) 372 for _ in range(n): 373 cell = [arr1d[i + k + 1] for k in range(arr1d[i])] 374 conn.append(cell) 375 i += arr1d[i] + 1 376 if i >= n: 377 break 378 379 return conn # cannot always make a numpy array of it! 380 381 def edges(self): 382 """Return an array containing the edges connectivity.""" 383 extractEdges = vtk.vtkExtractEdges() 384 extractEdges.SetInputData(self._data) 385 # eed.UseAllPointsOn() 386 extractEdges.Update() 387 lpoly = extractEdges.GetOutput() 388 389 arr1d = vtk2numpy(lpoly.GetLines().GetData()) 390 # [nids1, id0 ... idn, niids2, id0 ... idm, etc]. 391 392 i = 0 393 conn = [] 394 n = len(arr1d) 395 for _ in range(n): 396 cell = [arr1d[i + k + 1] for k in range(arr1d[i])] 397 conn.append(cell) 398 i += arr1d[i] + 1 399 if i >= n: 400 break 401 return conn # cannot always make a numpy array of it! 402 403 def texture( 404 self, 405 tname, 406 tcoords=None, 407 interpolate=True, 408 repeat=True, 409 edge_clamp=False, 410 scale=None, 411 ushift=None, 412 vshift=None, 413 seam_threshold=None, 414 ): 415 """ 416 Assign a texture to mesh from image file or predefined texture `tname`. 417 If tname is set to `None` texture is disabled. 418 Input tname can also be an array or a `vtkTexture`. 419 420 Arguments: 421 tname : (numpy.array, str, Picture, vtkTexture, None) 422 the input texture to be applied. Can be a numpy array, a path to an image file, 423 a vedo Picture. The None value disables texture. 424 tcoords : (numpy.array, str) 425 this is the (u,v) texture coordinate array. Can also be a string of an existing array 426 in the mesh. 427 interpolate : (bool) 428 turn on/off linear interpolation of the texture map when rendering. 429 repeat : (bool) 430 repeat of the texture when tcoords extend beyond the [0,1] range. 431 edge_clamp : (bool) 432 turn on/off the clamping of the texture map when 433 the texture coords extend beyond the [0,1] range. 434 Only used when repeat is False, and edge clamping is supported by the graphics card. 435 scale : (bool) 436 scale the texture image by this factor 437 ushift : (bool) 438 shift u-coordinates of texture by this amount 439 vshift : (bool) 440 shift v-coordinates of texture by this amount 441 seam_threshold : (float) 442 try to seal seams in texture by collapsing triangles 443 (test values around 1.0, lower values = stronger collapse) 444 445 Examples: 446 - [texturecubes.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/texturecubes.py) 447 448  449 """ 450 pd = self.polydata(False) 451 452 if tname is None: # disable texture 453 pd.GetPointData().SetTCoords(None) 454 pd.GetPointData().Modified() 455 return self ###################################### 456 457 if isinstance(tname, vtk.vtkTexture): 458 tu = tname 459 outimg = self._data 460 461 elif isinstance(tname, vedo.Picture): 462 tu = vtk.vtkTexture() 463 outimg = tname.inputdata() 464 465 elif is_sequence(tname): 466 tu = vtk.vtkTexture() 467 outimg = vedo.picture._get_img(tname) 468 469 elif isinstance(tname, str): 470 tu = vtk.vtkTexture() 471 472 if "https://" in tname: 473 try: 474 tname = vedo.io.download(tname, verbose=False) 475 except: 476 vedo.logger.error(f"texture {tname} could not be downloaded") 477 return self 478 479 fn = tname + ".jpg" 480 if os.path.exists(tname): 481 fn = tname 482 else: 483 vedo.logger.error(f"texture file {tname} does not exist") 484 return self 485 486 fnl = fn.lower() 487 if ".jpg" in fnl or ".jpeg" in fnl: 488 reader = vtk.vtkJPEGReader() 489 elif ".png" in fnl: 490 reader = vtk.vtkPNGReader() 491 elif ".bmp" in fnl: 492 reader = vtk.vtkBMPReader() 493 else: 494 vedo.logger.error( 495 "in texture() supported files are only PNG, BMP or JPG" 496 ) 497 return self 498 reader.SetFileName(fn) 499 reader.Update() 500 outimg = reader.GetOutput() 501 502 else: 503 vedo.logger.error(f"in texture() cannot understand input {type(tname)}") 504 return self 505 506 if tcoords is not None: 507 508 if isinstance(tcoords, str): 509 510 vtarr = pd.GetPointData().GetArray(tcoords) 511 512 else: 513 514 tcoords = np.asarray(tcoords) 515 if tcoords.ndim != 2: 516 vedo.logger.error("tcoords must be a 2-dimensional array") 517 return self 518 if tcoords.shape[0] != pd.GetNumberOfPoints(): 519 vedo.logger.error("nr of texture coords must match nr of points") 520 return self 521 if tcoords.shape[1] != 2: 522 vedo.logger.error("tcoords texture vector must have 2 components") 523 vtarr = numpy2vtk(tcoords) 524 vtarr.SetName("TCoordinates") 525 526 pd.GetPointData().SetTCoords(vtarr) 527 pd.GetPointData().Modified() 528 529 elif not pd.GetPointData().GetTCoords(): 530 531 # TCoords still void.. 532 # check that there are no texture-like arrays: 533 names = self.pointdata.keys() 534 candidate_arr = "" 535 for name in names: 536 vtarr = pd.GetPointData().GetArray(name) 537 if vtarr.GetNumberOfComponents() != 2: 538 continue 539 t0, t1 = vtarr.GetRange() 540 if t0 >= 0 and t1 <= 1: 541 candidate_arr = name 542 543 if candidate_arr: 544 545 vtarr = pd.GetPointData().GetArray(candidate_arr) 546 pd.GetPointData().SetTCoords(vtarr) 547 pd.GetPointData().Modified() 548 549 else: 550 # last resource is automatic mapping 551 tmapper = vtk.vtkTextureMapToPlane() 552 tmapper.AutomaticPlaneGenerationOn() 553 tmapper.SetInputData(pd) 554 tmapper.Update() 555 tc = tmapper.GetOutput().GetPointData().GetTCoords() 556 if scale or ushift or vshift: 557 ntc = vtk2numpy(tc) 558 if scale: 559 ntc *= scale 560 if ushift: 561 ntc[:, 0] += ushift 562 if vshift: 563 ntc[:, 1] += vshift 564 tc = numpy2vtk(tc) 565 pd.GetPointData().SetTCoords(tc) 566 pd.GetPointData().Modified() 567 568 tu.SetInputData(outimg) 569 tu.SetInterpolate(interpolate) 570 tu.SetRepeat(repeat) 571 tu.SetEdgeClamp(edge_clamp) 572 573 self.property.SetColor(1, 1, 1) 574 self._mapper.ScalarVisibilityOff() 575 self.SetTexture(tu) 576 577 if seam_threshold is not None: 578 tname = self._data.GetPointData().GetTCoords().GetName() 579 grad = self.gradient(tname) 580 ugrad, vgrad = np.split(grad, 2, axis=1) 581 ugradm, vgradm = vedo.utils.mag2(ugrad), vedo.utils.mag2(vgrad) 582 gradm = np.log(ugradm + vgradm) 583 largegrad_ids = np.arange(len(grad))[gradm > seam_threshold * 4] 584 uvmap = self.pointdata[tname] 585 # collapse triangles that have large gradient 586 new_points = self.points(transformed=False) 587 for f in self.faces(): 588 if np.isin(f, largegrad_ids).all(): 589 id1, id2, id3 = f 590 uv1, uv2, uv3 = uvmap[f] 591 d12 = vedo.mag2(uv1 - uv2) 592 d23 = vedo.mag2(uv2 - uv3) 593 d31 = vedo.mag2(uv3 - uv1) 594 idm = np.argmin([d12, d23, d31]) 595 if idm == 0: 596 new_points[id1] = new_points[id3] 597 new_points[id2] = new_points[id3] 598 elif idm == 1: 599 new_points[id2] = new_points[id1] 600 new_points[id3] = new_points[id1] 601 self.points(new_points) 602 603 self.Modified() 604 return self 605 606 @deprecated(reason=vedo.colors.red + "Please use compute_normals()" + vedo.colors.reset) 607 def computeNormals(self, points=True, cells=True, featureAngle=None, consistency=True): 608 "Deprecated. Please use compute_normals()" 609 return self.compute_normals(points, cells, featureAngle, consistency) 610 611 def compute_normals(self, points=True, cells=True, feature_angle=None, consistency=True): 612 """ 613 Compute cell and vertex normals for the mesh. 614 615 Arguments: 616 points : (bool) 617 do the computation for the vertices too 618 cells : (bool) 619 do the computation for the cells too 620 feature_angle : (float) 621 specify the angle that defines a sharp edge. 622 If the difference in angle across neighboring polygons is greater than this value, 623 the shared edge is considered "sharp" and it is split. 624 consistency : (bool) 625 turn on/off the enforcement of consistent polygon ordering. 626 627 .. warning:: 628 If feature_angle is set to a float the Mesh can be modified, and it 629 can have a different nr. of vertices from the original. 630 """ 631 poly = self.polydata(False) 632 pdnorm = vtk.vtkPolyDataNormals() 633 pdnorm.SetInputData(poly) 634 pdnorm.SetComputePointNormals(points) 635 pdnorm.SetComputeCellNormals(cells) 636 pdnorm.SetConsistency(consistency) 637 pdnorm.FlipNormalsOff() 638 if feature_angle: 639 pdnorm.SetSplitting(True) 640 pdnorm.SetFeatureAngle(feature_angle) 641 else: 642 pdnorm.SetSplitting(False) 643 # print(pdnorm.GetNonManifoldTraversal()) 644 pdnorm.Update() 645 return self._update(pdnorm.GetOutput()) 646 647 def reverse(self, cells=True, normals=False): 648 """ 649 Reverse the order of polygonal cells 650 and/or reverse the direction of point and cell normals. 651 Two flags are used to control these operations: 652 653 - `cells=True` reverses the order of the indices in the cell connectivity list. 654 If cell is a list of IDs only those cells will be reversed. 655 656 - `normals=True` reverses the normals by multiplying the normal vector by -1 657 (both point and cell normals, if present). 658 """ 659 poly = self.polydata(False) 660 661 if is_sequence(cells): 662 for cell in cells: 663 poly.ReverseCell(cell) 664 poly.GetCellData().Modified() 665 return self ############## 666 667 rev = vtk.vtkReverseSense() 668 if cells: 669 rev.ReverseCellsOn() 670 else: 671 rev.ReverseCellsOff() 672 if normals: 673 rev.ReverseNormalsOn() 674 else: 675 rev.ReverseNormalsOff() 676 rev.SetInputData(poly) 677 rev.Update() 678 out = self._update(rev.GetOutput()) 679 out.pipeline = OperationNode("reverse", parents=[self]) 680 return out 681 682 def wireframe(self, value=True): 683 """Set mesh's representation as wireframe or solid surface.""" 684 if value: 685 self.property.SetRepresentationToWireframe() 686 else: 687 self.property.SetRepresentationToSurface() 688 return self 689 690 def flat(self): 691 """Set surface interpolation to flat. 692 693 <img src="https://upload.wikimedia.org/wikipedia/commons/6/6b/Phong_components_version_4.png" width="700"> 694 """ 695 self.property.SetInterpolationToFlat() 696 return self 697 698 def phong(self): 699 """Set surface interpolation to "phong".""" 700 self.property.SetInterpolationToPhong() 701 return self 702 703 def backface_culling(self, value=True): 704 """Set culling of polygons based on orientation of normal with respect to camera.""" 705 self.property.SetBackfaceCulling(value) 706 return self 707 708 def render_lines_as_tubes(self, value=True): 709 """Wrap a fake tube around a simple line for visualization""" 710 self.property.SetRenderLinesAsTubes(value) 711 return self 712 713 def frontface_culling(self, value=True): 714 """Set culling of polygons based on orientation of normal with respect to camera.""" 715 self.property.SetFrontfaceCulling(value) 716 return self 717 718 @deprecated(reason=vedo.colors.red + "Please use backcolor()" + vedo.colors.reset) 719 def backColor(self, bc=None): 720 "Deprecated. Please use `backcolor()`" 721 return self.backcolor(bc) 722 723 def backcolor(self, bc=None): 724 """ 725 Set/get mesh's backface color. 726 """ 727 backProp = self.GetBackfaceProperty() 728 729 if bc is None: 730 if backProp: 731 return backProp.GetDiffuseColor() 732 return self 733 734 if self.property.GetOpacity() < 1: 735 return self 736 737 if not backProp: 738 backProp = vtk.vtkProperty() 739 740 backProp.SetDiffuseColor(get_color(bc)) 741 backProp.SetOpacity(self.property.GetOpacity()) 742 self.SetBackfaceProperty(backProp) 743 self._mapper.ScalarVisibilityOff() 744 return self 745 746 def bc(self, backcolor=False): 747 """Shortcut for `mesh.backcolor()`.""" 748 return self.backcolor(backcolor) 749 750 @deprecated(reason=vedo.colors.red + "Please use linewidth()" + vedo.colors.reset) 751 def lineWidth(self, lw): 752 "Deprecated: please use `linewidth()`" 753 return self.linewidth(lw) 754 755 def linewidth(self, lw=None): 756 """Set/get width of mesh edges. Same as `lw()`.""" 757 if lw is not None: 758 if lw == 0: 759 self.property.EdgeVisibilityOff() 760 self.property.SetRepresentationToSurface() 761 return self 762 self.property.EdgeVisibilityOn() 763 self.property.SetLineWidth(lw) 764 else: 765 return self.property.GetLineWidth() 766 return self 767 768 def lw(self, linewidth=None): 769 """Set/get width of mesh edges. Same as `linewidth()`.""" 770 return self.linewidth(linewidth) 771 772 def linecolor(self, lc=None): 773 """Set/get color of mesh edges. Same as `lc()`.""" 774 if lc is None: 775 return self.property.GetEdgeColor() 776 self.property.EdgeVisibilityOn() 777 self.property.SetEdgeColor(get_color(lc)) 778 return self 779 780 def lc(self, linecolor=None): 781 """Set/get color of mesh edges. Same as `linecolor()`.""" 782 return self.linecolor(linecolor) 783 784 def volume(self): 785 """Get/set the volume occupied by mesh.""" 786 mass = vtk.vtkMassProperties() 787 mass.SetGlobalWarningDisplay(0) 788 mass.SetInputData(self.polydata()) 789 mass.Update() 790 return mass.GetVolume() 791 792 def area(self): 793 """Get/set the surface area of mesh.""" 794 mass = vtk.vtkMassProperties() 795 mass.SetGlobalWarningDisplay(0) 796 mass.SetInputData(self.polydata()) 797 mass.Update() 798 return mass.GetSurfaceArea() 799 800 def is_closed(self): 801 """Return `True` if the mesh is watertight.""" 802 fe = vtk.vtkFeatureEdges() 803 fe.BoundaryEdgesOn() 804 fe.FeatureEdgesOff() 805 fe.NonManifoldEdgesOn() 806 fe.SetInputData(self.polydata(False)) 807 fe.Update() 808 ne = fe.GetOutput().GetNumberOfCells() 809 return not bool(ne) 810 811 def is_manifold(self): 812 """Return `True` if the mesh is manifold.""" 813 fe = vtk.vtkFeatureEdges() 814 fe.BoundaryEdgesOff() 815 fe.FeatureEdgesOff() 816 fe.NonManifoldEdgesOn() 817 fe.SetInputData(self.polydata(False)) 818 fe.Update() 819 ne = fe.GetOutput().GetNumberOfCells() 820 return not bool(ne) 821 822 def non_manifold_faces(self, remove=True, tol="auto"): 823 """ 824 Detect and (try to) remove non-manifold faces of a triangular mesh. 825 826 Set `remove` to `False` to mark cells without removing them. 827 Set `tol=0` for zero-tolerance, the result will be manifold but with holes. 828 Set `tol>0` to cut off non-manifold faces, and try to recover the good ones. 829 Set `tol="auto"` to make an automatic choice of the tolerance. 830 """ 831 # mark original point and cell ids 832 self.add_ids() 833 toremove = self.boundaries( 834 boundary_edges=False, 835 non_manifold_edges=True, 836 cell_edge=True, 837 return_cell_ids=True, 838 ) 839 if len(toremove) == 0: 840 return self 841 842 points = self.points() 843 faces = self.faces() 844 centers = self.cell_centers() 845 846 copy = self.clone() 847 copy.delete_cells(toremove).clean() 848 copy.compute_normals(cells=False) 849 normals = copy.normals() 850 deltas, deltas_i = [], [] 851 852 for i in vedo.utils.progressbar(toremove, delay=3, title="recover faces"): 853 pids = copy.closest_point(centers[i], n=3, return_point_id=True) 854 norms = normals[pids] 855 n = np.mean(norms, axis=0) 856 dn = np.linalg.norm(n) 857 if not dn: 858 continue 859 n = n / dn 860 861 p0, p1, p2 = points[faces[i]][:3] 862 v = np.cross(p1-p0, p2-p0) 863 lv = np.linalg.norm(v) 864 if not lv: 865 continue 866 v = v / lv 867 868 cosa = 1 - np.dot(n, v) 869 deltas.append(cosa) 870 deltas_i.append(i) 871 872 recover = [] 873 if len(deltas): 874 mean_delta = np.mean(deltas) 875 err_delta = np.std(deltas) 876 txt = "" 877 if tol == "auto": #automatic choice 878 tol = mean_delta / 5 879 txt = f"\n Automatic tol. : {tol: .4f}" 880 for i, cosa in zip(deltas_i, deltas): 881 if cosa < tol: 882 recover.append(i) 883 884 vedo.logger.info( 885 f"\n --------- Non manifold faces ---------" 886 f"\n Average tol. : {mean_delta: .4f} +- {err_delta: .4f}{txt}" 887 f"\n Removed faces : {len(toremove)}" 888 f"\n Recovered faces: {len(recover)}" 889 ) 890 891 toremove = list(set(toremove) - set(recover)) 892 893 if not remove: 894 mark = np.zeros(self.ncells, dtype=np.uint8) 895 mark[recover] = 1 896 mark[toremove] = 2 897 self.celldata["NonManifoldCell"] = mark 898 else: 899 self.delete_cells(toremove) 900 901 self.pipeline = OperationNode( 902 "non_manifold_faces", parents=[self], 903 comment=f"#cells {self._data.GetNumberOfCells()}", 904 ) 905 return self 906 907 def shrink(self, fraction=0.85): 908 """Shrink the triangle polydata in the representation of the input mesh. 909 910 Examples: 911 - [shrink.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/shrink.py) 912 913  914 """ 915 shrink = vtk.vtkShrinkPolyData() 916 shrink.SetInputData(self._data) 917 shrink.SetShrinkFactor(fraction) 918 shrink.Update() 919 self.point_locator = None 920 self.cell_locator = None 921 out = self._update(shrink.GetOutput()) 922 923 out.pipeline = OperationNode("shrink", parents=[self]) 924 return out 925 926 def stretch(self, q1, q2): 927 """ 928 Stretch mesh between points `q1` and `q2`. 929 930 Examples: 931 - [aspring.py](https://github.com/marcomusy/vedo/tree/master/examples/simulations/aspring.py) 932 933  934 935 .. note:: 936 for `Mesh` objects, two vectors `mesh.base`, and `mesh.top` must be defined. 937 """ 938 if self.base is None: 939 vedo.logger.error( 940 "in stretch() must define vectors mesh.base and mesh.top at creation" 941 ) 942 raise RuntimeError() 943 944 p1, p2 = self.base, self.top 945 q1, q2, z = np.asarray(q1), np.asarray(q2), np.array([0, 0, 1]) 946 a = p2 - p1 947 b = q2 - q1 948 plength = np.linalg.norm(a) 949 qlength = np.linalg.norm(b) 950 T = vtk.vtkTransform() 951 T.PostMultiply() 952 T.Translate(-p1) 953 cosa = np.dot(a, z) / plength 954 n = np.cross(a, z) 955 if np.linalg.norm(n): 956 T.RotateWXYZ(np.rad2deg(np.arccos(cosa)), n) 957 T.Scale(1, 1, qlength / plength) 958 959 cosa = np.dot(b, z) / qlength 960 n = np.cross(b, z) 961 if np.linalg.norm(n): 962 T.RotateWXYZ(-np.rad2deg(np.arccos(cosa)), n) 963 else: 964 if np.dot(b, z) < 0: 965 T.RotateWXYZ(180, [1,0,0]) 966 967 T.Translate(q1) 968 969 self.SetUserMatrix(T.GetMatrix()) 970 return self 971 972 def crop( 973 self, 974 top=None, 975 bottom=None, 976 right=None, 977 left=None, 978 front=None, 979 back=None, 980 ): 981 """ 982 Crop an `Mesh` object. 983 Use this method at creation (before moving the object). 984 985 Arguments: 986 top : (float) 987 fraction to crop from the top plane (positive z) 988 bottom : (float) 989 fraction to crop from the bottom plane (negative z) 990 front : (float) 991 fraction to crop from the front plane (positive y) 992 back : (float) 993 fraction to crop from the back plane (negative y) 994 right : (float) 995 fraction to crop from the right plane (positive x) 996 left : (float) 997 fraction to crop from the left plane (negative x) 998 999 Example: 1000 ```python 1001 from vedo import Sphere 1002 Sphere().crop(right=0.3, left=0.1).show() 1003 ``` 1004  1005 """ 1006 cu = vtk.vtkBox() 1007 pos = np.array(self.GetPosition()) 1008 x0, x1, y0, y1, z0, z1 = self.bounds() 1009 x0, y0, z0 = [x0, y0, z0] - pos 1010 x1, y1, z1 = [x1, y1, z1] - pos 1011 1012 dx, dy, dz = x1 - x0, y1 - y0, z1 - z0 1013 if top: 1014 z1 = z1 - top * dz 1015 if bottom: 1016 z0 = z0 + bottom * dz 1017 if front: 1018 y1 = y1 - front * dy 1019 if back: 1020 y0 = y0 + back * dy 1021 if right: 1022 x1 = x1 - right * dx 1023 if left: 1024 x0 = x0 + left * dx 1025 bounds = (x0, x1, y0, y1, z0, z1) 1026 1027 cu.SetBounds(bounds) 1028 1029 clipper = vtk.vtkClipPolyData() 1030 clipper.SetInputData(self._data) 1031 clipper.SetClipFunction(cu) 1032 clipper.InsideOutOn() 1033 clipper.GenerateClippedOutputOff() 1034 clipper.GenerateClipScalarsOff() 1035 clipper.SetValue(0) 1036 clipper.Update() 1037 self._update(clipper.GetOutput()) 1038 self.point_locator = None 1039 1040 self.pipeline = OperationNode( 1041 "crop", parents=[self], 1042 comment=f"#pts {self._data.GetNumberOfPoints()}", 1043 ) 1044 return self 1045 1046 1047 def cap(self, return_cap=False): 1048 """ 1049 Generate a "cap" on a clipped mesh, or caps sharp edges. 1050 1051 Examples: 1052 - [cut_and_cap.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/cut_and_cap.py) 1053 1054  1055 1056 See also: `join()`, `join_segments()`, `slice()`. 1057 """ 1058 poly = self._data 1059 1060 fe = vtk.vtkFeatureEdges() 1061 fe.SetInputData(poly) 1062 fe.BoundaryEdgesOn() 1063 fe.FeatureEdgesOff() 1064 fe.NonManifoldEdgesOff() 1065 fe.ManifoldEdgesOff() 1066 fe.Update() 1067 1068 stripper = vtk.vtkStripper() 1069 stripper.SetInputData(fe.GetOutput()) 1070 stripper.JoinContiguousSegmentsOn() 1071 stripper.Update() 1072 1073 boundaryPoly = vtk.vtkPolyData() 1074 boundaryPoly.SetPoints(stripper.GetOutput().GetPoints()) 1075 boundaryPoly.SetPolys(stripper.GetOutput().GetLines()) 1076 1077 rev = vtk.vtkReverseSense() 1078 rev.ReverseCellsOn() 1079 rev.SetInputData(boundaryPoly) 1080 rev.Update() 1081 1082 tf = vtk.vtkTriangleFilter() 1083 tf.SetInputData(rev.GetOutput()) 1084 tf.Update() 1085 1086 if return_cap: 1087 m = Mesh(tf.GetOutput()) 1088 # assign the same transformation to the copy 1089 m.SetOrigin(self.GetOrigin()) 1090 m.SetScale(self.GetScale()) 1091 m.SetOrientation(self.GetOrientation()) 1092 m.SetPosition(self.GetPosition()) 1093 1094 m.pipeline = OperationNode( 1095 "cap", parents=[self], 1096 comment=f"#pts {m._data.GetNumberOfPoints()}", 1097 ) 1098 return m 1099 1100 polyapp = vtk.vtkAppendPolyData() 1101 polyapp.AddInputData(poly) 1102 polyapp.AddInputData(tf.GetOutput()) 1103 polyapp.Update() 1104 out = self._update(polyapp.GetOutput()).clean() 1105 1106 out.pipeline = OperationNode( 1107 "capped", parents=[self], 1108 comment=f"#pts {out._data.GetNumberOfPoints()}", 1109 ) 1110 return out 1111 1112 def join(self, polys=True, reset=False): 1113 """ 1114 Generate triangle strips and/or polylines from 1115 input polygons, triangle strips, and lines. 1116 1117 Input polygons are assembled into triangle strips only if they are triangles; 1118 other types of polygons are passed through to the output and not stripped. 1119 Use mesh.triangulate() to triangulate non-triangular polygons prior to running 1120 this filter if you need to strip all the data. 1121 1122 Also note that if triangle strips or polylines are present in the input 1123 they are passed through and not joined nor extended. 1124 If you wish to strip these use mesh.triangulate() to fragment the input 1125 into triangles and lines prior to applying join(). 1126 1127 Arguments: 1128 polys : (bool) 1129 polygonal segments will be joined if they are contiguous 1130 reset : (bool) 1131 reset points ordering 1132 1133 Warning: 1134 If triangle strips or polylines exist in the input data 1135 they will be passed through to the output data. 1136 This filter will only construct triangle strips if triangle polygons 1137 are available; and will only construct polylines if lines are available. 1138 1139 Example: 1140 ```python 1141 from vedo import * 1142 c1 = Cylinder(pos=(0,0,0), r=2, height=3, axis=(1,.0,0), alpha=.1).triangulate() 1143 c2 = Cylinder(pos=(0,0,2), r=1, height=2, axis=(0,.3,1), alpha=.1).triangulate() 1144 intersect = c1.intersect_with(c2).join(reset=True) 1145 spline = Spline(intersect).c('blue').lw(5) 1146 show(c1, c2, spline, intersect.labels('id'), axes=1).close() 1147 ``` 1148  1149 """ 1150 sf = vtk.vtkStripper() 1151 sf.SetPassThroughCellIds(True) 1152 sf.SetPassThroughPointIds(True) 1153 sf.SetJoinContiguousSegments(polys) 1154 sf.SetInputData(self.polydata(False)) 1155 sf.Update() 1156 if reset: 1157 poly = sf.GetOutput() 1158 cpd = vtk.vtkCleanPolyData() 1159 cpd.PointMergingOn() 1160 cpd.ConvertLinesToPointsOn() 1161 cpd.ConvertPolysToLinesOn() 1162 cpd.ConvertStripsToPolysOn() 1163 cpd.SetInputData(poly) 1164 cpd.Update() 1165 poly = cpd.GetOutput() 1166 vpts = poly.GetCell(0).GetPoints().GetData() 1167 poly.GetPoints().SetData(vpts) 1168 return self._update(poly) 1169 1170 out = self._update(sf.GetOutput()) 1171 out.pipeline = OperationNode("join", parents=[self], 1172 comment=f"#pts {out._data.GetNumberOfPoints()}", 1173 ) 1174 return out 1175 1176 def join_segments(self, closed=True, tol=1e-03): 1177 """ 1178 Join line segments into contiguous lines. 1179 Useful to call with `triangulate()` method. 1180 1181 Returns: 1182 list of `shapes.Lines` 1183 1184 Example: 1185 ```python 1186 from vedo import * 1187 msh = Torus().alpha(0.1).wireframe() 1188 intersection = msh.intersect_with_plane(normal=[1,1,1]).c('purple5') 1189 slices = [s.triangulate() for s in intersection.join_segments()] 1190 show(msh, intersection, merge(slices), axes=1, viewup='z') 1191 ``` 1192  1193 """ 1194 vlines = [] 1195 for ipiece, outline in enumerate(self.split()): 1196 1197 outline.clean() 1198 pts = outline.points() 1199 if len(pts) < 3: 1200 continue 1201 avesize = outline.average_size() 1202 lines = outline.lines() 1203 # print("---lines", lines, "in piece", ipiece) 1204 tol = avesize / pts.shape[0] * tol 1205 1206 k = 0 1207 joinedpts = [pts[k]] 1208 for _ in range(len(pts)): 1209 pk = pts[k] 1210 for j, line in enumerate(lines): 1211 1212 id0, id1 = line[0], line[-1] 1213 p0, p1 = pts[id0], pts[id1] 1214 1215 if np.linalg.norm(p0 - pk) < tol: 1216 n = len(line) 1217 for m in range(1, n): 1218 joinedpts.append(pts[line[m]]) 1219 # joinedpts.append(p1) 1220 k = id1 1221 lines.pop(j) 1222 break 1223 1224 elif np.linalg.norm(p1 - pk) < tol: 1225 n = len(line) 1226 for m in reversed(range(0, n-1)): 1227 joinedpts.append(pts[line[m]]) 1228 # joinedpts.append(p0) 1229 k = id0 1230 lines.pop(j) 1231 break 1232 1233 if len(joinedpts) > 1: 1234 newline = vedo.shapes.Line(joinedpts, closed=closed) 1235 newline.clean() 1236 newline.SetProperty(self.GetProperty()) 1237 newline.property = self.GetProperty() 1238 newline.pipeline = OperationNode("join_segments", parents=[self], 1239 comment=f"#pts {newline._data.GetNumberOfPoints()}", 1240 ) 1241 vlines.append(newline) 1242 1243 return vlines 1244 1245 def slice(self, origin=[0,0,0], normal=[1,0,0]): 1246 """ 1247 Slice a mesh with a plane and fill the contour. 1248 1249 Example: 1250 ```python 1251 from vedo import * 1252 msh = Mesh(dataurl+"bunny.obj").alpha(0.1).wireframe() 1253 mslice = msh.slice(normal=[0,1,0.3], origin=[0,0.16,0]) 1254 mslice.c('purple5') 1255 show(msh, mslice, axes=1) 1256 ``` 1257  1258 1259 See also: `join()`, `join_segments()`, `cap()`, `cut_with_plane()`. 1260 """ 1261 intersection = self.intersect_with_plane(origin=origin, normal=normal) 1262 slices = [s.triangulate() for s in intersection.join_segments()] 1263 mslices = vedo.pointcloud.merge(slices) 1264 if mslices: 1265 mslices.name = "MeshSlice" 1266 mslices.pipeline = OperationNode( 1267 "slice", parents=[self], comment=f"normal = {normal}") 1268 return mslices 1269 1270 1271 def triangulate(self, verts=True, lines=True): 1272 """ 1273 Converts mesh polygons into triangles. 1274 1275 If the input mesh is only made of 2D lines (no faces) the output will be a triangulation 1276 that fills the internal area. The contours may be concave, and may even contain holes, 1277 i.e. a contour may contain an internal contour winding in the opposite 1278 direction to indicate that it is a hole. 1279 1280 Arguments: 1281 verts : (bool) 1282 if True, break input vertex cells into individual vertex cells 1283 (one point per cell). If False, the input vertex cells will be ignored. 1284 lines : (bool) 1285 if True, break input polylines into line segments. 1286 If False, input lines will be ignored and the output will have no lines. 1287 """ 1288 if self._data.GetNumberOfPolys() or self._data.GetNumberOfStrips(): 1289 tf = vtk.vtkTriangleFilter() 1290 tf.SetPassLines(lines) 1291 tf.SetPassVerts(verts) 1292 1293 elif self._data.GetNumberOfLines(): 1294 tf = vtk.vtkContourTriangulator() 1295 1296 else: 1297 vedo.logger.debug("input in triangulate() seems to be void! Skip.") 1298 return self 1299 1300 tf.SetInputData(self._data) 1301 tf.Update() 1302 out = self._update(tf.GetOutput()) 1303 1304 out.pipeline = OperationNode( 1305 "triangulate", parents=[self], comment=f"#cells {out._data.GetNumberOfCells()}", 1306 ) 1307 return out 1308 1309 1310 def compute_cell_area(self, name="Area"): 1311 """Add to this mesh a cell data array containing the areas of the polygonal faces""" 1312 csf = vtk.vtkCellSizeFilter() 1313 csf.SetInputData(self.polydata(False)) 1314 csf.SetComputeArea(True) 1315 csf.SetComputeVolume(False) 1316 csf.SetComputeLength(False) 1317 csf.SetComputeVertexCount(False) 1318 csf.SetAreaArrayName(name) 1319 csf.Update() 1320 return self._update(csf.GetOutput()) 1321 1322 def compute_cell_vertex_count(self, name="VertexCount"): 1323 """Add to this mesh a cell data array containing the nr of vertices 1324 that a polygonal face has.""" 1325 csf = vtk.vtkCellSizeFilter() 1326 csf.SetInputData(self.polydata(False)) 1327 csf.SetComputeArea(False) 1328 csf.SetComputeVolume(False) 1329 csf.SetComputeLength(False) 1330 csf.SetComputeVertexCount(True) 1331 csf.SetVertexCountArrayName(name) 1332 csf.Update() 1333 return self._update(csf.GetOutput()) 1334 1335 def compute_quality(self, metric=6): 1336 """ 1337 Calculate metrics of quality for the elements of a triangular mesh. 1338 This method adds to the mesh a cell array named "Quality". 1339 See class [vtkMeshQuality](https://vtk.org/doc/nightly/html/classvtkMeshQuality.html) 1340 for explanation. 1341 1342 Arguments: 1343 metric : (int) 1344 type of available estimators are: 1345 - EDGE RATIO, 0 1346 - ASPECT RATIO, 1 1347 - RADIUS RATIO, 2 1348 - ASPECT FROBENIUS, 3 1349 - MED ASPECT FROBENIUS, 4 1350 - MAX ASPECT FROBENIUS, 5 1351 - MIN_ANGLE, 6 1352 - COLLAPSE RATIO, 7 1353 - MAX ANGLE, 8 1354 - CONDITION, 9 1355 - SCALED JACOBIAN, 10 1356 - SHEAR, 11 1357 - RELATIVE SIZE SQUARED, 12 1358 - SHAPE, 13 1359 - SHAPE AND SIZE, 14 1360 - DISTORTION, 15 1361 - MAX EDGE RATIO, 16 1362 - SKEW, 17 1363 - TAPER, 18 1364 - VOLUME, 19 1365 - STRETCH, 20 1366 - DIAGONAL, 21 1367 - DIMENSION, 22 1368 - ODDY, 23 1369 - SHEAR AND SIZE, 24 1370 - JACOBIAN, 25 1371 - WARPAGE, 26 1372 - ASPECT GAMMA, 27 1373 - AREA, 28 1374 - ASPECT BETA, 29 1375 1376 Examples: 1377 - [meshquality.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/meshquality.py) 1378 1379  1380 """ 1381 qf = vtk.vtkMeshQuality() 1382 qf.SetInputData(self.polydata(False)) 1383 qf.SetTriangleQualityMeasure(metric) 1384 qf.SaveCellQualityOn() 1385 qf.Update() 1386 pd = qf.GetOutput() 1387 self._update(pd) 1388 self.pipeline = OperationNode("compute_quality", parents=[self]) 1389 return self 1390 1391 def check_validity(self, tol=0): 1392 """ 1393 Return an array of possible problematic faces following this convention: 1394 - Valid = 0 1395 - WrongNumberOfPoints = 1 1396 - IntersectingEdges = 2 1397 - IntersectingFaces = 4 1398 - NoncontiguousEdges = 8 1399 - Nonconvex = 10 1400 - OrientedIncorrectly = 20 1401 1402 Arguments: 1403 tol : (float) 1404 value is used as an epsilon for floating point 1405 equality checks throughout the cell checking process. 1406 """ 1407 vald = vtk.vtkCellValidator() 1408 if tol: 1409 vald.SetTolerance(tol) 1410 vald.SetInputData(self._data) 1411 vald.Update() 1412 varr = vald.GetOutput().GetCellData().GetArray("ValidityState") 1413 return vtk2numpy(varr) 1414 1415 1416 @deprecated(reason=vedo.colors.red + "Please use compute_curvature()" + vedo.colors.reset) 1417 def addCurvatureScalars(self, method=0): 1418 """Deprecated. Please use `compute_curvature()`""" 1419 return self.compute_curvature(method) 1420 1421 def compute_curvature(self, method=0): 1422 """ 1423 Add scalars to `Mesh` that contains the curvature calculated in three different ways. 1424 1425 Variable `method` can be: 1426 - 0 = gaussian 1427 - 1 = mean curvature 1428 - 2 = max curvature 1429 - 3 = min curvature 1430 1431 Example: 1432 ```python 1433 from vedo import Torus 1434 Torus().compute_curvature().add_scalarbar().show(axes=1).close() 1435 ``` 1436  1437 """ 1438 curve = vtk.vtkCurvatures() 1439 curve.SetInputData(self._data) 1440 curve.SetCurvatureType(method) 1441 curve.Update() 1442 self._update(curve.GetOutput()) 1443 self._mapper.ScalarVisibilityOn() 1444 return self 1445 1446 def compute_connectivity(self): 1447 """ 1448 Flag a mesh by connectivity: each disconnected region will receive a different Id. 1449 You can access the array of ids through `mesh.pointdata["RegionId"]`. 1450 """ 1451 cf = vtk.vtkConnectivityFilter() 1452 cf.SetInputData(self.polydata(False)) 1453 cf.SetExtractionModeToAllRegions() 1454 cf.ColorRegionsOn() 1455 cf.Update() 1456 return self._update(cf.GetOutput()) 1457 1458 def compute_elevation(self, low=(0, 0, 0), high=(0, 0, 1), vrange=(0, 1)): 1459 """ 1460 Add to `Mesh` a scalar array that contains distance along a specified direction. 1461 1462 Arguments: 1463 low : (list) 1464 one end of the line (small scalar values) 1465 high : (list) 1466 other end of the line (large scalar values) 1467 vrange : (list) 1468 set the range of the scalar 1469 1470 Example: 1471 ```python 1472 from vedo import Sphere 1473 s = Sphere().compute_elevation(low=(0,0,0), high=(1,1,1)) 1474 s.add_scalarbar().show(axes=1).close() 1475 ``` 1476  1477 """ 1478 ef = vtk.vtkElevationFilter() 1479 ef.SetInputData(self.polydata()) 1480 ef.SetLowPoint(low) 1481 ef.SetHighPoint(high) 1482 ef.SetScalarRange(vrange) 1483 ef.Update() 1484 self._update(ef.GetOutput()) 1485 self._mapper.ScalarVisibilityOn() 1486 return self 1487 1488 @deprecated(reason=vedo.colors.red + "Please use add_shadow()" + vedo.colors.reset) 1489 def addShadow(self, *a, **b): 1490 """Please use `add_shadow()`""" 1491 return self.add_shadow(*a, **b) 1492 1493 def add_shadow( 1494 self, 1495 plane, 1496 point, 1497 direction=None, 1498 c=(0.6, 0.6, 0.6), 1499 alpha=1, 1500 culling=0, 1501 ): 1502 """ 1503 Generate a shadow out of an `Mesh` on one of the three Cartesian planes. 1504 The output is a new `Mesh` representing the shadow. 1505 This new mesh is accessible through `mesh.shadow`. 1506 By default the shadow mesh is placed on the bottom wall of the bounding box. 1507 1508 See also `pointcloud.project_on_plane()`. 1509 1510 Arguments: 1511 plane : (str, Plane) 1512 if plane is `str`, plane can be one of `['x', 'y', 'z']`, 1513 represents x-plane, y-plane and z-plane, respectively. 1514 Otherwise, plane should be an instance of `vedo.shapes.Plane` 1515 point : (float, array) 1516 if plane is `str`, point should be a float represents the intercept. 1517 Otherwise, point is the camera point of perspective projection 1518 direction : (list) 1519 direction of oblique projection 1520 culling : (int) 1521 choose between front [1] or backface [-1] culling or None. 1522 1523 Examples: 1524 - [shadow1.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/shadow1.py) 1525 - [airplane1.py](https://github.com/marcomusy/vedo/tree/master/examples/simulations/airplane1.py) 1526 - [airplane2.py](https://github.com/marcomusy/vedo/tree/master/examples/simulations/airplane2.py) 1527 1528  1529 """ 1530 shad = self.clone() 1531 shad.name = "Shadow" 1532 pts = shad.points() 1533 if "x" == plane: 1534 # shad = shad.project_on_plane('x') 1535 # instead do it manually so in case of alpha<1 we dont see glitches due to coplanar points 1536 # we leave a small tolerance of 0.1% in thickness 1537 x0, x1 = self.xbounds() 1538 pts[:, 0] = (pts[:, 0] - (x0 + x1) / 2) / 1000 + self.GetOrigin()[0] 1539 shad.points(pts) 1540 shad.x(point) 1541 elif "y" == plane: 1542 x0, x1 = self.ybounds() 1543 pts[:, 1] = (pts[:, 1] - (x0 + x1) / 2) / 1000 + self.GetOrigin()[1] 1544 shad.points(pts) 1545 shad.y(point) 1546 elif "z" == plane: 1547 x0, x1 = self.zbounds() 1548 pts[:, 2] = (pts[:, 2] - (x0 + x1) / 2) / 1000 + self.GetOrigin()[2] 1549 shad.points(pts)