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