vedo.core

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