vedo.core

   1#!/usr/bin/env python3
   2# -*- coding: utf-8 -*-
   3import numpy as np
   4from typing import List, Union, Any
   5from typing_extensions import Self
   6
   7import vedo.vtkclasses as vtki
   8
   9import vedo
  10from vedo import colors
  11from vedo import utils
  12from vedo.transformations import LinearTransform, NonLinearTransform
  13
  14
  15__docformat__ = "google"
  16
  17__doc__ = """
  18Base classes providing functionality to different vedo objects.
  19
  20![](https://vedo.embl.es/images/feats/algorithms_illustration.png)
  21"""
  22
  23__all__ = [
  24    "DataArrayHelper",
  25    "CommonAlgorithms",
  26    "PointAlgorithms",
  27    "VolumeAlgorithms",
  28]
  29
  30warnings = dict(
  31    points_getter=(
  32        "WARNING: points() is deprecated, use vertices instead. Change:\n"
  33        "         mesh.points() -> mesh.vertices\n"
  34        "         (silence this with vedo.core.warnings['points_getter']=False)"
  35    ),
  36    points_setter=(
  37        "WARNING: points() is deprecated, use vertices instead. Change:\n"
  38        "         mesh.points([[x,y,z], ...]) -> mesh.vertices = [[x,y,z], ...]\n"
  39        "         (silence this with vedo.core.warnings['points_getter']=False)"
  40    ),
  41)
  42
  43###############################################################################
  44class DataArrayHelper:
  45    # Internal use only.
  46    # Helper class to manage data associated to either
  47    # points (or vertices) and cells (or faces).
  48    def __init__(self, obj, association):
  49
  50        self.obj = obj
  51        self.association = association
  52
  53    def __getitem__(self, key):
  54
  55        if self.association == 0:
  56            data = self.obj.dataset.GetPointData()
  57
  58        elif self.association == 1:
  59            data = self.obj.dataset.GetCellData()
  60
  61        elif self.association == 2:
  62            data = self.obj.dataset.GetFieldData()
  63
  64            varr = data.GetAbstractArray(key)
  65            if isinstance(varr, vtki.vtkStringArray):
  66                if isinstance(key, int):
  67                    key = data.GetArrayName(key)
  68                n = varr.GetNumberOfValues()
  69                narr = [varr.GetValue(i) for i in range(n)]
  70                return narr
  71                ###########
  72
  73        else:
  74            raise RuntimeError()
  75
  76        if isinstance(key, int):
  77            key = data.GetArrayName(key)
  78
  79        arr = data.GetArray(key)
  80        if not arr:
  81            return None
  82        return utils.vtk2numpy(arr)
  83
  84    def __setitem__(self, key, input_array):
  85
  86        if self.association == 0:
  87            data = self.obj.dataset.GetPointData()
  88            n = self.obj.dataset.GetNumberOfPoints()
  89            self.obj.mapper.SetScalarModeToUsePointData()
  90
  91        elif self.association == 1:
  92            data = self.obj.dataset.GetCellData()
  93            n = self.obj.dataset.GetNumberOfCells()
  94            self.obj.mapper.SetScalarModeToUseCellData()
  95
  96        elif self.association == 2:
  97            data = self.obj.dataset.GetFieldData()
  98            if not utils.is_sequence(input_array):
  99                input_array = [input_array]
 100
 101            if isinstance(input_array[0], str):
 102                varr = vtki.vtkStringArray()
 103                varr.SetName(key)
 104                varr.SetNumberOfComponents(1)
 105                varr.SetNumberOfTuples(len(input_array))
 106                for i, iarr in enumerate(input_array):
 107                    if isinstance(iarr, np.ndarray):
 108                        iarr = iarr.tolist()  # better format
 109                        # Note: a string k can be converted to numpy with
 110                        # import json; k = np.array(json.loads(k))
 111                    varr.InsertValue(i, str(iarr))
 112            else:
 113                try:
 114                    varr = utils.numpy2vtk(input_array, name=key)
 115                except TypeError as e:
 116                    vedo.logger.error(
 117                        f"cannot create metadata with input object:\n"
 118                        f"{input_array}"
 119                        f"\n\nAllowed content examples are:\n"
 120                        f"- flat list of strings ['a','b', 1, [1,2,3], ...]"
 121                        f" (first item must be a string in this case)\n"
 122                        f"  hint: use k = np.array(json.loads(k)) to convert strings\n"
 123                        f"- numpy arrays of any shape"
 124                    )
 125                    raise e
 126
 127            data.AddArray(varr)
 128            return  ############
 129
 130        else:
 131            raise RuntimeError()
 132
 133        if len(input_array) != n:
 134            vedo.logger.error(
 135                f"Error in point/cell data: length of input {len(input_array)}"
 136                f" !=  {n} nr. of elements"
 137            )
 138            raise RuntimeError()
 139
 140        input_array = np.asarray(input_array)
 141        varr = utils.numpy2vtk(input_array, name=key)
 142        data.AddArray(varr)
 143
 144        if len(input_array.shape) == 1:  # scalars
 145            data.SetActiveScalars(key)
 146            try:  # could be a volume mapper
 147                self.obj.mapper.SetScalarRange(data.GetScalars().GetRange())
 148            except AttributeError:
 149                pass
 150        elif len(input_array.shape) == 2 and input_array.shape[1] == 3:  # vectors
 151            if key.lower() == "normals":
 152                data.SetActiveNormals(key)
 153            else:
 154                data.SetActiveVectors(key)
 155
 156    def keys(self) -> List[str]:
 157        """Return the list of available data array names"""
 158        if self.association == 0:
 159            data = self.obj.dataset.GetPointData()
 160        elif self.association == 1:
 161            data = self.obj.dataset.GetCellData()
 162        elif self.association == 2:
 163            data = self.obj.dataset.GetFieldData()
 164        arrnames = []
 165        for i in range(data.GetNumberOfArrays()):
 166            name = ""
 167            if self.association == 2:
 168                name = data.GetAbstractArray(i).GetName()
 169            else:
 170                iarr = data.GetArray(i)
 171                if iarr:
 172                    name = iarr.GetName()
 173            if name:
 174                arrnames.append(name)
 175        return arrnames
 176
 177    def items(self) -> List:
 178        """Return the list of available data array `(names, values)`."""
 179        if self.association == 0:
 180            data = self.obj.dataset.GetPointData()
 181        elif self.association == 1:
 182            data = self.obj.dataset.GetCellData()
 183        elif self.association == 2:
 184            data = self.obj.dataset.GetFieldData()
 185        arrnames = []
 186        for i in range(data.GetNumberOfArrays()):
 187            if self.association == 2:
 188                name = data.GetAbstractArray(i).GetName()
 189            else:
 190                name = data.GetArray(i).GetName()
 191            if name:
 192                arrnames.append((name, self[name]))
 193        return arrnames
 194
 195    def todict(self) -> dict:
 196        """Return a dictionary of the available data arrays."""
 197        return dict(self.items())
 198
 199    def rename(self, oldname: str, newname: str) -> None:
 200        """Rename an array"""
 201        if self.association == 0:
 202            varr = self.obj.dataset.GetPointData().GetArray(oldname)
 203        elif self.association == 1:
 204            varr = self.obj.dataset.GetCellData().GetArray(oldname)
 205        elif self.association == 2:
 206            varr = self.obj.dataset.GetFieldData().GetAbstractArray(oldname)
 207        if varr:
 208            varr.SetName(newname)
 209        else:
 210            vedo.logger.warning(
 211                f"Cannot rename non existing array {oldname} to {newname}"
 212            )
 213
 214    def remove(self, key: Union[int, str]) -> None:
 215        """Remove a data array by name or number"""
 216        if self.association == 0:
 217            self.obj.dataset.GetPointData().RemoveArray(key)
 218        elif self.association == 1:
 219            self.obj.dataset.GetCellData().RemoveArray(key)
 220        elif self.association == 2:
 221            self.obj.dataset.GetFieldData().RemoveArray(key)
 222
 223    def clear(self) -> None:
 224        """Remove all data associated to this object"""
 225        if self.association == 0:
 226            data = self.obj.dataset.GetPointData()
 227        elif self.association == 1:
 228            data = self.obj.dataset.GetCellData()
 229        elif self.association == 2:
 230            data = self.obj.dataset.GetFieldData()
 231        for i in range(data.GetNumberOfArrays()):
 232            if self.association == 2:
 233                name = data.GetAbstractArray(i).GetName()
 234            else:
 235                name = data.GetArray(i).GetName()
 236            data.RemoveArray(name)
 237
 238    def select(self, key: Union[int, str]) -> Any:
 239        """Select one specific array by its name to make it the `active` one."""
 240        # Default (ColorModeToDefault): unsigned char scalars are treated as colors,
 241        # and NOT mapped through the lookup table, while everything else is.
 242        # ColorModeToDirectScalar extends ColorModeToDefault such that all integer
 243        # types are treated as colors with values in the range 0-255
 244        # and floating types are treated as colors with values in the range 0.0-1.0.
 245        # Setting ColorModeToMapScalars means that all scalar data will be mapped
 246        # through the lookup table.
 247        # (Note that for multi-component scalars, the particular component
 248        # to use for mapping can be specified using the SelectColorArray() method.)
 249        if self.association == 0:
 250            data = self.obj.dataset.GetPointData()
 251            self.obj.mapper.SetScalarModeToUsePointData()
 252        else:
 253            data = self.obj.dataset.GetCellData()
 254            self.obj.mapper.SetScalarModeToUseCellData()
 255
 256        if isinstance(key, int):
 257            key = data.GetArrayName(key)
 258
 259        arr = data.GetArray(key)
 260        if not arr:
 261            return self.obj
 262
 263        nc = arr.GetNumberOfComponents()
 264        # print("GetNumberOfComponents", nc)
 265        if nc == 1:
 266            data.SetActiveScalars(key)
 267        elif nc == 2:
 268            data.SetTCoords(arr)
 269        elif nc in (3, 4):
 270            if "rgb" in key.lower(): # type: ignore
 271                data.SetActiveScalars(key)
 272                try:
 273                    # could be a volume mapper
 274                    self.obj.mapper.SetColorModeToDirectScalars()
 275                    data.SetActiveVectors(None) # need this to fix bug in #1066
 276                    # print("SetColorModeToDirectScalars for", key)
 277                except AttributeError:
 278                    pass
 279            else:
 280                data.SetActiveVectors(key)
 281        elif nc == 9:
 282            data.SetActiveTensors(key)
 283        else:
 284            vedo.logger.error(f"Cannot select array {key} with {nc} components")
 285            return self.obj
 286
 287        try:
 288            # could be a volume mapper
 289            self.obj.mapper.SetArrayName(key)
 290            self.obj.mapper.ScalarVisibilityOn()
 291        except AttributeError:
 292            pass
 293
 294        return self.obj
 295
 296    def select_texture_coords(self, key: Union[int,str]) -> Any:
 297        """Select one specific array to be used as texture coordinates."""
 298        if self.association == 0:
 299            data = self.obj.dataset.GetPointData()
 300        else:
 301            vedo.logger.warning("texture coordinates are only available for point data")
 302            return
 303
 304        if isinstance(key, int):
 305            key = data.GetArrayName(key)
 306        data.SetTCoords(data.GetArray(key))
 307        return self.obj
 308
 309    def select_normals(self, key: Union[int,str]) -> Any:
 310        """Select one specific normal array by its name to make it the "active" one."""
 311        if self.association == 0:
 312            data = self.obj.dataset.GetPointData()
 313            self.obj.mapper.SetScalarModeToUsePointData()
 314        else:
 315            data = self.obj.dataset.GetCellData()
 316            self.obj.mapper.SetScalarModeToUseCellData()
 317
 318        if isinstance(key, int):
 319            key = data.GetArrayName(key)
 320        data.SetActiveNormals(key)
 321        return self.obj
 322
 323    def print(self, **kwargs) -> None:
 324        """Print the array names available to terminal"""
 325        colors.printc(self.keys(), **kwargs)
 326
 327    def __repr__(self) -> str:
 328        """Representation"""
 329
 330        def _get_str(pd, header):
 331            out = f"\x1b[2m\x1b[1m\x1b[7m{header}"
 332            if pd.GetNumberOfArrays():
 333                if self.obj.name:
 334                    out += f" in {self.obj.name}"
 335                out += f" contains {pd.GetNumberOfArrays()} array(s)\x1b[0m"
 336                for i in range(pd.GetNumberOfArrays()):
 337                    varr = pd.GetArray(i)
 338                    out += f"\n\x1b[1m\x1b[4mArray name    : {varr.GetName()}\x1b[0m"
 339                    out += "\nindex".ljust(15) + f": {i}"
 340                    t = varr.GetDataType()
 341                    if t in vtki.array_types:
 342                        out += "\ntype".ljust(15)
 343                        out += f": {vtki.array_types[t]}"
 344                    shape = (varr.GetNumberOfTuples(), varr.GetNumberOfComponents())
 345                    out += "\nshape".ljust(15) + f": {shape}"
 346                    out += "\nrange".ljust(15) + f": {np.array(varr.GetRange())}"
 347                    out += "\nmax id".ljust(15) + f": {varr.GetMaxId()}"
 348                    out += "\nlook up table".ljust(15) + f": {bool(varr.GetLookupTable())}"
 349                    out += "\nin-memory size".ljust(15) + f": {varr.GetActualMemorySize()} KB"
 350            else:
 351                out += " is empty.\x1b[0m"
 352            return out
 353
 354        if self.association == 0:
 355            out = _get_str(self.obj.dataset.GetPointData(), "Point Data")
 356        elif self.association == 1:
 357            out = _get_str(self.obj.dataset.GetCellData(), "Cell Data")
 358        elif self.association == 2:
 359            pd = self.obj.dataset.GetFieldData()
 360            if pd.GetNumberOfArrays():
 361                out = "\x1b[2m\x1b[1m\x1b[7mMeta Data"
 362                if self.obj.name:
 363                    out += f" in {self.obj.name}"
 364                out += f" contains {pd.GetNumberOfArrays()} entries\x1b[0m"
 365                for i in range(pd.GetNumberOfArrays()):
 366                    varr = pd.GetAbstractArray(i)
 367                    out += f"\n\x1b[1m\x1b[4mEntry name    : {varr.GetName()}\x1b[0m"
 368                    out += "\nindex".ljust(15) + f": {i}"
 369                    shape = (varr.GetNumberOfTuples(), varr.GetNumberOfComponents())
 370                    out += "\nshape".ljust(15) + f": {shape}"
 371
 372        return out
 373
 374
 375###############################################################################
 376class CommonAlgorithms:
 377    """Common algorithms."""
 378
 379    @property
 380    def pointdata(self):
 381        """
 382        Create and/or return a `numpy.array` associated to points (vertices).
 383        A data array can be indexed either as a string or by an integer number.
 384        E.g.:  `myobj.pointdata["arrayname"]`
 385
 386        Usage:
 387
 388            `myobj.pointdata.keys()` to return the available data array names
 389
 390            `myobj.pointdata.select(name)` to make this array the active one
 391
 392            `myobj.pointdata.remove(name)` to remove this array
 393        """
 394        return DataArrayHelper(self, 0)
 395
 396    @property
 397    def celldata(self):
 398        """
 399        Create and/or return a `numpy.array` associated to cells (faces).
 400        A data array can be indexed either as a string or by an integer number.
 401        E.g.:  `myobj.celldata["arrayname"]`
 402
 403        Usage:
 404
 405            `myobj.celldata.keys()` to return the available data array names
 406
 407            `myobj.celldata.select(name)` to make this array the active one
 408
 409            `myobj.celldata.remove(name)` to remove this array
 410        """
 411        return DataArrayHelper(self, 1)
 412
 413    @property
 414    def metadata(self):
 415        """
 416        Create and/or return a `numpy.array` associated to neither cells nor faces.
 417        A data array can be indexed either as a string or by an integer number.
 418        E.g.:  `myobj.metadata["arrayname"]`
 419
 420        Usage:
 421
 422            `myobj.metadata.keys()` to return the available data array names
 423
 424            `myobj.metadata.select(name)` to make this array the active one
 425
 426            `myobj.metadata.remove(name)` to remove this array
 427        """
 428        return DataArrayHelper(self, 2)
 429
 430    def memory_address(self) -> int:
 431        """
 432        Return a unique memory address integer which may serve as the ID of the
 433        object, or passed to c++ code.
 434        """
 435        # https://www.linkedin.com/pulse/speedup-your-code-accessing-python-vtk-objects-from-c-pletzer/
 436        # https://github.com/tfmoraes/polydata_connectivity
 437        return int(self.dataset.GetAddressAsString("")[5:], 16)
 438
 439    def memory_size(self) -> int:
 440        """Return the size in bytes of the object in memory."""
 441        return self.dataset.GetActualMemorySize()
 442
 443    def modified(self) -> Self:
 444        """Use in conjunction with `tonumpy()` to update any modifications to the image array."""
 445        self.dataset.GetPointData().Modified()
 446        scals = self.dataset.GetPointData().GetScalars()
 447        if scals:
 448            scals.Modified()
 449        return self
 450
 451    def box(self, scale=1, padding=0) -> "vedo.Mesh":
 452        """
 453        Return the bounding box as a new `Mesh` object.
 454
 455        Arguments:
 456            scale : (float)
 457                box size can be scaled by a factor
 458            padding : (float, list)
 459                a constant padding can be added (can be a list `[padx,pady,padz]`)
 460        """
 461        b = self.bounds()
 462        if not utils.is_sequence(padding):
 463            padding = [padding, padding, padding]
 464        length, width, height = b[1] - b[0], b[3] - b[2], b[5] - b[4]
 465        tol = (length + width + height) / 30000  # useful for boxing text
 466        pos = [(b[0] + b[1]) / 2, (b[3] + b[2]) / 2, (b[5] + b[4]) / 2 - tol]
 467        bx = vedo.shapes.Box(
 468            pos,
 469            length * scale + padding[0],
 470            width  * scale + padding[1],
 471            height * scale + padding[2],
 472            c="gray",
 473        )
 474        try:
 475            pr = vtki.vtkProperty()
 476            pr.DeepCopy(self.properties)
 477            bx.actor.SetProperty(pr)
 478            bx.properties = pr
 479        except (AttributeError, TypeError):
 480            pass
 481        bx.flat().lighting("off").wireframe(True)
 482        return bx
 483    
 484    def update_dataset(self, dataset, **kwargs) -> Self:
 485        """Update the dataset of the object with the provided VTK dataset."""
 486        self._update(dataset, **kwargs)
 487        return self
 488
 489    def bounds(self) -> np.ndarray:
 490        """
 491        Get the object bounds.
 492        Returns a list in format `[xmin,xmax, ymin,ymax, zmin,zmax]`.
 493        """
 494        try:  # this is very slow for large meshes
 495            pts = self.vertices
 496            xmin, ymin, zmin = np.min(pts, axis=0)
 497            xmax, ymax, zmax = np.max(pts, axis=0)
 498            return np.array([xmin, xmax, ymin, ymax, zmin, zmax])
 499        except (AttributeError, ValueError):
 500            return np.array(self.dataset.GetBounds())
 501
 502    def xbounds(self, i=None) -> np.ndarray:
 503        """Get the bounds `[xmin,xmax]`. Can specify upper or lower with i (0,1)."""
 504        b = self.bounds()
 505        if i is not None:
 506            return b[i]
 507        return np.array([b[0], b[1]])
 508
 509    def ybounds(self, i=None) -> np.ndarray:
 510        """Get the bounds `[ymin,ymax]`. Can specify upper or lower with i (0,1)."""
 511        b = self.bounds()
 512        if i == 0:
 513            return b[2]
 514        if i == 1:
 515            return b[3]
 516        return np.array([b[2], b[3]])
 517
 518    def zbounds(self, i=None) -> np.ndarray:
 519        """Get the bounds `[zmin,zmax]`. Can specify upper or lower with i (0,1)."""
 520        b = self.bounds()
 521        if i == 0:
 522            return b[4]
 523        if i == 1:
 524            return b[5]
 525        return np.array([b[4], b[5]])
 526
 527    def diagonal_size(self) -> float:
 528        """Get the length of the diagonal of the bounding box."""
 529        b = self.bounds()
 530        return np.sqrt((b[1] - b[0])**2 + (b[3] - b[2])**2 + (b[5] - b[4])**2)
 531
 532    def average_size(self) -> float:
 533        """
 534        Calculate and return the average size of the object.
 535        This is the mean of the vertex distances from the center of mass.
 536        """
 537        coords = self.vertices
 538        cm = np.mean(coords, axis=0)
 539        if coords.shape[0] == 0:
 540            return 0.0
 541        cc = coords - cm
 542        return np.mean(np.linalg.norm(cc, axis=1))
 543
 544    def center_of_mass(self) -> np.ndarray:
 545        """Get the center of mass of the object."""
 546        if isinstance(self, (vedo.RectilinearGrid, vedo.Volume)):
 547            return np.array(self.dataset.GetCenter())
 548        cmf = vtki.new("CenterOfMass")
 549        cmf.SetInputData(self.dataset)
 550        cmf.Update()
 551        c = cmf.GetCenter()
 552        return np.array(c)
 553
 554    def copy_data_from(self, obj: Any) -> Self:
 555        """Copy all data (point and cell data) from this input object"""
 556        self.dataset.GetPointData().PassData(obj.dataset.GetPointData())
 557        self.dataset.GetCellData().PassData(obj.dataset.GetCellData())
 558        self.pipeline = utils.OperationNode(
 559            "copy_data_from",
 560            parents=[self, obj],
 561            comment=f"{obj.__class__.__name__}",
 562            shape="note",
 563            c="#ccc5b9",
 564        )
 565        return self
 566
 567    def inputdata(self):
 568        """Obsolete, use `.dataset` instead."""
 569        colors.printc("WARNING: 'inputdata()' is obsolete, use '.dataset' instead.", c="y")
 570        return self.dataset
 571
 572    @property
 573    def npoints(self):
 574        """Retrieve the number of points (or vertices)."""
 575        return self.dataset.GetNumberOfPoints()
 576
 577    @property
 578    def nvertices(self):
 579        """Retrieve the number of vertices (or points)."""
 580        return self.dataset.GetNumberOfPoints()
 581
 582    @property
 583    def ncells(self):
 584        """Retrieve the number of cells."""
 585        return self.dataset.GetNumberOfCells()
 586
 587    def points(self, pts=None):
 588        """Obsolete, use `self.vertices` or `self.coordinates` instead."""
 589        if pts is None:  ### getter
 590
 591            if warnings["points_getter"]:
 592                colors.printc(warnings["points_getter"], c="y")
 593                warnings["points_getter"] = ""
 594            return self.vertices
 595
 596        else:  ### setter
 597
 598            if warnings["points_setter"]:
 599                colors.printc(warnings["points_setter"], c="y")
 600                warnings["points_setter"] = ""
 601
 602            pts = np.asarray(pts, dtype=np.float32)
 603
 604            if pts.ndim == 1:
 605                ### getter by point index ###################
 606                indices = pts.astype(int)
 607                vpts = self.dataset.GetPoints()
 608                arr = utils.vtk2numpy(vpts.GetData())
 609                return arr[indices]  ###########
 610
 611            ### setter ####################################
 612            if pts.shape[1] == 2:
 613                pts = np.c_[pts, np.zeros(pts.shape[0], dtype=np.float32)]
 614            arr = utils.numpy2vtk(pts, dtype=np.float32)
 615
 616            vpts = self.dataset.GetPoints()
 617            vpts.SetData(arr)
 618            vpts.Modified()
 619            # reset mesh to identity matrix position/rotation:
 620            self.point_locator = None
 621            self.cell_locator = None
 622            self.line_locator = None
 623            self.transform = LinearTransform()
 624            return self
 625
 626    @property
 627    def cell_centers(self):
 628        """
 629        Get the coordinates of the cell centers.
 630
 631        Examples:
 632            - [delaunay2d.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/delaunay2d.py)
 633        
 634        See also: `CellCenters()`.
 635        """
 636        vcen = vtki.new("CellCenters")
 637        vcen.CopyArraysOff()
 638        vcen.SetInputData(self.dataset)
 639        vcen.Update()
 640        return utils.vtk2numpy(vcen.GetOutput().GetPoints().GetData())
 641
 642    @property
 643    def lines(self):
 644        """
 645        Get lines connectivity ids as a python array
 646        formatted as `[[id0,id1], [id3,id4], ...]`
 647
 648        See also: `lines_as_flat_array()`.
 649        """
 650        # Get cell connettivity ids as a 1D array. The vtk format is:
 651        #    [nids1, id0 ... idn, niids2, id0 ... idm,  etc].
 652        try:
 653            arr1d = utils.vtk2numpy(self.dataset.GetLines().GetData())
 654        except AttributeError:
 655            return np.array([])
 656        i = 0
 657        conn = []
 658        n = len(arr1d)
 659        for _ in range(n):
 660            cell = [arr1d[i + k + 1] for k in range(arr1d[i])]
 661            conn.append(cell)
 662            i += arr1d[i] + 1
 663            if i >= n:
 664                break
 665
 666        return conn  # cannot always make a numpy array of it!
 667
 668    @property
 669    def lines_as_flat_array(self):
 670        """
 671        Get lines connectivity ids as a 1D numpy array.
 672        Format is e.g. [2,  10,20,  3, 10,11,12,  2, 70,80, ...]
 673
 674        See also: `lines()`.
 675        """
 676        try:
 677            return utils.vtk2numpy(self.dataset.GetLines().GetData())
 678        except AttributeError:
 679            return np.array([])
 680
 681    def mark_boundaries(self) -> Self:
 682        """
 683        Mark cells and vertices if they lie on a boundary.
 684        A new array called `BoundaryCells` is added to the object.
 685        """
 686        mb = vtki.new("MarkBoundaryFilter")
 687        mb.SetInputData(self.dataset)
 688        mb.Update()
 689        self.dataset.DeepCopy(mb.GetOutput())
 690        self.pipeline = utils.OperationNode("mark_boundaries", parents=[self])
 691        return self
 692
 693    def find_cells_in_bounds(self, xbounds=(), ybounds=(), zbounds=()) -> np.ndarray:
 694        """
 695        Find cells that are within the specified bounds.
 696        """
 697        try:
 698            xbounds = list(xbounds.bounds())
 699        except AttributeError:
 700            pass
 701
 702        if len(xbounds) == 6:
 703            bnds = xbounds
 704        else:
 705            bnds = list(self.bounds())
 706            if len(xbounds) == 2:
 707                bnds[0] = xbounds[0]
 708                bnds[1] = xbounds[1]
 709            if len(ybounds) == 2:
 710                bnds[2] = ybounds[0]
 711                bnds[3] = ybounds[1]
 712            if len(zbounds) == 2:
 713                bnds[4] = zbounds[0]
 714                bnds[5] = zbounds[1]
 715
 716        cell_ids = vtki.vtkIdList()
 717        if not self.cell_locator:
 718            self.cell_locator = vtki.new("CellTreeLocator")
 719            self.cell_locator.SetDataSet(self.dataset)
 720            self.cell_locator.BuildLocator()
 721        self.cell_locator.FindCellsWithinBounds(bnds, cell_ids)
 722        cids = []
 723        for i in range(cell_ids.GetNumberOfIds()):
 724            cid = cell_ids.GetId(i)
 725            cids.append(cid)
 726        return np.array(cids)
 727
 728    def find_cells_along_line(self, p0, p1, tol=0.001) -> np.ndarray:
 729        """
 730        Find cells that are intersected by a line segment.
 731        """
 732        cell_ids = vtki.vtkIdList()
 733        if not self.cell_locator:
 734            self.cell_locator = vtki.new("CellTreeLocator")
 735            self.cell_locator.SetDataSet(self.dataset)
 736            self.cell_locator.BuildLocator()
 737        self.cell_locator.FindCellsAlongLine(p0, p1, tol, cell_ids)
 738        cids = []
 739        for i in range(cell_ids.GetNumberOfIds()):
 740            cid = cell_ids.GetId(i)
 741            cids.append(cid)
 742        return np.array(cids)
 743
 744    def find_cells_along_plane(self, origin, normal, tol=0.001) -> np.ndarray:
 745        """
 746        Find cells that are intersected by a plane.
 747        """
 748        cell_ids = vtki.vtkIdList()
 749        if not self.cell_locator:
 750            self.cell_locator = vtki.new("CellTreeLocator")
 751            self.cell_locator.SetDataSet(self.dataset)
 752            self.cell_locator.BuildLocator()
 753        self.cell_locator.FindCellsAlongPlane(origin, normal, tol, cell_ids)
 754        cids = []
 755        for i in range(cell_ids.GetNumberOfIds()):
 756            cid = cell_ids.GetId(i)
 757            cids.append(cid)
 758        return np.array(cids)
 759
 760    def keep_cell_types(self, types=()):
 761        """
 762        Extract cells of a specific type.
 763
 764        Check the VTK cell types here:
 765        https://vtk.org/doc/nightly/html/vtkCellType_8h.html
 766        """
 767        fe = vtki.new("ExtractCellsByType")
 768        fe.SetInputData(self.dataset)
 769        for t in types:
 770            try:
 771                if utils.is_integer(t):
 772                    it = t
 773                else:
 774                    it = vtki.cell_types[t.upper()]
 775            except KeyError:
 776                vedo.logger.error(f"Cell type '{t}' not recognized")
 777                continue
 778            fe.AddCellType(it)
 779        fe.Update()
 780        self._update(fe.GetOutput())
 781        return self
 782
 783    def map_cells_to_points(self, arrays=(), move=False) -> Self:
 784        """
 785        Interpolate cell data (i.e., data specified per cell or face)
 786        into point data (i.e., data specified at each vertex).
 787        The method of transformation is based on averaging the data values
 788        of all cells using a particular point.
 789
 790        A custom list of arrays to be mapped can be passed in input.
 791
 792        Set `move=True` to delete the original `celldata` array.
 793        """
 794        c2p = vtki.new("CellDataToPointData")
 795        c2p.SetInputData(self.dataset)
 796        if not move:
 797            c2p.PassCellDataOn()
 798        if arrays:
 799            c2p.ClearCellDataArrays()
 800            c2p.ProcessAllArraysOff()
 801            for arr in arrays:
 802                c2p.AddCellDataArray(arr)
 803        else:
 804            c2p.ProcessAllArraysOn()
 805        c2p.Update()
 806        self._update(c2p.GetOutput(), reset_locators=False)
 807        self.mapper.SetScalarModeToUsePointData()
 808        self.pipeline = utils.OperationNode("map_cells_to_points", parents=[self])
 809        return self
 810
 811    @property
 812    def vertices(self):
 813        """Return the vertices (points) coordinates."""
 814        try:
 815            # for polydata and unstructured grid
 816            varr = self.dataset.GetPoints().GetData()
 817        except (AttributeError, TypeError):
 818            try:
 819                # for RectilinearGrid, StructuredGrid
 820                vpts = vtki.vtkPoints()
 821                self.dataset.GetPoints(vpts)
 822                varr = vpts.GetData()
 823            except (AttributeError, TypeError):
 824                try:
 825                    # for ImageData
 826                    v2p = vtki.new("ImageToPoints")
 827                    v2p.SetInputData(self.dataset)
 828                    v2p.Update()
 829                    varr = v2p.GetOutput().GetPoints().GetData()
 830                except AttributeError:
 831                    return np.array([])
 832
 833        return utils.vtk2numpy(varr)
 834
 835    # setter
 836    @vertices.setter
 837    def vertices(self, pts):
 838        """Set vertices (points) coordinates."""
 839        pts = utils.make3d(pts)
 840        arr = utils.numpy2vtk(pts, dtype=np.float32)
 841        try:
 842            vpts = self.dataset.GetPoints()
 843            vpts.SetData(arr)
 844            vpts.Modified()
 845        except (AttributeError, TypeError):
 846            vedo.logger.error(f"Cannot set vertices for {type(self)}")
 847            return self
 848        # reset mesh to identity matrix position/rotation:
 849        self.point_locator = None
 850        self.cell_locator = None
 851        self.line_locator = None
 852        self.transform = LinearTransform()
 853
 854    @property
 855    def coordinates(self):
 856        """Return the vertices (points) coordinates. Same as `vertices`."""
 857        return self.vertices
 858
 859    @coordinates.setter
 860    def coordinates(self, pts):
 861        """Set vertices (points) coordinates. Same as `vertices`."""
 862        self.vertices = pts
 863
 864    @property
 865    def cells_as_flat_array(self):
 866        """
 867        Get cell connectivity ids as a 1D numpy array.
 868        Format is e.g. [3,  10,20,30  4, 10,11,12,13  ...]
 869        """
 870        try:
 871            # valid for unstructured grid
 872            arr1d = utils.vtk2numpy(self.dataset.GetCells().GetData())
 873        except AttributeError:
 874            # valid for polydata
 875            arr1d = utils.vtk2numpy(self.dataset.GetPolys().GetData())
 876        return arr1d
 877
 878    @property
 879    def cells(self):
 880        """
 881        Get the cells connectivity ids as a numpy array.
 882
 883        The output format is: `[[id0 ... idn], [id0 ... idm],  etc]`.
 884        """
 885        try:
 886            # valid for unstructured grid
 887            arr1d = utils.vtk2numpy(self.dataset.GetCells().GetData())
 888        except AttributeError:
 889            try:
 890                # valid for polydata
 891                arr1d = utils.vtk2numpy(self.dataset.GetPolys().GetData())
 892            except AttributeError:
 893                vedo.logger.warning(f"Cannot get cells for {type(self)}")
 894                return np.array([])
 895
 896        # Get cell connettivity ids as a 1D array. vtk format is:
 897        # [nids1, id0 ... idn, niids2, id0 ... idm,  etc].
 898        i = 0
 899        conn = []
 900        n = len(arr1d)
 901        if n:
 902            while True:
 903                cell = [arr1d[i + k] for k in range(1, arr1d[i] + 1)]
 904                conn.append(cell)
 905                i += arr1d[i] + 1
 906                if i >= n:
 907                    break
 908        return conn
 909    
 910    def cell_edge_neighbors(self):
 911        """
 912        Get the cell neighbor indices of each cell.
 913
 914        Returns a python list of lists.
 915        """
 916
 917        def face_to_edges(face):
 918            edges = []
 919            size = len(face)
 920            for i in range(1, size + 1):
 921                if i == size:
 922                    edges.append([face[i - 1], face[0]])
 923                else:
 924                    edges.append([face[i - 1], face[i]])
 925            return edges
 926
 927        pd = self.dataset
 928        pd.BuildLinks()
 929
 930        neicells = []
 931        for i, cell in enumerate(self.cells):
 932            nn = []
 933            for edge in face_to_edges(cell):
 934                neighbors = vtki.vtkIdList()
 935                pd.GetCellEdgeNeighbors(i, edge[0], edge[1], neighbors)
 936                if neighbors.GetNumberOfIds() > 0:
 937                    neighbor = neighbors.GetId(0)
 938                    nn.append(neighbor)
 939            neicells.append(nn)
 940
 941        return neicells
 942
 943
 944    def map_points_to_cells(self, arrays=(), move=False) -> Self:
 945        """
 946        Interpolate point data (i.e., data specified per point or vertex)
 947        into cell data (i.e., data specified per cell).
 948        The method of transformation is based on averaging the data values
 949        of all points defining a particular cell.
 950
 951        A custom list of arrays to be mapped can be passed in input.
 952
 953        Set `move=True` to delete the original `pointdata` array.
 954
 955        Examples:
 956            - [mesh_map2cell.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/mesh_map2cell.py)
 957        """
 958        p2c = vtki.new("PointDataToCellData")
 959        p2c.SetInputData(self.dataset)
 960        if not move:
 961            p2c.PassPointDataOn()
 962        if arrays:
 963            p2c.ClearPointDataArrays()
 964            p2c.ProcessAllArraysOff()
 965            for arr in arrays:
 966                p2c.AddPointDataArray(arr)
 967        else:
 968            p2c.ProcessAllArraysOn()
 969        p2c.Update()
 970        self._update(p2c.GetOutput(), reset_locators=False)
 971        self.mapper.SetScalarModeToUseCellData()
 972        self.pipeline = utils.OperationNode("map_points_to_cells", parents=[self])
 973        return self
 974
 975    def resample_data_from(self, source, tol=None, categorical=False) -> Self:
 976        """
 977        Resample point and cell data from another dataset.
 978        The output has the same structure but its point data have
 979        the resampled values from target.
 980
 981        Use `tol` to set the tolerance used to compute whether
 982        a point in the source is in a cell of the current object.
 983        Points without resampled values, and their cells, are marked as blank.
 984        If the data is categorical, then the resulting data will be determined
 985        by a nearest neighbor interpolation scheme.
 986
 987        Example:
 988        ```python
 989        from vedo import *
 990        m1 = Mesh(dataurl+'bunny.obj')#.add_gaussian_noise(0.1)
 991        pts = m1.vertices
 992        ces = m1.cell_centers
 993        m1.pointdata["xvalues"] = np.power(pts[:,0], 3)
 994        m1.celldata["yvalues"]  = np.power(ces[:,1], 3)
 995        m2 = Mesh(dataurl+'bunny.obj')
 996        m2.resample_data_from(m1)
 997        # print(m2.pointdata["xvalues"])
 998        show(m1, m2 , N=2, axes=1)
 999        ```
1000        """
1001        rs = vtki.new("ResampleWithDataSet")
1002        rs.SetInputData(self.dataset)
1003        rs.SetSourceData(source.dataset)
1004
1005        rs.SetPassPointArrays(True)
1006        rs.SetPassCellArrays(True)
1007        rs.SetPassFieldArrays(True)
1008        rs.SetCategoricalData(categorical)
1009
1010        rs.SetComputeTolerance(True)
1011        if tol:
1012            rs.SetComputeTolerance(False)
1013            rs.SetTolerance(tol)
1014        rs.Update()
1015        self._update(rs.GetOutput(), reset_locators=False)
1016        self.pipeline = utils.OperationNode(
1017            "resample_data_from",
1018            comment=f"{source.__class__.__name__}",
1019            parents=[self, source],
1020        )
1021        return self
1022
1023    def interpolate_data_from(
1024        self,
1025        source,
1026        radius=None,
1027        n=None,
1028        kernel="shepard",
1029        exclude=("Normals",),
1030        on="points",
1031        null_strategy=1,
1032        null_value=0,
1033    ) -> Self:
1034        """
1035        Interpolate over source to port its data onto the current object using various kernels.
1036
1037        If n (number of closest points to use) is set then radius value is ignored.
1038
1039        Check out also:
1040            `probe()` which in many cases can be faster.
1041
1042        Arguments:
1043            kernel : (str)
1044                available kernels are [shepard, gaussian, linear]
1045            null_strategy : (int)
1046                specify a strategy to use when encountering a "null" point
1047                during the interpolation process. Null points occur when the local neighborhood
1048                (of nearby points to interpolate from) is empty.
1049
1050                - Case 0: an output array is created that marks points
1051                  as being valid (=1) or null (invalid =0), and the null_value is set as well
1052                - Case 1: the output data value(s) are set to the provided null_value
1053                - Case 2: simply use the closest point to perform the interpolation.
1054            null_value : (float)
1055                see above.
1056
1057        Examples:
1058            - [interpolate_scalar1.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/interpolate_scalar1.py)
1059            - [interpolate_scalar3.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/interpolate_scalar3.py)
1060            - [interpolate_scalar4.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/interpolate_scalar4.py)
1061            - [image_probe.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/image_probe.py)
1062
1063                ![](https://vedo.embl.es/images/advanced/interpolateMeshArray.png)
1064        """
1065        if radius is None and not n:
1066            vedo.logger.error("in interpolate_data_from(): please set either radius or n")
1067            raise RuntimeError
1068
1069        if on == "points":
1070            points = source.dataset
1071        elif on == "cells":
1072            c2p = vtki.new("CellDataToPointData")
1073            c2p.SetInputData(source.dataset)
1074            c2p.Update()
1075            points = c2p.GetOutput()
1076        else:
1077            vedo.logger.error("in interpolate_data_from(), on must be on points or cells")
1078            raise RuntimeError()
1079
1080        locator = vtki.new("PointLocator")
1081        locator.SetDataSet(points)
1082        locator.BuildLocator()
1083
1084        if kernel.lower() == "shepard":
1085            kern = vtki.new("ShepardKernel")
1086            kern.SetPowerParameter(2)
1087        elif kernel.lower() == "gaussian":
1088            kern = vtki.new("GaussianKernel")
1089            kern.SetSharpness(2)
1090        elif kernel.lower() == "linear":
1091            kern = vtki.new("LinearKernel")
1092        else:
1093            vedo.logger.error("available kernels are: [shepard, gaussian, linear]")
1094            raise RuntimeError()
1095
1096        if n:
1097            kern.SetNumberOfPoints(n)
1098            kern.SetKernelFootprintToNClosest()
1099        else:
1100            kern.SetRadius(radius)
1101            kern.SetKernelFootprintToRadius()
1102
1103        interpolator = vtki.new("PointInterpolator")
1104        interpolator.SetInputData(self.dataset)
1105        interpolator.SetSourceData(points)
1106        interpolator.SetKernel(kern)
1107        interpolator.SetLocator(locator)
1108        interpolator.PassFieldArraysOn()
1109        interpolator.SetNullPointsStrategy(null_strategy)
1110        interpolator.SetNullValue(null_value)
1111        interpolator.SetValidPointsMaskArrayName("ValidPointMask")
1112        for ex in exclude:
1113            interpolator.AddExcludedArray(ex)
1114
1115        # remove arrays that are already present in the source
1116        # this is because the interpolator will ignore them otherwise
1117        for i in range(self.dataset.GetPointData().GetNumberOfArrays()):
1118            name = self.dataset.GetPointData().GetArrayName(i)
1119            if name not in exclude:
1120                self.dataset.GetPointData().RemoveArray(name)
1121
1122        interpolator.Update()
1123
1124        if on == "cells":
1125            p2c = vtki.new("PointDataToCellData")
1126            p2c.SetInputData(interpolator.GetOutput())
1127            p2c.Update()
1128            cpoly = p2c.GetOutput()
1129        else:
1130            cpoly = interpolator.GetOutput()
1131
1132        self._update(cpoly, reset_locators=False)
1133
1134        self.pipeline = utils.OperationNode("interpolate_data_from", parents=[self, source])
1135        return self
1136
1137    def add_ids(self) -> Self:
1138        """
1139        Generate point and cell ids arrays.
1140
1141        Two new arrays are added to the mesh: `PointID` and `CellID`.
1142        """
1143        ids = vtki.new("IdFilter")
1144        ids.SetInputData(self.dataset)
1145        ids.PointIdsOn()
1146        ids.CellIdsOn()
1147        ids.FieldDataOff()
1148        ids.SetPointIdsArrayName("PointID")
1149        ids.SetCellIdsArrayName("CellID")
1150        ids.Update()
1151        self._update(ids.GetOutput(), reset_locators=False)
1152        self.pipeline = utils.OperationNode("add_ids", parents=[self])
1153        return self
1154
1155    def gradient(self, input_array=None, on="points", fast=False) -> np.ndarray:
1156        """
1157        Compute and return the gradiend of the active scalar field as a numpy array.
1158
1159        Arguments:
1160            input_array : (str)
1161                array of the scalars to compute the gradient,
1162                if None the current active array is selected
1163            on : (str)
1164                compute either on 'points' or 'cells' data
1165            fast : (bool)
1166                if True, will use a less accurate algorithm
1167                that performs fewer derivative calculations (and is therefore faster).
1168
1169        Examples:
1170            - [isolines.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/isolines.py)
1171
1172            ![](https://user-images.githubusercontent.com/32848391/72433087-f00a8780-3798-11ea-9778-991f0abeca70.png)
1173        """
1174        gra = vtki.new("GradientFilter")
1175        if on.startswith("p"):
1176            varr = self.dataset.GetPointData()
1177            tp = vtki.vtkDataObject.FIELD_ASSOCIATION_POINTS
1178        elif on.startswith("c"):
1179            varr = self.dataset.GetCellData()
1180            tp = vtki.vtkDataObject.FIELD_ASSOCIATION_CELLS
1181        else:
1182            vedo.logger.error(f"in gradient: unknown option {on}")
1183            raise RuntimeError
1184
1185        if input_array is None:
1186            if varr.GetScalars():
1187                input_array = varr.GetScalars().GetName()
1188            else:
1189                vedo.logger.error(f"in gradient: no scalars found for {on}")
1190                raise RuntimeError
1191
1192        gra.SetInputData(self.dataset)
1193        gra.SetInputScalars(tp, input_array)
1194        gra.SetResultArrayName("Gradient")
1195        gra.SetFasterApproximation(fast)
1196        gra.ComputeDivergenceOff()
1197        gra.ComputeVorticityOff()
1198        gra.ComputeGradientOn()
1199        gra.Update()
1200        # self._update(gra.GetOutput(), reset_locators=False)
1201        if on.startswith("p"):
1202            gvecs = utils.vtk2numpy(gra.GetOutput().GetPointData().GetArray("Gradient"))
1203        else:
1204            gvecs = utils.vtk2numpy(gra.GetOutput().GetCellData().GetArray("Gradient"))
1205        return gvecs
1206
1207    def divergence(self, array_name=None, on="points", fast=False) -> np.ndarray:
1208        """
1209        Compute and return the divergence of a vector field as a numpy array.
1210
1211        Arguments:
1212            array_name : (str)
1213                name of the array of vectors to compute the divergence,
1214                if None the current active array is selected
1215            on : (str)
1216                compute either on 'points' or 'cells' data
1217            fast : (bool)
1218                if True, will use a less accurate algorithm
1219                that performs fewer derivative calculations (and is therefore faster).
1220        """
1221        div = vtki.new("GradientFilter")
1222        if on.startswith("p"):
1223            varr = self.dataset.GetPointData()
1224            tp = vtki.vtkDataObject.FIELD_ASSOCIATION_POINTS
1225        elif on.startswith("c"):
1226            varr = self.dataset.GetCellData()
1227            tp = vtki.vtkDataObject.FIELD_ASSOCIATION_CELLS
1228        else:
1229            vedo.logger.error(f"in divergence(): unknown option {on}")
1230            raise RuntimeError
1231
1232        if array_name is None:
1233            if varr.GetVectors():
1234                array_name = varr.GetVectors().GetName()
1235            else:
1236                vedo.logger.error(f"in divergence(): no vectors found for {on}")
1237                raise RuntimeError
1238
1239        div.SetInputData(self.dataset)
1240        div.SetInputScalars(tp, array_name)
1241        div.ComputeDivergenceOn()
1242        div.ComputeGradientOff()
1243        div.ComputeVorticityOff()
1244        div.SetDivergenceArrayName("Divergence")
1245        div.SetFasterApproximation(fast)
1246        div.Update()
1247        # self._update(div.GetOutput(), reset_locators=False)
1248        if on.startswith("p"):
1249            dvecs = utils.vtk2numpy(div.GetOutput().GetPointData().GetArray("Divergence"))
1250        else:
1251            dvecs = utils.vtk2numpy(div.GetOutput().GetCellData().GetArray("Divergence"))
1252        return dvecs
1253
1254    def vorticity(self, array_name=None, on="points", fast=False) -> np.ndarray:
1255        """
1256        Compute and return the vorticity of a vector field as a numpy array.
1257
1258        Arguments:
1259            array_name : (str)
1260                name of the array to compute the vorticity,
1261                if None the current active array is selected
1262            on : (str)
1263                compute either on 'points' or 'cells' data
1264            fast : (bool)
1265                if True, will use a less accurate algorithm
1266                that performs fewer derivative calculations (and is therefore faster).
1267        """
1268        vort = vtki.new("GradientFilter")
1269        if on.startswith("p"):
1270            varr = self.dataset.GetPointData()
1271            tp = vtki.vtkDataObject.FIELD_ASSOCIATION_POINTS
1272        elif on.startswith("c"):
1273            varr = self.dataset.GetCellData()
1274            tp = vtki.vtkDataObject.FIELD_ASSOCIATION_CELLS
1275        else:
1276            vedo.logger.error(f"in vorticity(): unknown option {on}")
1277            raise RuntimeError
1278
1279        if array_name is None:
1280            if varr.GetVectors():
1281                array_name = varr.GetVectors().GetName()
1282            else:
1283                vedo.logger.error(f"in vorticity(): no vectors found for {on}")
1284                raise RuntimeError
1285
1286        vort.SetInputData(self.dataset)
1287        vort.SetInputScalars(tp, array_name)
1288        vort.ComputeDivergenceOff()
1289        vort.ComputeGradientOff()
1290        vort.ComputeVorticityOn()
1291        vort.SetVorticityArrayName("Vorticity")
1292        vort.SetFasterApproximation(fast)
1293        vort.Update()
1294        if on.startswith("p"):
1295            vvecs = utils.vtk2numpy(vort.GetOutput().GetPointData().GetArray("Vorticity"))
1296        else:
1297            vvecs = utils.vtk2numpy(vort.GetOutput().GetCellData().GetArray("Vorticity"))
1298        return vvecs
1299
1300    def probe(
1301            self,
1302            source,
1303            categorical=False,
1304            snap=False,
1305            tol=0,
1306        ) -> Self:
1307        """
1308        Takes a data set and probes its scalars at the specified points in space.
1309
1310        Note that a mask is also output with valid/invalid points which can be accessed
1311        with `mesh.pointdata['ValidPointMask']`.
1312
1313        Arguments:
1314            source : any dataset
1315                the data set to probe.
1316            categorical : bool
1317                control whether the source pointdata is to be treated as categorical.
1318            snap : bool
1319                snap to the cell with the closest point if no cell was found
1320            tol : float
1321                the tolerance to use when performing the probe.
1322
1323        Check out also:
1324            `interpolate_data_from()` and `tovolume()`
1325
1326        Examples:
1327            - [probe_points.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/probe_points.py)
1328
1329                ![](https://vedo.embl.es/images/volumetric/probePoints.png)
1330        """
1331        probe_filter = vtki.new("ProbeFilter")
1332        probe_filter.SetSourceData(source.dataset)
1333        probe_filter.SetInputData(self.dataset)
1334        probe_filter.PassCellArraysOn()
1335        probe_filter.PassFieldArraysOn()
1336        probe_filter.PassPointArraysOn()
1337        probe_filter.SetCategoricalData(categorical)
1338        probe_filter.ComputeToleranceOff()
1339        if tol:
1340            probe_filter.ComputeToleranceOn()
1341            probe_filter.SetTolerance(tol)
1342        probe_filter.SetSnapToCellWithClosestPoint(snap)
1343        probe_filter.Update()
1344        self._update(probe_filter.GetOutput(), reset_locators=False)
1345        self.pipeline = utils.OperationNode("probe", parents=[self, source])
1346        self.pointdata.rename("vtkValidPointMask", "ValidPointMask")
1347        return self
1348
1349    def compute_cell_size(self) -> Self:
1350        """
1351        Add to this object a cell data array
1352        containing the area, volume and edge length
1353        of the cells (when applicable to the object type).
1354
1355        Array names are: `Area`, `Volume`, `Length`.
1356        """
1357        csf = vtki.new("CellSizeFilter")
1358        csf.SetInputData(self.dataset)
1359        csf.SetComputeArea(1)
1360        csf.SetComputeVolume(1)
1361        csf.SetComputeLength(1)
1362        csf.SetComputeVertexCount(0)
1363        csf.SetAreaArrayName("Area")
1364        csf.SetVolumeArrayName("Volume")
1365        csf.SetLengthArrayName("Length")
1366        csf.Update()
1367        self._update(csf.GetOutput(), reset_locators=False)
1368        return self
1369
1370    def generate_random_data(self) -> Self:
1371        """Fill a dataset with random attributes"""
1372        gen = vtki.new("RandomAttributeGenerator")
1373        gen.SetInputData(self.dataset)
1374        gen.GenerateAllDataOn()
1375        gen.SetDataTypeToFloat()
1376        gen.GeneratePointNormalsOff()
1377        gen.GeneratePointTensorsOn()
1378        gen.GenerateCellScalarsOn()
1379        gen.Update()
1380        self._update(gen.GetOutput(), reset_locators=False)
1381        self.pipeline = utils.OperationNode("generate_random_data", parents=[self])
1382        return self
1383
1384    def integrate_data(self) -> dict:
1385        """
1386        Integrate point and cell data arrays while computing length,
1387        area or volume of the domain. It works for 1D, 2D or 3D cells.
1388    
1389        For volumetric datasets, this filter ignores all but 3D cells.
1390        It will not compute the volume contained in a closed surface.
1391
1392        Returns a dictionary with keys: `pointdata`, `celldata`, `metadata`,
1393        which contain the integration result for the corresponding attributes.
1394
1395        Examples:
1396            ```python
1397            from vedo import *
1398            surf = Sphere(res=100)
1399            surf.pointdata['scalars'] = np.ones(surf.npoints)
1400            data = surf.integrate_data()
1401            print(data['pointdata']['scalars'], "is equal to 4pi", 4*np.pi)
1402            ```
1403
1404            ```python
1405            from vedo import *
1406
1407            xcoords1 = np.arange(0, 2.2, 0.2)
1408            xcoords2 = sqrt(np.arange(0, 4.2, 0.2))
1409
1410            ycoords = np.arange(0, 1.2, 0.2)
1411
1412            surf1 = Grid(s=(xcoords1, ycoords)).rotate_y(-45).lw(2)
1413            surf2 = Grid(s=(xcoords2, ycoords)).rotate_y(-45).lw(2)
1414
1415            surf1.pointdata['scalars'] = surf1.vertices[:,2]
1416            surf2.pointdata['scalars'] = surf2.vertices[:,2]
1417
1418            data1 = surf1.integrate_data()
1419            data2 = surf2.integrate_data()
1420
1421            print(data1['pointdata']['scalars'],
1422                "is equal to",
1423                data2['pointdata']['scalars'],
1424                "even if the grids are different!",
1425                "(= the volume under the surface)"
1426            )
1427            show(surf1, surf2, N=2, axes=1).close()
1428            ```
1429        """
1430        vinteg = vtki.new("IntegrateAttributes")
1431        vinteg.SetInputData(self.dataset)
1432        vinteg.Update()
1433        ugrid = vedo.UnstructuredGrid(vinteg.GetOutput())
1434        data = dict(
1435            pointdata=ugrid.pointdata.todict(),
1436            celldata=ugrid.celldata.todict(),
1437            metadata=ugrid.metadata.todict(),
1438        )
1439        return data
1440
1441    def write(self, filename, binary=True) -> None:
1442        """Write object to file."""
1443        out = vedo.file_io.write(self, filename, binary)
1444        out.pipeline = utils.OperationNode(
1445            "write", parents=[self], comment=filename[:15], shape="folder", c="#8a817c"
1446        )
1447
1448    def tomesh(self, bounds=(), shrink=0) -> "vedo.Mesh":
1449        """
1450        Extract boundary geometry from dataset (or convert data to polygonal type).
1451
1452        Two new arrays are added to the mesh: `OriginalCellIds` and `OriginalPointIds`
1453        to keep track of the original mesh elements.
1454
1455        Arguments:
1456            bounds : (list)
1457                specify a sub-region to extract
1458            shrink : (float)
1459                shrink the cells to a fraction of their original size
1460        """
1461        geo = vtki.new("GeometryFilter")
1462
1463        if shrink:
1464            sf = vtki.new("ShrinkFilter")
1465            sf.SetInputData(self.dataset)
1466            sf.SetShrinkFactor(shrink)
1467            sf.Update()
1468            geo.SetInputData(sf.GetOutput())
1469        else:
1470            geo.SetInputData(self.dataset)
1471
1472        geo.SetPassThroughCellIds(1)
1473        geo.SetPassThroughPointIds(1)
1474        geo.SetOriginalCellIdsName("OriginalCellIds")
1475        geo.SetOriginalPointIdsName("OriginalPointIds")
1476        geo.SetNonlinearSubdivisionLevel(1)
1477        # geo.MergingOff() # crashes on StructuredGrids
1478        if bounds:
1479            geo.SetExtent(bounds)
1480            geo.ExtentClippingOn()
1481        geo.Update()
1482        msh = vedo.mesh.Mesh(geo.GetOutput())
1483        msh.pipeline = utils.OperationNode("tomesh", parents=[self], c="#9e2a2b")
1484        return msh
1485
1486    def signed_distance(self, dims=(20, 20, 20), bounds=None, invert=False, max_radius=None) -> "vedo.Volume":
1487        """
1488        Compute the `Volume` object whose voxels contains the signed distance from
1489        the object. The calling object must have "Normals" defined.
1490
1491        Arguments:
1492            bounds : (list, actor)
1493                bounding box sizes
1494            dims : (list)
1495                dimensions (nr. of voxels) of the output volume.
1496            invert : (bool)
1497                flip the sign
1498            max_radius : (float)
1499                specify how far out to propagate distance calculation
1500
1501        Examples:
1502            - [distance2mesh.py](https://github.com/marcomusy/vedo/blob/master/examples/basic/distance2mesh.py)
1503
1504                ![](https://vedo.embl.es/images/basic/distance2mesh.png)
1505        """
1506        if bounds is None:
1507            bounds = self.bounds()
1508        if max_radius is None:
1509            max_radius = self.diagonal_size() / 2
1510        dist = vtki.new("SignedDistance")
1511        dist.SetInputData(self.dataset)
1512        dist.SetRadius(max_radius)
1513        dist.SetBounds(bounds)
1514        dist.SetDimensions(dims)
1515        dist.Update()
1516        img = dist.GetOutput()
1517        if invert:
1518            mat = vtki.new("ImageMathematics")
1519            mat.SetInput1Data(img)
1520            mat.SetOperationToMultiplyByK()
1521            mat.SetConstantK(-1)
1522            mat.Update()
1523            img = mat.GetOutput()
1524
1525        vol = vedo.Volume(img)
1526        vol.name = "SignedDistanceVolume"
1527        vol.pipeline = utils.OperationNode(
1528            "signed_distance",
1529            parents=[self],
1530            comment=f"dims={tuple(vol.dimensions())}",
1531            c="#e9c46a:#0096c7",
1532        )
1533        return vol
1534    
1535    def unsigned_distance(
1536            self, dims=(25,25,25), bounds=(), max_radius=0, cap_value=0) -> "vedo.Volume":
1537        """
1538        Compute the `Volume` object whose voxels contains the unsigned distance. 
1539        """
1540        dist = vtki.new("UnsignedDistance")
1541        dist.SetInputData(self.dataset)
1542        dist.SetDimensions(dims)
1543
1544        if len(bounds) == 6:
1545            dist.SetBounds(bounds)
1546        # elif bounds == "auto":
1547        #     dist.AdjustBoundsOn()
1548        else:
1549            dist.SetBounds(self.bounds())
1550        if not max_radius:
1551            max_radius = self.diagonal_size() / 10
1552        dist.SetRadius(max_radius)
1553
1554        if self.point_locator:
1555            dist.SetLocator(self.point_locator)
1556        
1557        if cap_value is not None:
1558            dist.CappingOn()
1559            dist.SetCapValue(cap_value)
1560        dist.SetOutputScalarTypeToFloat()
1561        dist.Update()
1562        vol = vedo.Volume(dist.GetOutput())
1563        vol.name = "UnsignedDistanceVolume"
1564        vol.pipeline = utils.OperationNode(
1565            "unsigned_distance", parents=[self], c="#e9c46a:#0096c7")
1566        return vol
1567
1568    def smooth_data(self, 
1569            niter=10, relaxation_factor=0.1, strategy=0, mask=None,
1570            mode="distance2",
1571            exclude=("Normals", "TextureCoordinates"),
1572        ) -> Self:
1573        """
1574        Smooth point attribute data using distance weighted Laplacian kernel.
1575        The effect is to blur regions of high variation and emphasize low variation regions.
1576
1577        A central concept of this method is the point smoothing stencil.
1578        A smoothing stencil for a point p(i) is the list of points p(j) which connect to p(i) via an edge.
1579        To smooth the attributes of point p(i), p(i)'s attribute data a(i) are iteratively averaged using
1580        the distance weighted average of the attributes of a(j) (the weights w[j] sum to 1).
1581        This averaging process is repeated until the maximum number of iterations is reached.
1582
1583        The relaxation factor (R) is also important as the smoothing process proceeds in an iterative fashion.
1584        The a(i+1) attributes are determined from the a(i) attributes as follows:
1585            a(i+1) = (1-R)*a(i) + R*sum(w(j)*a(j))
1586    
1587        Convergence occurs faster for larger relaxation factors.
1588        Typically a small number of iterations is required for large relaxation factors,
1589        and in cases where only points adjacent to the boundary are being smoothed, a single iteration with R=1 may be
1590        adequate (i.e., just a distance weighted average is computed).
1591
1592        Warning:
1593            Certain data attributes cannot be correctly interpolated. For example, surface normals are expected to be |n|=1;
1594            after attribute smoothing this constraint is likely to be violated.
1595            Other vectors and tensors may suffer from similar issues.
1596            In such a situation, specify `exclude=...` which will not be smoothed (and simply passed through to the output).
1597            Distance weighting function is based on averaging, 1/r, or 1/r**2 weights, where r is the distance
1598            between the point to be smoothed and an edge connected neighbor (defined by the smoothing stencil).
1599            The weights are normalized so that sum(w(i))==1. When smoothing based on averaging, the weights are simply 1/n,
1600            where n is the number of connected points in the stencil.
1601            The smoothing process reduces high frequency information in the data attributes.
1602            With excessive smoothing (large numbers of iterations, and/or a large relaxation factor) important details may be lost,
1603            and the attributes will move towards an "average" value.
1604            While this filter will process any dataset type, if the input data is a 3D image volume, it's likely much faster to use
1605            an image-based algorithm to perform data smoothing.
1606            To determine boundary points in polygonal data, edges used by only one cell are considered boundary
1607            (and hence the associated points defining the edge). 
1608
1609        Arguments:
1610            niter : (int)
1611                number of iterations
1612            relaxation_factor : (float)
1613                relaxation factor controlling the amount of Laplacian smoothing applied
1614            strategy : (int)
1615                strategy to use for Laplacian smoothing
1616                    - 0: use all points, all point data attributes are smoothed
1617                    - 1: smooth all point attribute data except those on the boundary
1618                    - 2: only point data connected to a boundary point are smoothed
1619            mask : (str, np.ndarray)
1620                array to be used as a mask (ignore then the strategy keyword)
1621            mode : (str)
1622                smoothing mode, either "distance2", "distance" or "average"
1623                    - distance**2 weighted (i.e., 1/r**2 interpolation weights)
1624                    - distance weighted (i.e., 1/r) approach;
1625                    - simple average of all connected points in the stencil
1626            exclude : (list)
1627                list of arrays to be excluded from smoothing
1628        """
1629        try:
1630            saf = vtki.new("AttributeSmoothingFilter")
1631        except:
1632            vedo.logger.error("smooth_data() only avaialble in VTK>=9.3.0")
1633            return self
1634        saf.SetInputData(self.dataset)
1635        saf.SetRelaxationFactor(relaxation_factor)
1636        saf.SetNumberOfIterations(niter)
1637
1638        for ex in exclude:
1639            saf.AddExcludedArray(ex)
1640
1641        if mode == "distance":
1642            saf.SetWeightsTypeToDistance()
1643        elif mode == "distance2":
1644            saf.SetWeightsTypeToDistance2()
1645        elif mode == "average":
1646            saf.SetWeightsTypeToAverage()
1647        else:
1648            vedo.logger.error(f"smooth_data(): unknown mode {mode}")
1649            raise TypeError
1650
1651        saf.SetSmoothingStrategy(strategy)
1652        if mask is not None:
1653            saf.SetSmoothingStrategyToSmoothingMask()
1654            if isinstance(mask, str):
1655                mask_ = self.dataset.GetPointData().GetArray(mask)
1656                if not mask_:
1657                    vedo.logger.error(f"smooth_data(): mask array {mask} not found")
1658                    return self
1659                mask_array = vtki.vtkUnsignedCharArray()
1660                mask_array.ShallowCopy(mask_)
1661                mask_array.SetName(mask_.GetName())
1662            else:
1663                mask_array = utils.numpy2vtk(mask, dtype=np.uint8)
1664            saf.SetSmoothingMask(mask_array)
1665
1666        saf.Update()
1667
1668        self._update(saf.GetOutput())
1669        self.pipeline = utils.OperationNode(
1670            "smooth_data", comment=f"strategy {strategy}", parents=[self], c="#9e2a2b"
1671        )
1672        return self
1673        
1674    def compute_streamlines(
1675            self, 
1676            seeds: Any, 
1677            integrator="rk4",
1678            direction="forward",
1679            initial_step_size=None,
1680            max_propagation=None,
1681            max_steps=10000,
1682            step_length=0,
1683            surface_constrained=False,
1684            compute_vorticity=False,
1685        ) -> Union["vedo.Lines", None]:
1686        """
1687        Integrate a vector field to generate streamlines.
1688
1689        Arguments:
1690            seeds : (Mesh, Points, list)
1691                starting points of the streamlines
1692            integrator : (str)
1693                type of integration method to be used:
1694                    - "rk2" (Runge-Kutta 2)
1695                    - "rk4" (Runge-Kutta 4)
1696                    - "rk45" (Runge-Kutta 45)
1697            direction : (str)
1698                direction of integration, either "forward", "backward" or "both"
1699            initial_step_size : (float)
1700                initial step size used for line integration
1701            max_propagation : (float)
1702                maximum length of a streamline expressed in absolute units
1703            max_steps : (int)
1704                maximum number of steps for a streamline
1705            step_length : (float)
1706                maximum length of a step expressed in absolute units
1707            surface_constrained : (bool)
1708                whether to stop integrating when the streamline leaves the surface
1709            compute_vorticity : (bool)
1710                whether to compute the vorticity at each streamline point
1711        """
1712        b = self.dataset.GetBounds()
1713        size = (b[5]-b[4] + b[3]-b[2] + b[1]-b[0]) / 3
1714        if initial_step_size is None:
1715            initial_step_size = size / 1000.0
1716
1717        if max_propagation is None:
1718            max_propagation = size * 2
1719
1720        if utils.is_sequence(seeds):
1721            seeds = vedo.Points(seeds)
1722
1723        sti = vtki.new("StreamTracer")
1724        sti.SetSourceData(seeds.dataset)
1725        if isinstance(self, vedo.RectilinearGrid):
1726            sti.SetInputData(vedo.UnstructuredGrid(self.dataset).dataset)
1727        else:
1728            sti.SetInputDataObject(self.dataset)
1729
1730        sti.SetInitialIntegrationStep(initial_step_size)
1731        sti.SetComputeVorticity(compute_vorticity)
1732        sti.SetMaximumNumberOfSteps(max_steps)
1733        sti.SetMaximumPropagation(max_propagation)
1734        sti.SetSurfaceStreamlines(surface_constrained)
1735        if step_length:
1736            sti.SetMaximumIntegrationStep(step_length)
1737
1738        if "for" in direction:
1739            sti.SetIntegrationDirectionToForward()
1740        elif "back" in direction:
1741            sti.SetIntegrationDirectionToBackward()
1742        elif "both" in direction:
1743            sti.SetIntegrationDirectionToBoth()
1744        else:
1745            vedo.logger.error(f"in compute_streamlines(), unknown direction {direction}")
1746            return None
1747
1748        if integrator == "rk2":
1749            sti.SetIntegratorTypeToRungeKutta2()
1750        elif integrator == "rk4":
1751            sti.SetIntegratorTypeToRungeKutta4()
1752        elif integrator == "rk45":
1753            sti.SetIntegratorTypeToRungeKutta45()
1754        else:
1755            vedo.logger.error(f"in compute_streamlines(), unknown integrator {integrator}")
1756            return None
1757
1758        sti.Update()
1759
1760        stlines = vedo.shapes.Lines(sti.GetOutput(), lw=4)
1761        stlines.name = "StreamLines"
1762        self.pipeline = utils.OperationNode(
1763            "compute_streamlines", comment=f"{integrator}", parents=[self, seeds], c="#9e2a2b"
1764        )
1765        return stlines
1766
1767###############################################################################
1768class PointAlgorithms(CommonAlgorithms):
1769    """Methods for point clouds."""
1770
1771    def apply_transform(self, LT: Any, deep_copy=True) -> Self:
1772        """
1773        Apply a linear or non-linear transformation to the mesh polygonal data.
1774
1775        Example:
1776        ```python
1777        from vedo import Cube, show, settings
1778        settings.use_parallel_projection = True
1779        c1 = Cube().rotate_z(25).pos(2,1).mirror().alpha(0.5)
1780        T = c1.transform  # rotate by 5 degrees, place at (2,1)
1781        c2 = Cube().c('red4').wireframe().lw(10).lighting('off')
1782        c2.apply_transform(T)
1783        show(c1, c2, "The 2 cubes should overlap!", axes=1).close()
1784        ```
1785
1786        ![](https://vedo.embl.es/images/feats/apply_transform.png)
1787        """
1788        if self.dataset.GetNumberOfPoints() == 0:
1789            return self
1790
1791        if isinstance(LT, LinearTransform):
1792            LT_is_linear = True
1793            tr = LT.T
1794            if LT.is_identity():
1795                return self
1796        
1797        elif isinstance(LT, (vtki.vtkMatrix4x4, vtki.vtkLinearTransform)) or utils.is_sequence(LT):
1798            LT_is_linear = True
1799            LT = LinearTransform(LT)
1800            tr = LT.T
1801            if LT.is_identity():
1802                return self
1803
1804        elif isinstance(LT, NonLinearTransform):
1805            LT_is_linear = False
1806            tr = LT.T
1807            self.transform = LT  # reset
1808
1809        elif isinstance(LT, vtki.vtkThinPlateSplineTransform):
1810            LT_is_linear = False
1811            tr = LT
1812            self.transform = NonLinearTransform(LT)  # reset
1813
1814        else:
1815            vedo.logger.error(f"apply_transform(), unknown input type:\n{LT}")
1816            return self
1817
1818        ################
1819        if LT_is_linear:
1820            try:
1821                # self.transform might still not be linear
1822                self.transform.concatenate(LT)
1823            except AttributeError:
1824                # in that case reset it
1825                self.transform = LinearTransform()
1826
1827        ################
1828        if isinstance(self.dataset, vtki.vtkPolyData):
1829            tp = vtki.new("TransformPolyDataFilter")
1830        elif isinstance(self.dataset, vtki.vtkUnstructuredGrid):
1831            tp = vtki.new("TransformFilter")
1832            tp.TransformAllInputVectorsOn()
1833        # elif isinstance(self.dataset, vtki.vtkImageData):
1834        #     tp = vtki.new("ImageReslice")
1835        #     tp.SetInterpolationModeToCubic()
1836        #     tp.SetResliceTransform(tr)
1837        else:
1838            vedo.logger.error(f"apply_transform(), unknown input type: {[self.dataset]}")
1839            return self
1840
1841        tp.SetTransform(tr)
1842        tp.SetInputData(self.dataset)
1843        tp.Update()
1844        out = tp.GetOutput()
1845
1846        if deep_copy:
1847            self.dataset.DeepCopy(out)
1848        else:
1849            self.dataset.ShallowCopy(out)
1850
1851        # reset the locators
1852        self.point_locator = None
1853        self.cell_locator = None
1854        self.line_locator = None
1855        return self
1856
1857    def apply_transform_from_actor(self) -> LinearTransform:
1858        """
1859        Apply the current transformation of the actor to the data.
1860        Useful when manually moving an actor (eg. when pressing "a").
1861        Returns the `LinearTransform` object.
1862
1863        Note that this method is automatically called when the window is closed,
1864        or the interactor style is changed.
1865        """
1866        M = self.actor.GetMatrix()
1867        self.apply_transform(M)
1868        iden = vtki.vtkMatrix4x4()
1869        self.actor.PokeMatrix(iden)
1870        return LinearTransform(M)
1871
1872    def pos(self, x=None, y=None, z=None) -> Self:
1873        """Set/Get object position."""
1874        if x is None:  # get functionality
1875            return self.transform.position
1876
1877        if z is None and y is None:  # assume x is of the form (x,y,z)
1878            if len(x) == 3:
1879                x, y, z = x
1880            else:
1881                x, y = x
1882                z = 0
1883        elif z is None:  # assume x,y is of the form x, y
1884            z = 0
1885
1886        q = self.transform.position
1887        delta = [x, y, z] - q
1888        if delta[0] == delta[1] == delta[2] == 0:
1889            return self
1890        LT = LinearTransform().translate(delta)
1891        return self.apply_transform(LT)
1892
1893    def shift(self, dx=0, dy=0, dz=0) -> Self:
1894        """Add a vector to the current object position."""
1895        if utils.is_sequence(dx):
1896            dx, dy, dz = utils.make3d(dx)
1897        if dx == dy == dz == 0:
1898            return self
1899        LT = LinearTransform().translate([dx, dy, dz])
1900        return self.apply_transform(LT)
1901
1902    def x(self, val=None) -> Self:
1903        """Set/Get object position along x axis."""
1904        p = self.transform.position
1905        if val is None:
1906            return p[0]
1907        self.pos(val, p[1], p[2])
1908        return self
1909
1910    def y(self, val=None)-> Self:
1911        """Set/Get object position along y axis."""
1912        p = self.transform.position
1913        if val is None:
1914            return p[1]
1915        self.pos(p[0], val, p[2])
1916        return self
1917
1918    def z(self, val=None) -> Self:
1919        """Set/Get object position along z axis."""
1920        p = self.transform.position
1921        if val is None:
1922            return p[2]
1923        self.pos(p[0], p[1], val)
1924        return self
1925
1926    def rotate(self, angle: float, axis=(1, 0, 0), point=(0, 0, 0), rad=False) -> Self:
1927        """
1928        Rotate around an arbitrary `axis` passing through `point`.
1929
1930        Example:
1931        ```python
1932        from vedo import *
1933        c1 = Cube()
1934        c2 = c1.clone().c('violet').alpha(0.5) # copy of c1
1935        v = vector(0.2,1,0)
1936        p = vector(1,0,0)  # axis passes through this point
1937        c2.rotate(90, axis=v, point=p)
1938        l = Line(-v+p, v+p).lw(3).c('red')
1939        show(c1, l, c2, axes=1).close()
1940        ```
1941
1942        ![](https://vedo.embl.es/images/feats/rotate_axis.png)
1943        """
1944        LT = LinearTransform()
1945        LT.rotate(angle, axis, point, rad)
1946        return self.apply_transform(LT)
1947
1948    def rotate_x(self, angle: float, rad=False, around=None) -> Self:
1949        """
1950        Rotate around x-axis. If angle is in radians set `rad=True`.
1951
1952        Use `around` to define a pivoting point.
1953        """
1954        if angle == 0:
1955            return self
1956        LT = LinearTransform().rotate_x(angle, rad, around)
1957        return self.apply_transform(LT)
1958
1959    def rotate_y(self, angle: float, rad=False, around=None) -> Self:
1960        """
1961        Rotate around y-axis. If angle is in radians set `rad=True`.
1962
1963        Use `around` to define a pivoting point.
1964        """
1965        if angle == 0:
1966            return self
1967        LT = LinearTransform().rotate_y(angle, rad, around)
1968        return self.apply_transform(LT)
1969
1970    def rotate_z(self, angle: float, rad=False, around=None) -> Self:
1971        """
1972        Rotate around z-axis. If angle is in radians set `rad=True`.
1973
1974        Use `around` to define a pivoting point.
1975        """
1976        if angle == 0:
1977            return self
1978        LT = LinearTransform().rotate_z(angle, rad, around)
1979        return self.apply_transform(LT)
1980
1981    def reorient(self, initaxis, newaxis, rotation=0, rad=False, xyplane=False) -> Self:
1982        """
1983        Reorient the object to point to a new direction from an initial one.
1984        If `initaxis` is None, the object will be assumed in its "default" orientation.
1985        If `xyplane` is True, the object will be rotated to lie on the xy plane.
1986
1987        Use `rotation` to first rotate the object around its `initaxis`.
1988        """
1989        q = self.transform.position
1990        LT = LinearTransform()
1991        LT.reorient(initaxis, newaxis, q, rotation, rad, xyplane)
1992        return self.apply_transform(LT)
1993
1994    def scale(self, s=None, reset=False, origin=True) -> Union[Self, np.array]:
1995        """
1996        Set/get object's scaling factor.
1997
1998        Arguments:
1999            s : (list, float)
2000                scaling factor(s).
2001            reset : (bool)
2002                if True previous scaling factors are ignored.
2003            origin : (bool)
2004                if True scaling is applied with respect to object's position,
2005                otherwise is applied respect to (0,0,0).
2006
2007        Note:
2008            use `s=(sx,sy,sz)` to scale differently in the three coordinates.
2009        """
2010        if s is None:
2011            return np.array(self.transform.T.GetScale())
2012
2013        if not utils.is_sequence(s):
2014            s = [s, s, s]
2015
2016        LT = LinearTransform()
2017        if reset:
2018            old_s = np.array(self.transform.T.GetScale())
2019            LT.scale(s / old_s)
2020        else:
2021            if origin is True:
2022                LT.scale(s, origin=self.transform.position)
2023            elif origin is False:
2024                LT.scale(s, origin=False)
2025            else:
2026                LT.scale(s, origin=origin)
2027
2028        return self.apply_transform(LT)
2029
2030
2031###############################################################################
2032class VolumeAlgorithms(CommonAlgorithms):
2033    """Methods for Volume objects."""
2034
2035    def bounds(self) -> np.ndarray:
2036        """
2037        Get the object bounds.
2038        Returns a list in format `[xmin,xmax, ymin,ymax, zmin,zmax]`.
2039        """
2040        # OVERRIDE CommonAlgorithms.bounds() which is too slow
2041        return np.array(self.dataset.GetBounds())
2042
2043    def isosurface(self, value=None, flying_edges=False) -> "vedo.mesh.Mesh":
2044        """
2045        Return an `Mesh` isosurface extracted from the `Volume` object.
2046
2047        Set `value` as single float or list of values to draw the isosurface(s).
2048        Use flying_edges for faster results (but sometimes can interfere with `smooth()`).
2049
2050        Examples:
2051            - [isosurfaces1.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/isosurfaces1.py)
2052
2053                ![](https://vedo.embl.es/images/volumetric/isosurfaces.png)
2054        """
2055        scrange = self.dataset.GetScalarRange()
2056
2057        if flying_edges:
2058            cf = vtki.new("FlyingEdges3D")
2059            cf.InterpolateAttributesOn()
2060        else:
2061            cf = vtki.new("ContourFilter")
2062            cf.UseScalarTreeOn()
2063
2064        cf.SetInputData(self.dataset)
2065        cf.ComputeNormalsOn()
2066
2067        if utils.is_sequence(value):
2068            cf.SetNumberOfContours(len(value))
2069            for i, t in enumerate(value):
2070                cf.SetValue(i, t)
2071        else:
2072            if value is None:
2073                value = (2 * scrange[0] + scrange[1]) / 3.0
2074                # print("automatic isosurface value =", value)
2075            cf.SetValue(0, value)
2076
2077        cf.Update()
2078        poly = cf.GetOutput()
2079
2080        out = vedo.mesh.Mesh(poly, c=None).phong()
2081        out.mapper.SetScalarRange(scrange[0], scrange[1])
2082
2083        out.pipeline = utils.OperationNode(
2084            "isosurface",
2085            parents=[self],
2086            comment=f"#pts {out.dataset.GetNumberOfPoints()}",
2087            c="#4cc9f0:#e9c46a",
2088        )
2089        return out
2090    
2091    def isosurface_discrete(
2092            self, values, background_label=None, internal_boundaries=True, use_quads=False, nsmooth=0,
2093        ) -> "vedo.mesh.Mesh":
2094        """
2095        Create boundary/isocontour surfaces from a label map (e.g., a segmented image) using a threaded,
2096        3D version of the multiple objects/labels Surface Nets algorithm.
2097        The input is a 3D image (i.e., volume) where each voxel is labeled
2098        (integer labels are preferred to real values), and the output data is a polygonal mesh separating
2099        labeled regions / objects.
2100        (Note that on output each region [corresponding to a different segmented object] will share
2101        points/edges on a common boundary, i.e., two neighboring objects will share the boundary that separates them).
2102        
2103        Besides output geometry defining the surface net, the filter outputs a two-component celldata array indicating
2104        the labels on either side of the polygons composing the output Mesh. 
2105        (This can be used for advanced operations like extracting shared/contacting boundaries between two objects.
2106        The name of this celldata array is "BoundaryLabels").
2107
2108        Arguments:
2109            value : (float, list)
2110                single value or list of values to draw the isosurface(s).
2111            background_label : (float)
2112                this value specifies the label value to use when referencing the background
2113                region outside of any of the specified regions.
2114            boundaries : (bool, list)
2115                if True, the output will only contain the boundary surface. Internal surfaces will be removed.
2116                If a list of integers is provided, only the boundaries between the specified labels will be extracted.
2117            use_quads : (bool)
2118                if True, the output polygons will be quads. If False, the output polygons will be triangles.
2119            nsmooth : (int)
2120                number of iterations of smoothing (0 means no smoothing).
2121
2122        Examples:
2123            - [isosurfaces2.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/isosurfaces2.py)
2124        """
2125        logger = vtki.get_class("Logger")
2126        logger.SetStderrVerbosity(logger.VERBOSITY_ERROR)
2127
2128        snets = vtki.new("SurfaceNets3D")
2129        snets.SetInputData(self.dataset)
2130
2131        if nsmooth:
2132            snets.SmoothingOn()
2133            snets.AutomaticSmoothingConstraintsOn()
2134            snets.GetSmoother().SetNumberOfIterations(nsmooth)
2135            # snets.GetSmoother().SetRelaxationFactor(relaxation_factor)
2136            # snets.GetSmoother().SetConstraintDistance(constraint_distance)
2137        else:
2138            snets.SmoothingOff()
2139
2140        if internal_boundaries is False:
2141            snets.SetOutputStyleToBoundary()
2142        elif internal_boundaries is True:
2143            snets.SetOutputStyleToDefault()
2144        elif utils.is_sequence(internal_boundaries):
2145            snets.SetOutputStyleToSelected()
2146            snets.InitializeSelectedLabelsList()
2147            for val in internal_boundaries:
2148                snets.AddSelectedLabel(val)
2149        else:
2150            vedo.logger.error("isosurface_discrete(): unknown boundaries option")
2151
2152        n = len(values)
2153        snets.SetNumberOfContours(n)
2154        snets.SetNumberOfLabels(n)
2155
2156        if background_label is not None:
2157            snets.SetBackgroundLabel(background_label)
2158
2159        for i, val in enumerate(values):
2160            snets.SetValue(i, val)
2161        
2162        if use_quads:
2163            snets.SetOutputMeshTypeToQuads()
2164        else:
2165            snets.SetOutputMeshTypeToTriangles()
2166        snets.Update()
2167
2168        out = vedo.mesh.Mesh(snets.GetOutput())
2169        out.pipeline = utils.OperationNode(
2170            "isosurface_discrete",
2171            parents=[self],
2172            comment=f"#pts {out.dataset.GetNumberOfPoints()}",
2173            c="#4cc9f0:#e9c46a",
2174        )
2175
2176        logger.SetStderrVerbosity(logger.VERBOSITY_INFO)
2177        return out
2178
2179
2180    def legosurface(
2181        self,
2182        vmin=None,
2183        vmax=None,
2184        invert=False,
2185        boundary=True,
2186        array_name="input_scalars",
2187    ) -> "vedo.mesh.Mesh":
2188        """
2189        Represent an object - typically a `Volume` - as lego blocks (voxels).
2190        By default colors correspond to the volume's scalar.
2191        Returns an `Mesh` object.
2192
2193        Arguments:
2194            vmin : (float)
2195                the lower threshold, voxels below this value are not shown.
2196            vmax : (float)
2197                the upper threshold, voxels above this value are not shown.
2198            boundary : (bool)
2199                controls whether to include cells that are partially inside
2200            array_name : (int, str)
2201                name or index of the scalar array to be considered
2202
2203        Examples:
2204            - [legosurface.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/legosurface.py)
2205
2206                ![](https://vedo.embl.es/images/volumetric/56820682-da40e500-684c-11e9-8ea3-91cbcba24b3a.png)
2207        """
2208        imp_dataset = vtki.new("ImplicitDataSet")
2209        imp_dataset.SetDataSet(self.dataset)
2210        window = vtki.new("ImplicitWindowFunction")
2211        window.SetImplicitFunction(imp_dataset)
2212
2213        srng = list(self.dataset.GetScalarRange())
2214        if vmin is not None:
2215            srng[0] = vmin
2216        if vmax is not None:
2217            srng[1] = vmax
2218        if not boundary:
2219            tol = 0.00001 * (srng[1] - srng[0])
2220            srng[0] -= tol
2221            srng[1] += tol
2222        window.SetWindowRange(srng)
2223        # print("legosurface window range:", srng)
2224
2225        extract = vtki.new("ExtractGeometry")
2226        extract.SetInputData(self.dataset)
2227        extract.SetImplicitFunction(window)
2228        extract.SetExtractInside(invert)
2229        extract.SetExtractBoundaryCells(boundary)
2230        extract.Update()
2231
2232        gf = vtki.new("GeometryFilter")
2233        gf.SetInputData(extract.GetOutput())
2234        gf.Update()
2235
2236        m = vedo.mesh.Mesh(gf.GetOutput()).lw(0.1).flat()
2237        m.map_points_to_cells()
2238        m.celldata.select(array_name)
2239
2240        m.pipeline = utils.OperationNode(
2241            "legosurface",
2242            parents=[self],
2243            comment=f"array: {array_name}",
2244            c="#4cc9f0:#e9c46a",
2245        )
2246        return m
2247
2248    def tomesh(self, fill=True, shrink=1.0) -> "vedo.mesh.Mesh":
2249        """
2250        Build a polygonal Mesh from the current object.
2251
2252        If `fill=True`, the interior faces of all the cells are created.
2253        (setting a `shrink` value slightly smaller than the default 1.0
2254        can avoid flickering due to internal adjacent faces).
2255
2256        If `fill=False`, only the boundary faces will be generated.
2257        """
2258        gf = vtki.new("GeometryFilter")
2259        if fill:
2260            sf = vtki.new("ShrinkFilter")
2261            sf.SetInputData(self.dataset)
2262            sf.SetShrinkFactor(shrink)
2263            sf.Update()
2264            gf.SetInputData(sf.GetOutput())
2265            gf.Update()
2266            poly = gf.GetOutput()
2267            if shrink == 1.0:
2268                clean_poly = vtki.new("CleanPolyData")
2269                clean_poly.PointMergingOn()
2270                clean_poly.ConvertLinesToPointsOn()
2271                clean_poly.ConvertPolysToLinesOn()
2272                clean_poly.ConvertStripsToPolysOn()
2273                clean_poly.SetInputData(poly)
2274                clean_poly.Update()
2275                poly = clean_poly.GetOutput()
2276        else:
2277            gf.SetInputData(self.dataset)
2278            gf.Update()
2279            poly = gf.GetOutput()
2280
2281        msh = vedo.mesh.Mesh(poly).flat()
2282        msh.scalarbar = self.scalarbar
2283        lut = utils.ctf2lut(self)
2284        if lut:
2285            msh.mapper.SetLookupTable(lut)
2286
2287        msh.pipeline = utils.OperationNode(
2288            "tomesh", parents=[self], comment=f"fill={fill}", c="#9e2a2b:#e9c46a"
2289        )
2290        return msh