vedo.core

   1#!/usr/bin/env python3
   2# -*- coding: utf-8 -*-
   3import numpy as np
   4
   5import vedo.vtkclasses as vtk
   6
   7import vedo
   8from vedo import colors
   9from vedo import utils
  10from vedo.transformations import LinearTransform, NonLinearTransform
  11
  12
  13__docformat__ = "google"
  14
  15__doc__ = """
  16Base classes providing functionality to different vedo objects.
  17
  18![](https://vedo.embl.es/images/feats/algorithms_illustration.png)
  19"""
  20
  21__all__ = [
  22    "DataArrayHelper",
  23    "CommonAlgorithms",
  24    "PointAlgorithms",
  25    "VolumeAlgorithms",
  26]
  27
  28warnings = dict(
  29    points_getter=(
  30        "WARNING: points() is deprecated, use vertices instead. Change:\n"
  31        "         mesh.points() -> mesh.vertices\n"
  32        "         (silence this with vedo.core.warnings['points_getter']=False)"
  33    ),
  34    points_setter=(
  35        "WARNING: points() is deprecated, use vertices instead. Change:\n"
  36        "         mesh.points([[x,y,z], ...]) -> mesh.vertices = [[x,y,z], ...]\n"
  37        "         (silence this with vedo.core.warnings['points_getter']=False)"
  38    ),
  39)
  40
  41###############################################################################
  42class DataArrayHelper:
  43    # Internal use only.
  44    # Helper class to manage data associated to either
  45    # points (or vertices) and cells (or faces).
  46    def __init__(self, obj, association):
  47
  48        self.obj = obj
  49        self.association = association
  50
  51    def __getitem__(self, key):
  52
  53        if self.association == 0:
  54            data = self.obj.dataset.GetPointData()
  55
  56        elif self.association == 1:
  57            data = self.obj.dataset.GetCellData()
  58
  59        elif self.association == 2:
  60            data = self.obj.dataset.GetFieldData()
  61
  62            varr = data.GetAbstractArray(key)
  63            if isinstance(varr, vtk.vtkStringArray):
  64                if isinstance(key, int):
  65                    key = data.GetArrayName(key)
  66                n = varr.GetNumberOfValues()
  67                narr = [varr.GetValue(i) for i in range(n)]
  68                return narr
  69                ###########
  70
  71        else:
  72            raise RuntimeError()
  73
  74        if isinstance(key, int):
  75            key = data.GetArrayName(key)
  76
  77        arr = data.GetArray(key)
  78        if not arr:
  79            return None
  80        return utils.vtk2numpy(arr)
  81
  82    def __setitem__(self, key, input_array):
  83
  84        if self.association == 0:
  85            data = self.obj.dataset.GetPointData()
  86            n = self.obj.dataset.GetNumberOfPoints()
  87            self.obj.mapper.SetScalarModeToUsePointData()
  88
  89        elif self.association == 1:
  90            data = self.obj.dataset.GetCellData()
  91            n = self.obj.dataset.GetNumberOfCells()
  92            self.obj.mapper.SetScalarModeToUseCellData()
  93
  94        elif self.association == 2:
  95            data = self.obj.dataset.GetFieldData()
  96            if not utils.is_sequence(input_array):
  97                input_array = [input_array]
  98
  99            if isinstance(input_array[0], str):
 100                varr = vtk.vtkStringArray()
 101                varr.SetName(key)
 102                varr.SetNumberOfComponents(1)
 103                varr.SetNumberOfTuples(len(input_array))
 104                for i, iarr in enumerate(input_array):
 105                    if isinstance(iarr, np.ndarray):
 106                        iarr = iarr.tolist()  # better format
 107                        # Note: a string k can be converted to numpy with
 108                        # import json; k = np.array(json.loads(k))
 109                    varr.InsertValue(i, str(iarr))
 110            else:
 111                try:
 112                    varr = utils.numpy2vtk(input_array, name=key)
 113                except TypeError as e:
 114                    vedo.logger.error(
 115                        f"cannot create metadata with input object:\n"
 116                        f"{input_array}"
 117                        f"\n\nAllowed content examples are:\n"
 118                        f"- flat list of strings ['a','b', 1, [1,2,3], ...]"
 119                        f" (first item must be a string in this case)\n"
 120                        f"  hint: use k = np.array(json.loads(k)) to convert strings\n"
 121                        f"- numpy arrays of any shape"
 122                    )
 123                    raise e
 124
 125            data.AddArray(varr)
 126            return  ############
 127
 128        else:
 129            raise RuntimeError()
 130
 131        if len(input_array) != n:
 132            vedo.logger.error(
 133                f"Error in point/cell data: length of input {len(input_array)}"
 134                f" !=  {n} nr. of elements"
 135            )
 136            raise RuntimeError()
 137
 138        input_array = np.asarray(input_array)
 139        varr = utils.numpy2vtk(input_array, name=key)
 140        data.AddArray(varr)
 141
 142        if len(input_array.shape) == 1:  # scalars
 143            data.SetActiveScalars(key)
 144        elif len(input_array.shape) == 2 and input_array.shape[1] == 3:  # vectors
 145            if key.lower() == "normals":
 146                data.SetActiveNormals(key)
 147            else:
 148                data.SetActiveVectors(key)
 149
 150    def keys(self):
 151        """Return the list of available data array names"""
 152        if self.association == 0:
 153            data = self.obj.dataset.GetPointData()
 154        elif self.association == 1:
 155            data = self.obj.dataset.GetCellData()
 156        elif self.association == 2:
 157            data = self.obj.dataset.GetFieldData()
 158        arrnames = []
 159        for i in range(data.GetNumberOfArrays()):
 160            if self.association == 2:
 161                name = data.GetAbstractArray(i).GetName()
 162            else:
 163                name = data.GetArray(i).GetName()
 164            if name:
 165                arrnames.append(name)
 166        return arrnames
 167    
 168    def items(self):
 169        """Return the list of available data array `(names, values)`."""
 170        if self.association == 0:
 171            data = self.obj.dataset.GetPointData()
 172        elif self.association == 1:
 173            data = self.obj.dataset.GetCellData()
 174        elif self.association == 2:
 175            data = self.obj.dataset.GetFieldData()
 176        arrnames = []
 177        for i in range(data.GetNumberOfArrays()):
 178            if self.association == 2:
 179                name = data.GetAbstractArray(i).GetName()
 180            else:
 181                name = data.GetArray(i).GetName()
 182            if name:
 183                arrnames.append((name, self[name]))
 184        return arrnames
 185    
 186    def todict(self):
 187        """Return a dictionary of the available data arrays."""
 188        return dict(self.items())
 189
 190    def rename(self, oldname, newname):
 191        """Rename an array"""
 192        if self.association == 0:
 193            varr = self.obj.dataset.GetPointData().GetArray(oldname)
 194        elif self.association == 1:
 195            varr = self.obj.dataset.GetCellData().GetArray(oldname)
 196        elif self.association == 2:
 197            varr = self.obj.dataset.GetFieldData().GetAbstractArray(oldname)
 198        if varr:
 199            varr.SetName(newname)
 200        else:
 201            vedo.logger.warning(
 202                f"Cannot rename non existing array {oldname} to {newname}"
 203            )
 204
 205    def remove(self, key):
 206        """Remove a data array by name or number"""
 207        if self.association == 0:
 208            self.obj.dataset.GetPointData().RemoveArray(key)
 209        elif self.association == 1:
 210            self.obj.dataset.GetCellData().RemoveArray(key)
 211        elif self.association == 2:
 212            self.obj.dataset.GetFieldData().RemoveArray(key)
 213
 214    def clear(self):
 215        """Remove all data associated to this object"""
 216        if self.association == 0:
 217            data = self.obj.dataset.GetPointData()
 218        elif self.association == 1:
 219            data = self.obj.dataset.GetCellData()
 220        elif self.association == 2:
 221            data = self.obj.dataset.GetFieldData()
 222        for i in range(data.GetNumberOfArrays()):
 223            if self.association == 2:
 224                name = data.GetAbstractArray(i).GetName()
 225            else:
 226                name = data.GetArray(i).GetName()
 227            data.RemoveArray(name)
 228
 229    def select(self, key):
 230        """Select one specific array by its name to make it the `active` one."""
 231        # Default (ColorModeToDefault): unsigned char scalars are treated as colors,
 232        # and NOT mapped through the lookup table, while everything else is.
 233        # ColorModeToDirectScalar extends ColorModeToDefault such that all integer
 234        # types are treated as colors with values in the range 0-255
 235        # and floating types are treated as colors with values in the range 0.0-1.0.
 236        # Setting ColorModeToMapScalars means that all scalar data will be mapped
 237        # through the lookup table.
 238        # (Note that for multi-component scalars, the particular component 
 239        # to use for mapping can be specified using the SelectColorArray() method.)
 240        if self.association == 0:
 241            data = self.obj.dataset.GetPointData()
 242            self.obj.mapper.SetScalarModeToUsePointData()
 243        else:
 244            data = self.obj.dataset.GetCellData()
 245            self.obj.mapper.SetScalarModeToUseCellData()
 246
 247        if isinstance(key, int):
 248            key = data.GetArrayName(key)
 249
 250        arr = data.GetArray(key)
 251        if not arr:
 252            return self.obj
 253
 254        nc = arr.GetNumberOfComponents()
 255        # print("GetNumberOfComponents", nc)
 256        if nc == 1:
 257            data.SetActiveScalars(key)
 258        elif nc == 2:
 259            data.SetTCoords(arr)
 260        elif nc == 3 or nc == 4:
 261            if "rgb" in key.lower():
 262                data.SetActiveScalars(key)
 263                try:
 264                    # could be a volume mapper
 265                    self.obj.mapper.SetColorModeToDirectScalars()
 266                except AttributeError:
 267                    pass
 268            else:
 269                data.SetActiveVectors(key)
 270        elif nc == 9:
 271            data.SetActiveTensors(key)
 272        else:
 273            vedo.logger.error(f"Cannot select array {key} with {nc} components")
 274            return self.obj
 275
 276        try:
 277            # could be a volume mapper
 278            self.obj.mapper.SetArrayName(key)
 279            self.obj.mapper.ScalarVisibilityOn()
 280        except AttributeError:
 281            pass
 282
 283        return self.obj
 284    
 285    def select_texture_coords(self, key):
 286        """Select one specific array to be used as texture coordinates."""
 287        if self.association == 0:
 288            data = self.obj.dataset.GetPointData()
 289        else:
 290            vedo.logger.warning("texture coordinates are only available for point data")
 291            return
 292
 293        if isinstance(key, int):
 294            key = data.GetArrayName(key)
 295        data.SetTCoords(data.GetArray(key))
 296        return self.obj
 297
 298    def select_normals(self, key):
 299        """Select one specific normal array by its name to make it the "active" one."""
 300        if self.association == 0:
 301            data = self.obj.dataset.GetPointData()
 302            self.obj.mapper.SetScalarModeToUsePointData()
 303        else:
 304            data = self.obj.dataset.GetCellData()
 305            self.obj.mapper.SetScalarModeToUseCellData()
 306
 307        if isinstance(key, int):
 308            key = data.GetArrayName(key)
 309        data.SetActiveNormals(key)
 310        return self.obj
 311
 312    def print(self, **kwargs):
 313        """Print the array names available to terminal"""
 314        colors.printc(self.keys(), **kwargs)
 315
 316    def __repr__(self) -> str:
 317        """Representation"""
 318
 319        def _get_str(pd, header):
 320            out = f"\x1b[2m\x1b[1m\x1b[7m{header}"
 321            if pd.GetNumberOfArrays():
 322                if self.obj.name:
 323                    out += f" in {self.obj.name}"
 324                out += f" contains {pd.GetNumberOfArrays()} array(s)\x1b[0m"
 325                for i in range(pd.GetNumberOfArrays()):
 326                    varr = pd.GetArray(i)
 327                    out += f"\n\x1b[1m\x1b[4mArray name    : {varr.GetName()}\x1b[0m"
 328                    out += "\nindex".ljust(15) + f": {i}"
 329                    t = varr.GetDataType()
 330                    if t in vedo.utils.array_types:
 331                        out += "\ntype".ljust(15)
 332                        out += f": {vedo.utils.array_types[t][1]} ({vedo.utils.array_types[t][0]})"
 333                    shape = (varr.GetNumberOfTuples(), varr.GetNumberOfComponents())
 334                    out += "\nshape".ljust(15) + f": {shape}"
 335                    out += "\nrange".ljust(15) + f": {np.array(varr.GetRange())}"
 336                    out += "\nmax id".ljust(15) + f": {varr.GetMaxId()}"
 337                    out += "\nlook up table".ljust(15) + f": {bool(varr.GetLookupTable())}"
 338                    out += "\nin-memory size".ljust(15) + f": {varr.GetActualMemorySize()} KB"
 339            else:
 340                out += " is empty.\x1b[0m"
 341            return out
 342
 343        if self.association == 0:
 344            out = _get_str(self.obj.dataset.GetPointData(), "Point Data")
 345        elif self.association == 1:
 346            out = _get_str(self.obj.dataset.GetCellData(), "Cell Data")
 347        elif self.association == 2:
 348            pd = self.obj.dataset.GetFieldData()
 349            if pd.GetNumberOfArrays():
 350                out = "\x1b[2m\x1b[1m\x1b[7mMeta Data"
 351                if self.obj.name:
 352                    out += f" in {self.obj.name}"
 353                out += f" contains {pd.GetNumberOfArrays()} entries\x1b[0m"
 354                for i in range(pd.GetNumberOfArrays()):
 355                    varr = pd.GetAbstractArray(i)
 356                    out += f"\n\x1b[1m\x1b[4mEntry name    : {varr.GetName()}\x1b[0m"
 357                    out += "\nindex".ljust(15) + f": {i}"
 358                    shape = (varr.GetNumberOfTuples(), varr.GetNumberOfComponents())
 359                    out += "\nshape".ljust(15) + f": {shape}"
 360
 361        return out
 362
 363
 364###############################################################################
 365class CommonAlgorithms:
 366    """Common algorithms."""
 367
 368    def __init__(self):
 369        # print("init CommonAlgorithms")
 370        self.dataset = None
 371        self.pipeline = None
 372        self.name = ""
 373        self.filename = ""
 374        self.time = 0
 375
 376
 377    @property
 378    def pointdata(self):
 379        """
 380        Create and/or return a `numpy.array` associated to points (vertices).
 381        A data array can be indexed either as a string or by an integer number.
 382        E.g.:  `myobj.pointdata["arrayname"]`
 383
 384        Usage:
 385
 386            `myobj.pointdata.keys()` to return the available data array names
 387
 388            `myobj.pointdata.select(name)` to make this array the active one
 389
 390            `myobj.pointdata.remove(name)` to remove this array
 391        """
 392        return DataArrayHelper(self, 0)
 393
 394    @property
 395    def celldata(self):
 396        """
 397        Create and/or return a `numpy.array` associated to cells (faces).
 398        A data array can be indexed either as a string or by an integer number.
 399        E.g.:  `myobj.celldata["arrayname"]`
 400
 401        Usage:
 402
 403            `myobj.celldata.keys()` to return the available data array names
 404
 405            `myobj.celldata.select(name)` to make this array the active one
 406
 407            `myobj.celldata.remove(name)` to remove this array
 408        """
 409        return DataArrayHelper(self, 1)
 410
 411    @property
 412    def metadata(self):
 413        """
 414        Create and/or return a `numpy.array` associated to neither cells nor faces.
 415        A data array can be indexed either as a string or by an integer number.
 416        E.g.:  `myobj.metadata["arrayname"]`
 417
 418        Usage:
 419
 420            `myobj.metadata.keys()` to return the available data array names
 421
 422            `myobj.metadata.select(name)` to make this array the active one
 423
 424            `myobj.metadata.remove(name)` to remove this array
 425        """
 426        return DataArrayHelper(self, 2)
 427
 428    def memory_address(self):
 429        """
 430        Return a unique memory address integer which may serve as the ID of the
 431        object, or passed to c++ code.
 432        """
 433        # https://www.linkedin.com/pulse/speedup-your-code-accessing-python-vtk-objects-from-c-pletzer/
 434        # https://github.com/tfmoraes/polydata_connectivity
 435        return int(self.dataset.GetAddressAsString("")[5:], 16)
 436
 437    def memory_size(self):
 438        """Return the size in bytes of the object in memory."""
 439        return self.dataset.GetActualMemorySize()
 440
 441    def modified(self):
 442        """Use in conjunction with `tonumpy()` to update any modifications to the image array."""
 443        self.dataset.GetPointData().Modified()
 444        self.dataset.GetPointData().GetScalars().Modified()
 445        return self
 446
 447    def box(self, scale=1, padding=0):
 448        """
 449        Return the bounding box as a new `Mesh` object.
 450
 451        Arguments:
 452            scale : (float)
 453                box size can be scaled by a factor
 454            padding : (float, list)
 455                a constant padding can be added (can be a list `[padx,pady,padz]`)
 456        """
 457        b = self.bounds()
 458        if not utils.is_sequence(padding):
 459            padding = [padding, padding, padding]
 460        length, width, height = b[1] - b[0], b[3] - b[2], b[5] - b[4]
 461        tol = (length + width + height) / 30000  # useful for boxing text
 462        pos = [(b[0] + b[1])/2, (b[3] + b[2])/2, (b[5] + b[4])/2 - tol]
 463        bx = vedo.shapes.Box(
 464            pos,
 465            length * scale + padding[0],
 466            width  * scale + padding[1],
 467            height * scale + padding[2],
 468            c="gray",
 469        )
 470        try:
 471            pr = vtk.vtkProperty()
 472            pr.DeepCopy(self.properties)
 473            bx.SetProperty(pr)
 474            bx.properties = pr
 475        except (AttributeError, TypeError):
 476            pass
 477        bx.flat().lighting("off").wireframe(True)
 478        return bx
 479
 480    def bounds(self):
 481        """
 482        Get the object bounds.
 483        Returns a list in format `[xmin,xmax, ymin,ymax, zmin,zmax]`.
 484        """
 485        try: # this is very slow for large meshes
 486            pts = self.vertices
 487            xmin, ymin, zmin = np.min(pts, axis=0)
 488            xmax, ymax, zmax = np.max(pts, axis=0)
 489            return np.array([xmin, xmax, ymin, ymax, zmin, zmax])
 490        except (AttributeError, ValueError):
 491            return np.array(self.dataset.GetBounds())
 492
 493    def xbounds(self, i=None):
 494        """Get the bounds `[xmin,xmax]`. Can specify upper or lower with i (0,1)."""
 495        b = self.bounds()
 496        if i is not None:
 497            return b[i]
 498        return np.array([b[0], b[1]])
 499
 500    def ybounds(self, i=None):
 501        """Get the bounds `[ymin,ymax]`. Can specify upper or lower with i (0,1)."""
 502        b = self.bounds()
 503        if i == 0:
 504            return b[2]
 505        if i == 1:
 506            return b[3]
 507        return np.array([b[2], b[3]])
 508
 509    def zbounds(self, i=None):
 510        """Get the bounds `[zmin,zmax]`. Can specify upper or lower with i (0,1)."""
 511        b = self.bounds()
 512        if i == 0:
 513            return b[4]
 514        if i == 1:
 515            return b[5]
 516        return np.array([b[4], b[5]])
 517
 518    def diagonal_size(self):
 519        """Get the length of the diagonal of the bounding box."""
 520        b = self.bounds()
 521        return np.sqrt((b[1] - b[0])**2 + (b[3] - b[2])**2 + (b[5] - b[4])**2)
 522
 523    def average_size(self):
 524        """
 525        Calculate and return the average size of the object.
 526        This is the mean of the vertex distances from the center of mass.
 527        """
 528        coords = self.vertices
 529        cm = np.mean(coords, axis=0)
 530        if coords.shape[0] == 0:
 531            return 0.0
 532        cc = coords - cm
 533        return np.mean(np.linalg.norm(cc, axis=1))
 534
 535    def center_of_mass(self):
 536        """Get the center of mass of the object."""
 537        cmf = vtk.new("CenterOfMass")
 538        cmf.SetInputData(self.dataset)
 539        cmf.Update()
 540        c = cmf.GetCenter()
 541        return np.array(c)
 542
 543    def copy_data_from(self, obj):
 544        """Copy all data (point and cell data) from this input object"""
 545        self.dataset.GetPointData().PassData(obj.dataset.GetPointData())
 546        self.dataset.GetCellData().PassData(obj.dataset.GetCellData())
 547        self.pipeline = utils.OperationNode(
 548            "copy_data_from",
 549            parents=[self, obj],
 550            comment=f"{obj.__class__.__name__}",
 551            shape="note",
 552            c="#ccc5b9",
 553        )
 554        return self
 555
 556    def inputdata(self):
 557        """Obsolete, use `.dataset` instead."""
 558        colors.printc("WARNING: 'inputdata()' is obsolete, use '.dataset' instead.", c="y")
 559        return self.dataset
 560
 561    @property
 562    def npoints(self):
 563        """Retrieve the number of points (or vertices)."""
 564        return self.dataset.GetNumberOfPoints()
 565
 566    @property
 567    def nvertices(self):
 568        """Retrieve the number of vertices (or points)."""
 569        return self.dataset.GetNumberOfPoints()
 570
 571    @property
 572    def ncells(self):
 573        """Retrieve the number of cells."""
 574        return self.dataset.GetNumberOfCells()
 575
 576    def points(self, pts=None):
 577        """Obsolete, use `self.vertices` or `self.coordinates` instead."""
 578        if pts is None:  ### getter
 579
 580            if warnings["points_getter"]:
 581                colors.printc(warnings["points_getter"], c="y")
 582                warnings["points_getter"] = ""
 583            return self.vertices
 584
 585        else:  ### setter
 586
 587            if warnings["points_setter"]:
 588                colors.printc(warnings["points_setter"], c="y")
 589                warnings["points_setter"] = ""
 590
 591            pts = np.asarray(pts, dtype=np.float32)
 592
 593            if pts.ndim == 1:
 594                ### getter by point index ###################
 595                indices = pts.astype(int)
 596                vpts = self.dataset.GetPoints()
 597                arr = utils.vtk2numpy(vpts.GetData())
 598                return arr[indices]  ###########
 599
 600            ### setter ####################################
 601            if pts.shape[1] == 2:
 602                pts = np.c_[pts, np.zeros(pts.shape[0], dtype=np.float32)]
 603            arr = utils.numpy2vtk(pts, dtype=np.float32)
 604
 605            vpts = self.dataset.GetPoints()
 606            vpts.SetData(arr)
 607            vpts.Modified()
 608            # reset mesh to identity matrix position/rotation:
 609            self.point_locator = None
 610            self.cell_locator = None
 611            self.line_locator = None
 612            self.transform = LinearTransform()
 613            return self
 614
 615    @property
 616    def cell_centers(self):
 617        """
 618        Get the coordinates of the cell centers.
 619
 620        Examples:
 621            - [delaunay2d.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/delaunay2d.py)
 622        """
 623        vcen = vtk.new("CellCenters")
 624        vcen.SetInputData(self.dataset)
 625        vcen.Update()
 626        return utils.vtk2numpy(vcen.GetOutput().GetPoints().GetData())
 627
 628    @property
 629    def lines(self):
 630        """
 631        Get lines connectivity ids as a python array
 632        formatted as `[[id0,id1], [id3,id4], ...]`
 633
 634        See also: `lines_as_flat_array()`.
 635        """
 636        # Get cell connettivity ids as a 1D array. The vtk format is:
 637        #    [nids1, id0 ... idn, niids2, id0 ... idm,  etc].
 638        arr1d = utils.vtk2numpy(self.dataset.GetLines().GetData())
 639        i = 0
 640        conn = []
 641        n = len(arr1d)
 642        for _ in range(n):
 643            cell = [arr1d[i + k + 1] for k in range(arr1d[i])]
 644            conn.append(cell)
 645            i += arr1d[i] + 1
 646            if i >= n:
 647                break
 648
 649        return conn  # cannot always make a numpy array of it!
 650
 651    @property
 652    def lines_as_flat_array(self):
 653        """
 654        Get lines connectivity ids as a 1D numpy array.
 655        Format is e.g. [2,  10,20,  3, 10,11,12,  2, 70,80, ...]
 656
 657        See also: `lines()`.
 658        """
 659        return utils.vtk2numpy(self.dataset.GetLines().GetData())
 660
 661    def mark_boundaries(self):
 662        """
 663        Mark cells and vertices if they lie on a boundary.
 664        A new array called `BoundaryCells` is added to the object.
 665        """
 666        mb = vtk.new("MarkBoundaryFilter")
 667        mb.SetInputData(self.dataset)
 668        mb.Update()
 669        self.dataset.DeepCopy(mb.GetOutput())
 670        self.pipeline = utils.OperationNode("mark_boundaries", parents=[self])
 671        return self
 672
 673    def find_cells_in_bounds(self, xbounds=(), ybounds=(), zbounds=()):
 674        """
 675        Find cells that are within the specified bounds.
 676        """
 677        try:
 678            xbounds = list(xbounds.bounds())
 679        except AttributeError:
 680            pass
 681
 682        if len(xbounds) == 6:
 683            bnds = xbounds
 684        else:
 685            bnds = list(self.bounds())
 686            if len(xbounds) == 2:
 687                bnds[0] = xbounds[0]
 688                bnds[1] = xbounds[1]
 689            if len(ybounds) == 2:
 690                bnds[2] = ybounds[0]
 691                bnds[3] = ybounds[1]
 692            if len(zbounds) == 2:
 693                bnds[4] = zbounds[0]
 694                bnds[5] = zbounds[1]
 695
 696        cell_ids = vtk.vtkIdList()
 697        if not self.cell_locator:
 698            self.cell_locator = vtk.new("CellTreeLocator")
 699            self.cell_locator.SetDataSet(self.dataset)
 700            self.cell_locator.BuildLocator()
 701        self.cell_locator.FindCellsWithinBounds(bnds, cell_ids)
 702        cids = []
 703        for i in range(cell_ids.GetNumberOfIds()):
 704            cid = cell_ids.GetId(i)
 705            cids.append(cid)
 706        return np.array(cids)
 707
 708    def find_cells_along_line(self, p0, p1, tol=0.001):
 709        """
 710        Find cells that are intersected by a line segment.
 711        """
 712        cell_ids = vtk.vtkIdList()
 713        if not self.cell_locator:
 714            self.cell_locator = vtk.new("CellTreeLocator")
 715            self.cell_locator.SetDataSet(self.dataset)
 716            self.cell_locator.BuildLocator()
 717        self.cell_locator.FindCellsAlongLine(p0, p1, tol, cell_ids)   
 718        cids = []
 719        for i in range(cell_ids.GetNumberOfIds()):
 720            cid = cell_ids.GetId(i)
 721            cids.append(cid)
 722        return np.array(cids)
 723
 724    def find_cells_along_plane(self, origin, normal, tol=0.001):
 725        """
 726        Find cells that are intersected by a plane.
 727        """
 728        cell_ids = vtk.vtkIdList()
 729        if not self.cell_locator:
 730            self.cell_locator = vtk.new("CellTreeLocator")
 731            self.cell_locator.SetDataSet(self.dataset)
 732            self.cell_locator.BuildLocator()
 733        self.cell_locator.FindCellsAlongPlane(origin, normal, tol, cell_ids)   
 734        cids = []
 735        for i in range(cell_ids.GetNumberOfIds()):
 736            cid = cell_ids.GetId(i)
 737            cids.append(cid)
 738        return np.array(cids)
 739
 740    def delete_cells_by_point_index(self, indices):
 741        """
 742        Delete a list of vertices identified by any of their vertex index.
 743
 744        See also `delete_cells()`.
 745
 746        Examples:
 747            - [delete_mesh_pts.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/delete_mesh_pts.py)
 748
 749                ![](https://vedo.embl.es/images/basic/deleteMeshPoints.png)
 750        """
 751        cell_ids = vtk.vtkIdList()
 752        self.dataset.BuildLinks()
 753        n = 0
 754        for i in np.unique(indices):
 755            self.dataset.GetPointCells(i, cell_ids)
 756            for j in range(cell_ids.GetNumberOfIds()):
 757                self.dataset.DeleteCell(cell_ids.GetId(j))  # flag cell
 758                n += 1
 759
 760        self.dataset.RemoveDeletedCells()
 761        self.dataset.Modified()
 762        self.pipeline = utils.OperationNode(
 763            "delete_cells_by_point_index", parents=[self])
 764        return self
 765
 766    def map_cells_to_points(self, arrays=(), move=False):
 767        """
 768        Interpolate cell data (i.e., data specified per cell or face)
 769        into point data (i.e., data specified at each vertex).
 770        The method of transformation is based on averaging the data values
 771        of all cells using a particular point.
 772
 773        A custom list of arrays to be mapped can be passed in input.
 774
 775        Set `move=True` to delete the original `celldata` array.
 776        """
 777        c2p = vtk.new("CellDataToPointData")
 778        c2p.SetInputData(self.dataset)
 779        if not move:
 780            c2p.PassCellDataOn()
 781        if arrays:
 782            c2p.ClearCellDataArrays()
 783            c2p.ProcessAllArraysOff()
 784            for arr in arrays:
 785                c2p.AddCellDataArray(arr)
 786        else:
 787            c2p.ProcessAllArraysOn()
 788        c2p.Update()
 789        self._update(c2p.GetOutput(), reset_locators=False)
 790        self.mapper.SetScalarModeToUsePointData()
 791        self.pipeline = utils.OperationNode("map_cells_to_points", parents=[self])
 792        return self
 793
 794    @property
 795    def vertices(self):
 796        """Return the vertices (points) coordinates."""
 797        try:
 798            # valid for polydata and unstructured grid
 799            varr = self.dataset.GetPoints().GetData()
 800
 801        except AttributeError:
 802            try:
 803                # valid for rectilinear/structured grid, image data
 804                v2p = vtk.new("ImageToPoints")
 805                v2p.SetInputData(self.dataset)
 806                v2p.Update()
 807                varr = v2p.GetOutput().GetPoints().GetData()
 808            except AttributeError:
 809                return np.array([])
 810        
 811        narr = utils.vtk2numpy(varr)
 812        return narr
 813
 814    # setter
 815    @vertices.setter
 816    def vertices(self, pts):
 817        """Set vertices (points) coordinates."""
 818        pts = utils.make3d(pts)
 819        arr = utils.numpy2vtk(pts, dtype=np.float32)
 820        try:
 821            vpts = self.dataset.GetPoints()
 822            vpts.SetData(arr)
 823            vpts.Modified()
 824        except AttributeError:
 825            vedo.logger.error(f"Cannot set vertices for object {type(self)}")
 826            return self
 827        # reset mesh to identity matrix position/rotation:
 828        self.point_locator = None
 829        self.cell_locator = None
 830        self.line_locator = None
 831        self.transform = LinearTransform()
 832        return self
 833
 834
 835    @property
 836    def coordinates(self):
 837        """Return the vertices (points) coordinates. Same as `vertices`."""
 838        return self.vertices
 839
 840    @coordinates.setter
 841    def coordinates(self, pts):
 842        """Set vertices (points) coordinates. Same as `vertices`."""
 843        self.vertices = pts
 844    
 845    @property
 846    def cells_as_flat_array(self):
 847        """
 848        Get cell connectivity ids as a 1D numpy array.
 849        Format is e.g. [3,  10,20,30  4, 10,11,12,13  ...]
 850        """
 851        try:
 852            # valid for unstructured grid
 853            arr1d = utils.vtk2numpy(self.dataset.GetCells().GetData())
 854        except AttributeError:
 855            # valid for polydata
 856            arr1d = utils.vtk2numpy(self.dataset.GetPolys().GetData())
 857            # if arr1d.size == 0:
 858            #     arr1d = utils.vtk2numpy(self.dataset.GetStrips().GetData())
 859        return arr1d
 860
 861    @property
 862    def cells(self):
 863        """
 864        Get the cells connectivity ids as a numpy array.
 865
 866        The output format is: `[[id0 ... idn], [id0 ... idm],  etc]`.
 867        """
 868        try:
 869            # valid for unstructured grid
 870            arr1d = utils.vtk2numpy(self.dataset.GetCells().GetData())
 871        except AttributeError:
 872            # valid for polydata
 873            arr1d = utils.vtk2numpy(self.dataset.GetPolys().GetData())
 874            # if arr1d.size == 0:
 875            #     arr1d = utils.vtk2numpy(self.dataset.GetStrips().GetData())
 876
 877        # Get cell connettivity ids as a 1D array. vtk format is:
 878        # [nids1, id0 ... idn, niids2, id0 ... idm,  etc].
 879        i = 0
 880        conn = []
 881        n = len(arr1d)
 882        if n:
 883            while True:
 884                cell = [arr1d[i + k] for k in range(1, arr1d[i] + 1)]
 885                conn.append(cell)
 886                i += arr1d[i] + 1
 887                if i >= n:
 888                    break
 889        return conn
 890
 891    def map_points_to_cells(self, arrays=(), move=False):
 892        """
 893        Interpolate point data (i.e., data specified per point or vertex)
 894        into cell data (i.e., data specified per cell).
 895        The method of transformation is based on averaging the data values
 896        of all points defining a particular cell.
 897
 898        A custom list of arrays to be mapped can be passed in input.
 899
 900        Set `move=True` to delete the original `pointdata` array.
 901
 902        Examples:
 903            - [mesh_map2cell.py](https://github.com/marcomusy/vedo/tree/master/examples/basic/mesh_map2cell.py)
 904        """
 905        p2c = vtk.new("PointDataToCellData")
 906        p2c.SetInputData(self.dataset)
 907        if not move:
 908            p2c.PassPointDataOn()
 909        if arrays:
 910            p2c.ClearPointDataArrays()
 911            p2c.ProcessAllArraysOff()
 912            for arr in arrays:
 913                p2c.AddPointDataArray(arr)
 914        else:
 915            p2c.ProcessAllArraysOn()
 916        p2c.Update()
 917        self._update(p2c.GetOutput(), reset_locators=False)
 918        self.mapper.SetScalarModeToUseCellData()
 919        self.pipeline = utils.OperationNode("map_points_to_cells", parents=[self])
 920        return self
 921
 922    def resample_data_from(self, source, tol=None, categorical=False):
 923        """
 924        Resample point and cell data from another dataset.
 925        The output has the same structure but its point data have
 926        the resampled values from target.
 927
 928        Use `tol` to set the tolerance used to compute whether
 929        a point in the source is in a cell of the current object.
 930        Points without resampled values, and their cells, are marked as blank.
 931        If the data is categorical, then the resulting data will be determined
 932        by a nearest neighbor interpolation scheme.
 933
 934        Example:
 935        ```python
 936        from vedo import *
 937        m1 = Mesh(dataurl+'bunny.obj')#.add_gaussian_noise(0.1)
 938        pts = m1.vertices
 939        ces = m1.cell_centers
 940        m1.pointdata["xvalues"] = np.power(pts[:,0], 3)
 941        m1.celldata["yvalues"]  = np.power(ces[:,1], 3)
 942        m2 = Mesh(dataurl+'bunny.obj')
 943        m2.resample_arrays_from(m1)
 944        # print(m2.pointdata["xvalues"])
 945        show(m1, m2 , N=2, axes=1)
 946        ```
 947        """
 948        rs = vtk.new("ResampleWithDataSet")
 949        rs.SetInputData(self.dataset)
 950        rs.SetSourceData(source.dataset)
 951
 952        rs.SetPassPointArrays(True)
 953        rs.SetPassCellArrays(True)
 954        rs.SetPassFieldArrays(True)
 955        rs.SetCategoricalData(categorical)
 956
 957        rs.SetComputeTolerance(True)
 958        if tol:
 959            rs.SetComputeTolerance(False)
 960            rs.SetTolerance(tol)
 961        rs.Update()
 962        self._update(rs.GetOutput(), reset_locators=False)
 963        self.pipeline = utils.OperationNode(
 964            "resample_data_from",
 965            comment=f"{source.__class__.__name__}",
 966            parents=[self, source]
 967        )
 968        return self
 969
 970    def interpolate_data_from(
 971        self,
 972        source,
 973        radius=None,
 974        n=None,
 975        kernel="shepard",
 976        exclude=("Normals",),
 977        on="points",
 978        null_strategy=1,
 979        null_value=0,
 980    ):
 981        """
 982        Interpolate over source to port its data onto the current object using various kernels.
 983
 984        If n (number of closest points to use) is set then radius value is ignored.
 985
 986        Check out also:
 987            `probe()` which in many cases can be faster.
 988
 989        Arguments:
 990            kernel : (str)
 991                available kernels are [shepard, gaussian, linear]
 992            null_strategy : (int)
 993                specify a strategy to use when encountering a "null" point
 994                during the interpolation process. Null points occur when the local neighborhood
 995                (of nearby points to interpolate from) is empty.
 996
 997                - Case 0: an output array is created that marks points
 998                  as being valid (=1) or null (invalid =0), and the null_value is set as well
 999                - Case 1: the output data value(s) are set to the provided null_value
1000                - Case 2: simply use the closest point to perform the interpolation.
1001            null_value : (float)
1002                see above.
1003
1004        Examples:
1005            - [interpolate_scalar1.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/interpolate_scalar1.py)
1006            - [interpolate_scalar3.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/interpolate_scalar3.py)
1007            - [interpolate_scalar4.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/interpolate_scalar4.py)
1008            - [image_probe.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/image_probe.py)
1009
1010                ![](https://vedo.embl.es/images/advanced/interpolateMeshArray.png)
1011        """
1012        if radius is None and not n:
1013            vedo.logger.error("in interpolate_data_from(): please set either radius or n")
1014            raise RuntimeError
1015
1016        if on == "points":
1017            points = source.dataset
1018        elif on == "cells":
1019            c2p = vtk.new("CellDataToPointData")
1020            c2p.SetInputData(source.dataset)
1021            c2p.Update()
1022            points = c2p.GetOutput()
1023        else:
1024            vedo.logger.error("in interpolate_data_from(), on must be on points or cells")
1025            raise RuntimeError()
1026
1027        locator = vtk.new("PointLocator")
1028        locator.SetDataSet(points)
1029        locator.BuildLocator()
1030
1031        if kernel.lower() == "shepard":
1032            kern = vtk.new("ShepardKernel")
1033            kern.SetPowerParameter(2)
1034        elif kernel.lower() == "gaussian":
1035            kern = vtk.new("GaussianKernel")
1036            kern.SetSharpness(2)
1037        elif kernel.lower() == "linear":
1038            kern = vtk.new("LinearKernel")
1039        else:
1040            vedo.logger.error("available kernels are: [shepard, gaussian, linear]")
1041            raise RuntimeError()
1042
1043        if n:
1044            kern.SetNumberOfPoints(n)
1045            kern.SetKernelFootprintToNClosest()
1046        else:
1047            kern.SetRadius(radius)
1048            kern.SetKernelFootprintToRadius()
1049
1050        interpolator = vtk.new("PointInterpolator")
1051        interpolator.SetInputData(self.dataset)
1052        interpolator.SetSourceData(points)
1053        interpolator.SetKernel(kern)
1054        interpolator.SetLocator(locator)
1055        interpolator.PassFieldArraysOff()
1056        interpolator.SetNullPointsStrategy(null_strategy)
1057        interpolator.SetNullValue(null_value)
1058        interpolator.SetValidPointsMaskArrayName("ValidPointMask")
1059        for ex in exclude:
1060            interpolator.AddExcludedArray(ex)
1061        interpolator.Update()
1062
1063        if on == "cells":
1064            p2c = vtk.new("PointDataToCellData")
1065            p2c.SetInputData(interpolator.GetOutput())
1066            p2c.Update()
1067            cpoly = p2c.GetOutput()
1068        else:
1069            cpoly = interpolator.GetOutput()
1070
1071        self._update(cpoly, reset_locators=False)
1072
1073        self.pipeline = utils.OperationNode("interpolate_data_from", parents=[self, source])
1074        return self
1075
1076    def add_ids(self):
1077        """
1078        Generate point and cell ids arrays.
1079        
1080        Two new arrays are added to the mesh: `PointID` and `CellID`.
1081        """
1082        ids = vtk.new("IdFilter")
1083        ids.SetInputData(self.dataset)
1084        ids.PointIdsOn()
1085        ids.CellIdsOn()
1086        ids.FieldDataOff()
1087        ids.SetPointIdsArrayName("PointID")
1088        ids.SetCellIdsArrayName("CellID")
1089        ids.Update()
1090        self._update(ids.GetOutput(), reset_locators=False)
1091        self.pipeline = utils.OperationNode("add_ids", parents=[self])
1092        return self
1093
1094    def gradient(self, input_array=None, on="points", fast=False):
1095        """
1096        Compute and return the gradiend of the active scalar field as a numpy array.
1097
1098        Arguments:
1099            input_array : (str)
1100                array of the scalars to compute the gradient,
1101                if None the current active array is selected
1102            on : (str)
1103                compute either on 'points' or 'cells' data
1104            fast : (bool)
1105                if True, will use a less accurate algorithm
1106                that performs fewer derivative calculations (and is therefore faster).
1107
1108        Examples:
1109            - [isolines.py](https://github.com/marcomusy/vedo/tree/master/examples/advanced/isolines.py)
1110
1111            ![](https://user-images.githubusercontent.com/32848391/72433087-f00a8780-3798-11ea-9778-991f0abeca70.png)
1112        """
1113        gra = vtk.new("GradientFilter")
1114        if on.startswith("p"):
1115            varr = self.dataset.GetPointData()
1116            tp = vtk.vtkDataObject.FIELD_ASSOCIATION_POINTS
1117        elif on.startswith("c"):
1118            varr = self.dataset.GetCellData()
1119            tp = vtk.vtkDataObject.FIELD_ASSOCIATION_CELLS
1120        else:
1121            vedo.logger.error(f"in gradient: unknown option {on}")
1122            raise RuntimeError
1123
1124        if input_array is None:
1125            if varr.GetScalars():
1126                input_array = varr.GetScalars().GetName()
1127            else:
1128                vedo.logger.error(f"in gradient: no scalars found for {on}")
1129                raise RuntimeError
1130
1131        gra.SetInputData(self.dataset)
1132        gra.SetInputScalars(tp, input_array)
1133        gra.SetResultArrayName("Gradient")
1134        gra.SetFasterApproximation(fast)
1135        gra.ComputeDivergenceOff()
1136        gra.ComputeVorticityOff()
1137        gra.ComputeGradientOn()
1138        gra.Update()
1139        if on.startswith("p"):
1140            gvecs = utils.vtk2numpy(gra.GetOutput().GetPointData().GetArray("Gradient"))
1141        else:
1142            gvecs = utils.vtk2numpy(gra.GetOutput().GetCellData().GetArray("Gradient"))
1143        return gvecs
1144
1145    def divergence(self, array_name=None, on="points", fast=False):
1146        """
1147        Compute and return the divergence of a vector field as a numpy array.
1148
1149        Arguments:
1150            array_name : (str)
1151                name of the array of vectors to compute the divergence,
1152                if None the current active array is selected
1153            on : (str)
1154                compute either on 'points' or 'cells' data
1155            fast : (bool)
1156                if True, will use a less accurate algorithm
1157                that performs fewer derivative calculations (and is therefore faster).
1158        """
1159        div = vtk.new("GradientFilter")
1160        if on.startswith("p"):
1161            varr = self.dataset.GetPointData()
1162            tp = vtk.vtkDataObject.FIELD_ASSOCIATION_POINTS
1163        elif on.startswith("c"):
1164            varr = self.dataset.GetCellData()
1165            tp = vtk.vtkDataObject.FIELD_ASSOCIATION_CELLS
1166        else:
1167            vedo.logger.error(f"in divergence(): unknown option {on}")
1168            raise RuntimeError
1169
1170        if array_name is None:
1171            if varr.GetVectors():
1172                array_name = varr.GetVectors().GetName()
1173            else:
1174                vedo.logger.error(f"in divergence(): no vectors found for {on}")
1175                raise RuntimeError
1176
1177        div.SetInputData(self.dataset)
1178        div.SetInputScalars(tp, array_name)
1179        div.ComputeDivergenceOn()
1180        div.ComputeGradientOff()
1181        div.ComputeVorticityOff()
1182        div.SetDivergenceArrayName("Divergence")
1183        div.SetFasterApproximation(fast)
1184        div.Update()
1185        if on.startswith("p"):
1186            dvecs = utils.vtk2numpy(div.GetOutput().GetPointData().GetArray("Divergence"))
1187        else:
1188            dvecs = utils.vtk2numpy(div.GetOutput().GetCellData().GetArray("Divergence"))
1189        return dvecs
1190
1191    def vorticity(self, array_name=None, on="points", fast=False):
1192        """
1193        Compute and return the vorticity of a vector field as a numpy array.
1194
1195        Arguments:
1196            array_name : (str)
1197                name of the array to compute the vorticity,
1198                if None the current active array is selected
1199            on : (str)
1200                compute either on 'points' or 'cells' data
1201            fast : (bool)
1202                if True, will use a less accurate algorithm
1203                that performs fewer derivative calculations (and is therefore faster).
1204        """
1205        vort = vtk.new("GradientFilter")
1206        if on.startswith("p"):
1207            varr = self.dataset.GetPointData()
1208            tp = vtk.vtkDataObject.FIELD_ASSOCIATION_POINTS
1209        elif on.startswith("c"):
1210            varr = self.dataset.GetCellData()
1211            tp = vtk.vtkDataObject.FIELD_ASSOCIATION_CELLS
1212        else:
1213            vedo.logger.error(f"in vorticity(): unknown option {on}")
1214            raise RuntimeError
1215
1216        if array_name is None:
1217            if varr.GetVectors():
1218                array_name = varr.GetVectors().GetName()
1219            else:
1220                vedo.logger.error(f"in vorticity(): no vectors found for {on}")
1221                raise RuntimeError
1222
1223        vort.SetInputData(self.dataset)
1224        vort.SetInputScalars(tp, array_name)
1225        vort.ComputeDivergenceOff()
1226        vort.ComputeGradientOff()
1227        vort.ComputeVorticityOn()
1228        vort.SetVorticityArrayName("Vorticity")
1229        vort.SetFasterApproximation(fast)
1230        vort.Update()
1231        if on.startswith("p"):
1232            vvecs = utils.vtk2numpy(vort.GetOutput().GetPointData().GetArray("Vorticity"))
1233        else:
1234            vvecs = utils.vtk2numpy(vort.GetOutput().GetCellData().GetArray("Vorticity"))
1235        return vvecs
1236
1237    def probe(self, source):
1238        """
1239        Takes a `Volume` (or any other data set)
1240        and probes its scalars at the specified points in space.
1241
1242        Note that a mask is also output with valid/invalid points which can be accessed
1243        with `mesh.pointdata['ValidPointMask']`.
1244
1245        Check out also:
1246            `interpolate_data_from()`
1247
1248        Examples:
1249            - [probe_points.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/probe_points.py)
1250
1251                ![](https://vedo.embl.es/images/volumetric/probePoints.png)
1252        """
1253        probe_filter = vtk.new("ProbeFilter")
1254        probe_filter.SetSourceData(source.dataset)
1255        probe_filter.SetInputData(self.dataset)
1256        probe_filter.Update()
1257        self._update(probe_filter.GetOutput(), reset_locators=False)
1258        self.pipeline = utils.OperationNode("probe", parents=[self, source])
1259        self.pointdata.rename("vtkValidPointMask", "ValidPointMask")
1260        return self
1261
1262    def compute_cell_size(self):
1263        """
1264        Add to this object a cell data array 
1265        containing the area, volume and edge length
1266        of the cells (when applicable to the object type).
1267
1268        Array names are: `Area`, `Volume`, `Length`.
1269        """
1270        csf = vtk.new("CellSizeFilter")
1271        csf.SetInputData(self.dataset)
1272        csf.SetComputeArea(1)
1273        csf.SetComputeVolume(1)
1274        csf.SetComputeLength(1)
1275        csf.SetComputeVertexCount(0)
1276        csf.SetAreaArrayName("Area")
1277        csf.SetVolumeArrayName("Volume")
1278        csf.SetLengthArrayName("Length")
1279        csf.Update()
1280        self._update(csf.GetOutput(), reset_locators=False)
1281        return self
1282    
1283    def integrate_arrays_over_domain(self):
1284        """
1285        Integrate point and cell data arrays while computing length, area or volume
1286        of the domain. It works for 1D, 2D or 3D cells.
1287        For volumetric datasets, this filter ignores all but 3D cells.
1288        It will not compute the volume contained in a closed surface. 
1289        
1290        Returns a dictionary with keys: `pointdata`, `celldata`, `metadata`,
1291        which contain the integration result for the corresponding attributes.
1292
1293        Examples:
1294            ```python
1295            from vedo import *
1296            surf = Sphere(res=100)
1297            surf.pointdata['scalars'] = np.ones(surf.npoints)
1298            data = surf.integrate_arrays_over_domain()
1299            print(data['pointdata']['scalars'], "is equal to 4pi", 4*np.pi)
1300            ```
1301
1302            ```python
1303            from vedo import *
1304
1305            xcoords1 = np.arange(0, 2.2, 0.2)
1306            xcoords2 = sqrt(np.arange(0, 4.2, 0.2))
1307
1308            ycoords = np.arange(0, 1.2, 0.2)
1309
1310            surf1 = Grid(s=(xcoords1, ycoords)).rotate_y(-45).lw(2)
1311            surf2 = Grid(s=(xcoords2, ycoords)).rotate_y(-45).lw(2)
1312
1313            surf1.pointdata['scalars'] = surf1.vertices[:,2]
1314            surf2.pointdata['scalars'] = surf2.vertices[:,2]
1315
1316            data1 = surf1.integrate_arrays_over_domain()
1317            data2 = surf2.integrate_arrays_over_domain()
1318
1319            print(data1['pointdata']['scalars'],
1320                "is equal to", 
1321                data2['pointdata']['scalars'],
1322                "even if the grids are different!",
1323                "(= the volume under the surface)"
1324            )
1325            show(surf1, surf2, N=2, axes=1).close()
1326        ```
1327        """
1328        vinteg = vtk.new("IntegrateAttributes")
1329        vinteg.SetInputData(self.dataset)
1330        vinteg.Update()
1331        ugrid = vedo.UnstructuredGrid(vinteg.GetOutput())
1332        data = dict(
1333            pointdata=ugrid.pointdata.todict(), 
1334            celldata=ugrid.celldata.todict(), 
1335            metadata=ugrid.metadata.todict(),
1336        )
1337        return data
1338
1339    def write(self, filename, binary=True):
1340        """Write object to file."""
1341        out = vedo.file_io.write(self, filename, binary)
1342        out.pipeline = utils.OperationNode(
1343            "write", parents=[self], comment=filename[:15], shape="folder", c="#8a817c"
1344        )
1345        return out
1346
1347    def tomesh(self, bounds=()):
1348        """
1349        Extract boundary geometry from dataset (or convert data to polygonal type).
1350        
1351        Two new arrays are added to the mesh: `OriginalCellIds` and `OriginalPointIds`
1352        to keep track of the original mesh elements.
1353        """
1354        geo = vtk.new("GeometryFilter")
1355        geo.SetInputData(self.dataset)
1356        geo.SetPassThroughCellIds(1)
1357        geo.SetPassThroughPointIds(1)
1358        geo.SetOriginalCellIdsName("OriginalCellIds")
1359        geo.SetOriginalPointIdsName("OriginalPointIds")
1360        geo.SetNonlinearSubdivisionLevel(1)
1361        geo.MergingOff()
1362        if bounds:
1363            geo.SetExtent(bounds)
1364            geo.ExtentClippingOn()
1365        geo.Update()
1366        msh = vedo.mesh.Mesh(geo.GetOutput())
1367        msh.pipeline = utils.OperationNode("tomesh", parents=[self], c="#9e2a2b")
1368        return msh
1369
1370    def shrink(self, fraction=0.8):
1371        """
1372        Shrink the individual cells to improve visibility.
1373
1374        ![](https://vedo.embl.es/images/feats/shrink_hex.png)
1375        """
1376        sf = vtk.new("ShrinkFilter")
1377        sf.SetInputData(self.dataset)
1378        sf.SetShrinkFactor(fraction)
1379        sf.Update()
1380        self._update(sf.GetOutput())
1381        self.pipeline = utils.OperationNode(
1382            "shrink", comment=f"by {fraction}", parents=[self], c="#9e2a2b"
1383        )
1384        return self
1385
1386
1387###############################################################################
1388class PointAlgorithms(CommonAlgorithms):
1389    """Methods for point clouds."""
1390
1391    def __init__(self):
1392        # print('init PointAlgorithms')
1393        super().__init__()
1394
1395        self.transform = None
1396        self.point_locator = None
1397        self.cell_locator = None
1398        self.line_locator = None
1399
1400    def apply_transform(self, LT, concatenate=True, deep_copy=True):
1401        """
1402        Apply a linear or non-linear transformation to the mesh polygonal data.
1403
1404        Example:
1405        ```python
1406        from vedo import Cube, show, settings
1407        settings.use_parallel_projection = True
1408        c1 = Cube().rotate_z(25).pos(2,1).mirror().alpha(0.5)
1409        T = c1.transform  # rotate by 5 degrees, place at (2,1)
1410        c2 = Cube().c('red4').wireframe().lw(10).lighting('off')
1411        c2.apply_transform(T)
1412        show(c1, c2, "The 2 cubes should overlap!", axes=1).close()
1413        ```
1414
1415        ![](https://vedo.embl.es/images/feats/apply_transform.png)
1416        """
1417        if self.dataset.GetNumberOfPoints() == 0:
1418            return self
1419        
1420        if isinstance(LT, LinearTransform):
1421            LT_is_linear = True
1422            tr = LT.T
1423            if LT.is_identity():
1424                return self
1425        
1426        elif isinstance(LT, (vtk.vtkMatrix4x4, vtk.vtkLinearTransform)) or utils.is_sequence(LT):
1427            LT_is_linear = True
1428            LT = LinearTransform(LT)
1429            tr = LT.T
1430            if LT.is_identity():
1431                return self
1432        
1433        elif isinstance(LT, NonLinearTransform):
1434            LT_is_linear = False
1435            tr = LT.T
1436            self.transform = LT # reset
1437
1438        elif isinstance(LT, vtk.vtkThinPlateSplineTransform):
1439            LT_is_linear = False
1440            tr = LT
1441            self.transform = NonLinearTransform(LT) # reset
1442        
1443        else:
1444            vedo.logger.error(f"apply_transform(), unknown input type:\n{LT}")
1445            return self
1446        
1447        ################
1448        if LT_is_linear:
1449            if concatenate:
1450                try:
1451                    # self.transform might still not be linear
1452                    self.transform.concatenate(LT)
1453                except AttributeError:
1454                    # in that case reset it
1455                    self.transform = LinearTransform()
1456
1457        ################
1458        if isinstance(self.dataset, vtk.vtkPolyData):
1459            tp = vtk.new("TransformPolyDataFilter")
1460        elif isinstance(self.dataset, vtk.vtkUnstructuredGrid):
1461            tp = vtk.new("TransformFilter")
1462            tp.TransformAllInputVectorsOn()
1463        # elif isinstance(self.dataset, vtk.vtkImageData):
1464        #     tp = vtk.new("ImageReslice")
1465        #     tp.SetInterpolationModeToCubic()
1466        #     tp.SetResliceTransform(tr)
1467        else:
1468            vedo.logger.error(
1469                f"apply_transform(), unknown input type: {[self.dataset]}")
1470            return self
1471        tp.SetTransform(tr)
1472        tp.SetInputData(self.dataset)
1473        tp.Update()
1474        out = tp.GetOutput()
1475
1476        if deep_copy:
1477            self.dataset.DeepCopy(out)
1478        else:
1479            self.dataset.ShallowCopy(out)
1480
1481        # reset the locators
1482        self.point_locator = None
1483        self.cell_locator = None
1484        self.line_locator = None
1485        return self
1486
1487    def apply_transform_from_actor(self):
1488        """
1489        Apply the current transformation of the actor to the data.
1490        Useful when manually moving an actor (eg. when pressing "a").
1491        Returns the `LinearTransform` object.
1492        """
1493        M = self.actor.GetMatrix()
1494        self.apply_transform(M)
1495        iden = vtk.vtkMatrix4x4()
1496        self.actor.PokeMatrix(iden)
1497        return LinearTransform(M)
1498
1499    def pos(self, x=None, y=None, z=None):
1500        """Set/Get object position."""
1501        if x is None:  # get functionality
1502            return self.transform.position
1503
1504        if z is None and y is None:  # assume x is of the form (x,y,z)
1505            if len(x) == 3:
1506                x, y, z = x
1507            else:
1508                x, y = x
1509                z = 0
1510        elif z is None:  # assume x,y is of the form x, y
1511            z = 0
1512
1513        q = self.transform.position
1514        delta = [x, y, z] - q
1515        if delta[0] == delta[1] == delta[2] == 0:
1516            return self
1517        LT = LinearTransform().translate(delta)
1518        return self.apply_transform(LT)
1519
1520    def shift(self, dx=0, dy=0, dz=0):
1521        """Add a vector to the current object position."""
1522        if utils.is_sequence(dx):
1523            utils.make3d(dx)
1524            dx, dy, dz = dx
1525        if dx == dy == dz == 0:
1526            return self
1527        LT = LinearTransform().translate([dx, dy, dz])
1528        return self.apply_transform(LT)
1529
1530    def x(self, val=None):
1531        """Set/Get object position along x axis."""
1532        p = self.transform.position
1533        if val is None:
1534            return p[0]
1535        self.pos(val, p[1], p[2])
1536        return self
1537
1538    def y(self, val=None):
1539        """Set/Get object position along y axis."""
1540        p = self.transform.position
1541        if val is None:
1542            return p[1]
1543        self.pos(p[0], val, p[2])
1544        return self
1545
1546    def z(self, val=None):
1547        """Set/Get object position along z axis."""
1548        p = self.transform.position
1549        if val is None:
1550            return p[2]
1551        self.pos(p[0], p[1], val)
1552        return self
1553
1554    def rotate(self, angle, axis=(1, 0, 0), point=(0, 0, 0), rad=False):
1555        """
1556        Rotate around an arbitrary `axis` passing through `point`.
1557
1558        Example:
1559        ```python
1560        from vedo import *
1561        c1 = Cube()
1562        c2 = c1.clone().c('violet').alpha(0.5) # copy of c1
1563        v = vector(0.2,1,0)
1564        p = vector(1,0,0)  # axis passes through this point
1565        c2.rotate(90, axis=v, point=p)
1566        l = Line(-v+p, v+p).lw(3).c('red')
1567        show(c1, l, c2, axes=1).close()
1568        ```
1569
1570        ![](https://vedo.embl.es/images/feats/rotate_axis.png)
1571        """
1572        LT = LinearTransform()
1573        LT.rotate(angle, axis, point, rad)
1574        return self.apply_transform(LT)
1575
1576    def rotate_x(self, angle, rad=False, around=None):
1577        """
1578        Rotate around x-axis. If angle is in radians set `rad=True`.
1579
1580        Use `around` to define a pivoting point.
1581        """
1582        if angle == 0:
1583            return self
1584        LT = LinearTransform().rotate_x(angle, rad, around)
1585        return self.apply_transform(LT)
1586
1587    def rotate_y(self, angle, rad=False, around=None):
1588        """
1589        Rotate around y-axis. If angle is in radians set `rad=True`.
1590
1591        Use `around` to define a pivoting point.
1592        """
1593        if angle == 0:
1594            return self
1595        LT = LinearTransform().rotate_y(angle, rad, around)
1596        return self.apply_transform(LT)
1597
1598    def rotate_z(self, angle, rad=False, around=None):
1599        """
1600        Rotate around z-axis. If angle is in radians set `rad=True`.
1601
1602        Use `around` to define a pivoting point.
1603        """
1604        if angle == 0:
1605            return self
1606        LT = LinearTransform().rotate_z(angle, rad, around)
1607        return self.apply_transform(LT)
1608
1609    def reorient(self, initaxis, newaxis, rotation=0, rad=False, xyplane=False):
1610        """
1611        Reorient the object to point to a new direction from an initial one.
1612        If `initaxis` is None, the object will be assumed in its "default" orientation.
1613        If `xyplane` is True, the object will be rotated to lie on the xy plane.
1614
1615        Use `rotation` to first rotate the object around its `initaxis`.
1616        """
1617        q = self.transform.position
1618        LT = LinearTransform()
1619        LT.reorient(initaxis, newaxis, q, rotation, rad, xyplane)
1620        return self.apply_transform(LT)
1621
1622    def scale(self, s=None, reset=False, origin=True):
1623        """
1624        Set/get object's scaling factor.
1625
1626        Arguments:
1627            s : (list, float)
1628                scaling factor(s).
1629            reset : (bool)
1630                if True previous scaling factors are ignored.
1631            origin : (bool)
1632                if True scaling is applied with respect to object's position,
1633                otherwise is applied respect to (0,0,0).
1634
1635        Note:
1636            use `s=(sx,sy,sz)` to scale differently in the three coordinates.
1637        """
1638        if s is None:
1639            return np.array(self.transform.T.GetScale())
1640
1641        if not utils.is_sequence(s):
1642            s = [s, s, s]
1643
1644        LT = LinearTransform()
1645        if reset:
1646            old_s = np.array(self.transform.T.GetScale())
1647            LT.scale(s / old_s)
1648        else:
1649            if origin is True:
1650                LT.scale(s, origin=self.transform.position)
1651            elif origin is False:
1652                LT.scale(s, origin=False)
1653            else:
1654                LT.scale(s, origin=origin)
1655
1656        return self.apply_transform(LT)
1657
1658
1659###############################################################################
1660class VolumeAlgorithms(CommonAlgorithms):
1661    """Methods for Volume objects."""
1662
1663    def __init__(self):
1664       super().__init__()
1665
1666    def bounds(self):
1667        """
1668        Get the object bounds.
1669        Returns a list in format `[xmin,xmax, ymin,ymax, zmin,zmax]`.
1670        """
1671        # OVERRIDE CommonAlgorithms.bounds() which is too slow
1672        return np.array(self.dataset.GetBounds())
1673
1674    def isosurface(self, value=None, flying_edges=True):
1675        """
1676        Return an `Mesh` isosurface extracted from the `Volume` object.
1677
1678        Set `value` as single float or list of values to draw the isosurface(s).
1679        Use flying_edges for faster results (but sometimes can interfere with `smooth()`).
1680
1681        Examples:
1682            - [isosurfaces.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/isosurfaces.py)
1683
1684                ![](https://vedo.embl.es/images/volumetric/isosurfaces.png)
1685        """
1686        scrange = self.dataset.GetScalarRange()
1687
1688        if flying_edges:
1689            cf = vtk.new("FlyingEdges3D")
1690            cf.InterpolateAttributesOn()
1691        else:
1692            cf = vtk.new("ContourFilter")
1693            cf.UseScalarTreeOn()
1694
1695        cf.SetInputData(self.dataset)
1696        cf.ComputeNormalsOn()
1697
1698        if utils.is_sequence(value):
1699            cf.SetNumberOfContours(len(value))
1700            for i, t in enumerate(value):
1701                cf.SetValue(i, t)
1702        else:
1703            if value is None:
1704                value = (2 * scrange[0] + scrange[1]) / 3.0
1705                # print("automatic isosurface value =", value)
1706            cf.SetValue(0, value)
1707
1708        cf.Update()
1709        poly = cf.GetOutput()
1710
1711        out = vedo.mesh.Mesh(poly, c=None).phong()
1712        out.mapper.SetScalarRange(scrange[0], scrange[1])
1713
1714        out.pipeline = utils.OperationNode(
1715            "isosurface",
1716            parents=[self],
1717            comment=f"#pts {out.dataset.GetNumberOfPoints()}",
1718            c="#4cc9f0:#e9c46a",
1719        )
1720        return out
1721
1722    def legosurface(
1723        self,
1724        vmin=None,
1725        vmax=None,
1726        invert=False,
1727        boundary=False,
1728        array_name="input_scalars",
1729    ):
1730        """
1731        Represent an object - typically a `Volume` - as lego blocks (voxels).
1732        By default colors correspond to the volume's scalar.
1733        Returns an `Mesh` object.
1734
1735        Arguments:
1736            vmin : (float)
1737                the lower threshold, voxels below this value are not shown.
1738            vmax : (float)
1739                the upper threshold, voxels above this value are not shown.
1740            boundary : (bool)
1741                controls whether to include cells that are partially inside
1742            array_name : (int, str)
1743                name or index of the scalar array to be considered
1744
1745        Examples:
1746            - [legosurface.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/legosurface.py)
1747
1748                ![](https://vedo.embl.es/images/volumetric/56820682-da40e500-684c-11e9-8ea3-91cbcba24b3a.png)
1749        """
1750        imp_dataset = vtk.new("ImplicitDataSet")
1751        imp_dataset.SetDataSet(self.dataset)
1752        window = vtk.new("ImplicitWindowFunction")
1753        window.SetImplicitFunction(imp_dataset)
1754
1755        srng = list(self.dataset.GetScalarRange())
1756        if vmin is not None:
1757            srng[0] = vmin
1758        if vmax is not None:
1759            srng[1] = vmax
1760        tol = 0.00001 * (srng[1] - srng[0])
1761        srng[0] -= tol
1762        srng[1] += tol
1763        window.SetWindowRange(srng)
1764
1765        extract = vtk.new("ExtractGeometry")
1766        extract.SetInputData(self.dataset)
1767        extract.SetImplicitFunction(window)
1768        extract.SetExtractInside(invert)
1769        extract.SetExtractBoundaryCells(boundary)
1770        extract.Update()
1771
1772        gf = vtk.new("GeometryFilter")
1773        gf.SetInputData(extract.GetOutput())
1774        gf.Update()
1775
1776        m = vedo.mesh.Mesh(gf.GetOutput()).lw(0.1).flat()
1777        m.map_points_to_cells()
1778        m.celldata.select(array_name)
1779
1780        m.pipeline = utils.OperationNode(
1781            "legosurface",
1782            parents=[self],
1783            comment=f"array: {array_name}",
1784            c="#4cc9f0:#e9c46a",
1785        )
1786        return m
1787
1788    def tomesh(self, fill=True, shrink=1.0):
1789        """
1790        Build a polygonal Mesh from the current object.
1791
1792        If `fill=True`, the interior faces of all the cells are created.
1793        (setting a `shrink` value slightly smaller than the default 1.0
1794        can avoid flickering due to internal adjacent faces).
1795
1796        If `fill=False`, only the boundary faces will be generated.
1797        """
1798        gf = vtk.new("GeometryFilter")
1799        if fill:
1800            sf = vtk.new("ShrinkFilter")
1801            sf.SetInputData(self.dataset)
1802            sf.SetShrinkFactor(shrink)
1803            sf.Update()
1804            gf.SetInputData(sf.GetOutput())
1805            gf.Update()
1806            poly = gf.GetOutput()
1807            if shrink == 1.0:
1808                clean_poly = vtk.new("CleanPolyData")
1809                clean_poly.PointMergingOn()
1810                clean_poly.ConvertLinesToPointsOn()
1811                clean_poly.ConvertPolysToLinesOn()
1812                clean_poly.ConvertStripsToPolysOn()
1813                clean_poly.SetInputData(poly)
1814                clean_poly.Update()
1815                poly = clean_poly.GetOutput()
1816        else:
1817            gf.SetInputData(self.dataset)
1818            gf.Update()
1819            poly = gf.GetOutput()
1820
1821        msh = vedo.mesh.Mesh(poly).flat()
1822        msh.scalarbar = self.scalarbar
1823        lut = utils.ctf2lut(self)
1824        if lut:
1825            msh.mapper.SetLookupTable(lut)
1826
1827        msh.pipeline = utils.OperationNode(
1828            "tomesh", parents=[self], comment=f"fill={fill}", c="#9e2a2b:#e9c46a"
1829        )
1830        return msh