vedo.volume

Work with volumetric datasets (voxel data).

   1import glob
   2import os
   3import time
   4from weakref import ref as weak_ref_to
   5from typing import Union, List, Iterable
   6from typing_extensions import Self
   7
   8import numpy as np
   9
  10import vedo.vtkclasses as vtki
  11
  12import vedo
  13from vedo import transformations
  14from vedo import utils
  15from vedo.mesh import Mesh
  16from vedo.core import VolumeAlgorithms
  17from vedo.visual import VolumeVisual
  18
  19
  20__docformat__ = "google"
  21
  22__doc__ = """
  23Work with volumetric datasets (voxel data).
  24
  25![](https://vedo.embl.es/images/volumetric/slicePlane2.png)
  26"""
  27
  28__all__ = ["Volume"]
  29
  30
  31##########################################################################
  32class Volume(VolumeAlgorithms, VolumeVisual):
  33    """
  34    Class to describe dataset that are defined on "voxels",
  35    the 3D equivalent of 2D pixels.
  36    """
  37    def __init__(
  38        self,
  39        inputobj=None,
  40        dims=None,
  41        origin=None,
  42        spacing=None,
  43    ) -> None:
  44        """
  45        This class can be initialized with a numpy object,
  46        a `vtkImageData` or a list of 2D bmp files.
  47
  48        Arguments:
  49            origin : (list)
  50                set volume origin coordinates
  51            spacing : (list)
  52                voxel dimensions in x, y and z.
  53            dims : (list)
  54                specify the dimensions of the volume.
  55
  56        Example:
  57            ```python
  58            from vedo import Volume
  59            vol = Volume("path/to/mydata/rec*.bmp")
  60            vol.show()
  61            ```
  62
  63        Examples:
  64            - [numpy2volume1.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/numpy2volume1.py)
  65
  66                ![](https://vedo.embl.es/images/volumetric/numpy2volume1.png)
  67
  68            - [read_volume2.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/read_volume2.py)
  69
  70                ![](https://vedo.embl.es/images/volumetric/read_volume2.png)
  71
  72        .. note::
  73            if a `list` of values is used for `alphas` this is interpreted
  74            as a transfer function along the range of the scalar.
  75        """
  76        super().__init__()
  77
  78        self.name = "Volume"
  79        self.filename = ""
  80        self.file_size = ""
  81
  82        self.info = {}
  83        self.time =  time.time()
  84
  85        self.actor = vtki.vtkVolume()
  86        self.actor.retrieve_object = weak_ref_to(self)
  87        self.properties = self.actor.GetProperty()
  88
  89        self.transform = None
  90        self.point_locator = None
  91        self.cell_locator = None
  92        self.line_locator = None
  93
  94        ###################
  95        if isinstance(inputobj, str):
  96            if "https://" in inputobj:
  97                inputobj = vedo.file_io.download(inputobj, verbose=False)  # fpath
  98            elif os.path.isfile(inputobj):
  99                self.filename = inputobj
 100            else:
 101                inputobj = sorted(glob.glob(inputobj))
 102
 103        ###################
 104        inputtype = str(type(inputobj))
 105
 106        # print('Volume inputtype', inputtype, c='b')
 107
 108        if inputobj is None:
 109            img = vtki.vtkImageData()
 110
 111        elif utils.is_sequence(inputobj):
 112
 113            if isinstance(inputobj[0], str) and ".bmp" in inputobj[0].lower():
 114                # scan sequence of BMP files
 115                ima = vtki.new("ImageAppend")
 116                ima.SetAppendAxis(2)
 117                pb = utils.ProgressBar(0, len(inputobj))
 118                for i in pb.range():
 119                    f = inputobj[i]
 120                    if "_rec_spr" in f: # OPT specific
 121                        continue
 122                    picr = vtki.new("BMPReader")
 123                    picr.SetFileName(f)
 124                    picr.Update()
 125                    mgf = vtki.new("ImageMagnitude")
 126                    mgf.SetInputData(picr.GetOutput())
 127                    mgf.Update()
 128                    ima.AddInputData(mgf.GetOutput())
 129                    pb.print("loading...")
 130                ima.Update()
 131                img = ima.GetOutput()
 132
 133            else:
 134
 135                if len(inputobj.shape) == 1:
 136                    varr = utils.numpy2vtk(inputobj)
 137                else:
 138                    varr = utils.numpy2vtk(inputobj.ravel(order="F"))
 139                varr.SetName("input_scalars")
 140
 141                img = vtki.vtkImageData()
 142                if dims is not None:
 143                    img.SetDimensions(dims[2], dims[1], dims[0])
 144                else:
 145                    if len(inputobj.shape) == 1:
 146                        vedo.logger.error("must set dimensions (dims keyword) in Volume")
 147                        raise RuntimeError()
 148                    img.SetDimensions(inputobj.shape)
 149                img.GetPointData().AddArray(varr)
 150                img.GetPointData().SetActiveScalars(varr.GetName())
 151
 152        elif isinstance(inputobj, vtki.vtkImageData):
 153            img = inputobj
 154
 155        elif isinstance(inputobj, str):
 156            if "https://" in inputobj:
 157                inputobj = vedo.file_io.download(inputobj, verbose=False)
 158            img = vedo.file_io.loadImageData(inputobj)
 159            self.filename = inputobj
 160
 161        else:
 162            vedo.logger.error(f"cannot understand input type {inputtype}")
 163            return
 164
 165        if dims is not None:
 166            img.SetDimensions(dims)
 167
 168        if origin is not None:
 169            img.SetOrigin(origin)
 170
 171        if spacing is not None:
 172            img.SetSpacing(spacing)
 173
 174        self.dataset = img
 175        self.transform = None
 176
 177        #####################################
 178        mapper = vtki.new("SmartVolumeMapper")
 179        mapper.SetInputData(img)
 180        self.actor.SetMapper(mapper)
 181
 182        if img.GetPointData().GetScalars():
 183            if img.GetPointData().GetScalars().GetNumberOfComponents() == 1:
 184                self.properties.SetShade(True)
 185                self.properties.SetInterpolationType(1)
 186                self.cmap("RdBu_r")
 187                self.alpha([0.0, 0.0, 0.2, 0.4, 0.8, 1.0])
 188                self.alpha_gradient(None)
 189                self.properties.SetScalarOpacityUnitDistance(1.0)
 190
 191        self.pipeline = utils.OperationNode(
 192            "Volume", comment=f"dims={tuple(self.dimensions())}", c="#4cc9f0"
 193        )
 194        #######################################################################
 195
 196    @property
 197    def mapper(self):
 198        """Return the underlying `vtkMapper` object."""
 199        return self.actor.GetMapper()
 200    
 201    @mapper.setter
 202    def mapper(self, mapper):
 203        """
 204        Set the underlying `vtkMapper` object.
 205        
 206        Arguments:
 207            mapper : (str, vtkMapper)
 208                either 'gpu', 'opengl_gpu', 'fixed' or 'smart'
 209        """
 210        if isinstance(mapper, 
 211            (vtki.get_class("Mapper"), vtki.get_class("ImageResliceMapper"))
 212        ):
 213            pass
 214        elif mapper is None:
 215            mapper = vtki.new("SmartVolumeMapper")
 216        elif "gpu" in mapper:
 217            mapper = vtki.new("GPUVolumeRayCastMapper")
 218        elif "opengl_gpu" in mapper:
 219            mapper = vtki.new("OpenGLGPUVolumeRayCastMapper")
 220        elif "smart" in mapper:
 221            mapper = vtki.new("SmartVolumeMapper")
 222        elif "fixed" in mapper:
 223            mapper = vtki.new("FixedPointVolumeRayCastMapper")
 224        else:
 225            print("Error unknown mapper type", [mapper])
 226            raise RuntimeError()
 227        self.actor.SetMapper(mapper)
 228
 229    def c(self, *args, **kwargs) -> Self:
 230        """Deprecated. Use `Volume.cmap()` instead."""
 231        vedo.logger.warning("Volume.c() is deprecated, use Volume.cmap() instead")
 232        return self.cmap(*args, **kwargs)
 233
 234    def _update(self, data, reset_locators=False):
 235        # reset_locators here is dummy
 236        self.dataset = data
 237        self.mapper.SetInputData(data)
 238        self.dataset.GetPointData().Modified()
 239        self.mapper.Modified()
 240        self.mapper.Update()
 241        return self
 242
 243    def __str__(self):
 244        """Print a summary for the `Volume` object."""
 245        module = self.__class__.__module__
 246        name = self.__class__.__name__
 247        out = vedo.printc(
 248            f"{module}.{name} at ({hex(self.memory_address())})".ljust(75),
 249            c="c", bold=True, invert=True, return_string=True,
 250        )
 251        out += "\x1b[0m\x1b[36;1m"
 252
 253        out+= "name".ljust(14) + ": " + str(self.name) + "\n"
 254        if self.filename:
 255            out+= "filename".ljust(14) + ": " + str(self.filename) + "\n"
 256
 257        out+= "dimensions".ljust(14) + ": " + str(self.shape) + "\n"
 258
 259        out+= "origin".ljust(14) + ": "
 260        out+= utils.precision(self.origin(), 6) + "\n"
 261
 262        out+= "center".ljust(14) + ": "
 263        out+= utils.precision(self.center(), 6) + "\n"
 264
 265        out+= "spacing".ljust(14)    + ": "
 266        out+= utils.precision(self.spacing(), 6) + "\n"
 267
 268        bnds = self.bounds()
 269        bx1, bx2 = utils.precision(bnds[0], 3), utils.precision(bnds[1], 3)
 270        by1, by2 = utils.precision(bnds[2], 3), utils.precision(bnds[3], 3)
 271        bz1, bz2 = utils.precision(bnds[4], 3), utils.precision(bnds[5], 3)
 272        out+= "bounds".ljust(14) + ":"
 273        out+= " x=(" + bx1 + ", " + bx2 + "),"
 274        out+= " y=(" + by1 + ", " + by2 + "),"
 275        out+= " z=(" + bz1 + ", " + bz2 + ")\n"
 276
 277        out+= "memory size".ljust(14) + ": "
 278        out+= str(int(self.dataset.GetActualMemorySize()/1024+0.5))+" MB\n"
 279
 280        st = self.dataset.GetScalarTypeAsString()
 281        out+= "scalar size".ljust(14) + ": "
 282        out+= str(self.dataset.GetScalarSize()) + f" bytes ({st})\n"
 283        out+= "scalar range".ljust(14) + ": "
 284        out+= str(self.dataset.GetScalarRange()) + "\n"
 285
 286        #utils.print_histogram(self, logscale=True, bins=8, height=15, c="b", bold=True)
 287        return out.rstrip() + "\x1b[0m"
 288
 289    def _repr_html_(self):
 290        """
 291        HTML representation of the Volume object for Jupyter Notebooks.
 292
 293        Returns:
 294            HTML text with the image and some properties.
 295        """
 296        import io
 297        import base64
 298        from PIL import Image
 299
 300        library_name = "vedo.volume.Volume"
 301        help_url = "https://vedo.embl.es/docs/vedo/volume.html"
 302
 303        arr = self.thumbnail(azimuth=0, elevation=-60, zoom=1.4, axes=True)
 304
 305        im = Image.fromarray(arr)
 306        buffered = io.BytesIO()
 307        im.save(buffered, format="PNG", quality=100)
 308        encoded = base64.b64encode(buffered.getvalue()).decode("utf-8")
 309        url = "data:image/png;base64," + encoded
 310        image = f"<img src='{url}'></img>"
 311
 312        # statisitics
 313        bounds = "<br/>".join(
 314            [
 315                utils.precision(min_x, 4) + " ... " + utils.precision(max_x, 4)
 316                for min_x, max_x in zip(self.bounds()[::2], self.bounds()[1::2])
 317            ]
 318        )
 319
 320        help_text = ""
 321        if self.name:
 322            help_text += f"<b> {self.name}: &nbsp&nbsp</b>"
 323        help_text += '<b><a href="' + help_url + '" target="_blank">' + library_name + "</a></b>"
 324        if self.filename:
 325            dots = ""
 326            if len(self.filename) > 30:
 327                dots = "..."
 328            help_text += f"<br/><code><i>({dots}{self.filename[-30:]})</i></code>"
 329
 330        pdata = ""
 331        if self.dataset.GetPointData().GetScalars():
 332            if self.dataset.GetPointData().GetScalars().GetName():
 333                name = self.dataset.GetPointData().GetScalars().GetName()
 334                pdata = "<tr><td><b> point data array </b></td><td>" + name + "</td></tr>"
 335
 336        cdata = ""
 337        if self.dataset.GetCellData().GetScalars():
 338            if self.dataset.GetCellData().GetScalars().GetName():
 339                name = self.dataset.GetCellData().GetScalars().GetName()
 340                cdata = "<tr><td><b> voxel data array </b></td><td>" + name + "</td></tr>"
 341
 342        img = self.dataset
 343
 344        allt = [
 345            "<table>",
 346            "<tr>",
 347            "<td>",
 348            image,
 349            "</td>",
 350            "<td style='text-align: center; vertical-align: center;'><br/>",
 351            help_text,
 352            "<table>",
 353            "<tr><td><b> bounds </b> <br/> (x/y/z) </td><td>" + str(bounds) + "</td></tr>",
 354            "<tr><td><b> dimensions </b></td><td>" + str(img.GetDimensions()) + "</td></tr>",
 355            "<tr><td><b> voxel spacing </b></td><td>"
 356            + utils.precision(img.GetSpacing(), 3)
 357            + "</td></tr>",
 358            "<tr><td><b> in memory size </b></td><td>"
 359            + str(int(img.GetActualMemorySize() / 1024))
 360            + "MB</td></tr>",
 361            pdata,
 362            cdata,
 363            "<tr><td><b> scalar range </b></td><td>"
 364            + utils.precision(img.GetScalarRange(), 4)
 365            + "</td></tr>",
 366            "</table>",
 367            "</table>",
 368        ]
 369        return "\n".join(allt)
 370
 371    def copy(self, deep=True) -> "Volume":
 372        """Return a copy of the Volume. Alias of `clone()`."""
 373        return self.clone(deep=deep)
 374
 375    def clone(self, deep=True) -> "Volume":
 376        """Return a clone copy of the Volume. Alias of `copy()`."""
 377        if deep:
 378            newimg = vtki.vtkImageData()
 379            newimg.CopyStructure(self.dataset)
 380            newimg.CopyAttributes(self.dataset)
 381            newvol = Volume(newimg)
 382        else:
 383            newvol = Volume(self.dataset)
 384
 385        prop = vtki.vtkVolumeProperty()
 386        prop.DeepCopy(self.properties)
 387        newvol.actor.SetProperty(prop)
 388        newvol.properties = prop
 389
 390        newvol.pipeline = utils.OperationNode("clone", parents=[self], c="#bbd0ff", shape="diamond")
 391        return newvol
 392
 393    def astype(self, dtype: Union[str, int]) -> Self:
 394        """
 395        Reset the type of the scalars array.
 396
 397        Arguments:
 398            dtype : (str)
 399                the type of the scalars array in
 400                ["int8", "uint8", "int16", "uint16", "int32", "uint32", "float32", "float64"]
 401        """
 402        if dtype in ["int8", "uint8", "int16", "uint16", "int32", "uint32", "float32", "float64"]:
 403            caster = vtki.new("ImageCast")
 404            caster.SetInputData(self.dataset)
 405            caster.SetOutputScalarType(int(vtki.array_types[dtype]))
 406            caster.ClampOverflowOn()
 407            caster.Update()
 408            self._update(caster.GetOutput())
 409            self.pipeline = utils.OperationNode(f"astype({dtype})", parents=[self], c="#4cc9f0")
 410        else:
 411            vedo.logger.error(f"astype(): unknown type {dtype}")
 412            raise ValueError()
 413        return self
 414
 415    def component_weight(self, i: int, weight: float) -> Self:
 416        """Set the scalar component weight in range [0,1]."""
 417        self.properties.SetComponentWeight(i, weight)
 418        return self
 419
 420    def xslice(self, i: int) -> Mesh:
 421        """Extract the slice at index `i` of volume along x-axis."""
 422        vslice = vtki.new("ImageDataGeometryFilter")
 423        vslice.SetInputData(self.dataset)
 424        nx, ny, nz = self.dataset.GetDimensions()
 425        if i > nx - 1:
 426            i = nx - 1
 427        vslice.SetExtent(i, i, 0, ny, 0, nz)
 428        vslice.Update()
 429        m = Mesh(vslice.GetOutput())
 430        m.pipeline = utils.OperationNode(f"xslice {i}", parents=[self], c="#4cc9f0:#e9c46a")
 431        return m
 432
 433    def yslice(self, j: int) -> Mesh:
 434        """Extract the slice at index `j` of volume along y-axis."""
 435        vslice = vtki.new("ImageDataGeometryFilter")
 436        vslice.SetInputData(self.dataset)
 437        nx, ny, nz = self.dataset.GetDimensions()
 438        if j > ny - 1:
 439            j = ny - 1
 440        vslice.SetExtent(0, nx, j, j, 0, nz)
 441        vslice.Update()
 442        m = Mesh(vslice.GetOutput())
 443        m.pipeline = utils.OperationNode(f"yslice {j}", parents=[self], c="#4cc9f0:#e9c46a")
 444        return m
 445
 446    def zslice(self, k: int) -> Mesh:
 447        """Extract the slice at index `i` of volume along z-axis."""
 448        vslice = vtki.new("ImageDataGeometryFilter")
 449        vslice.SetInputData(self.dataset)
 450        nx, ny, nz = self.dataset.GetDimensions()
 451        if k > nz - 1:
 452            k = nz - 1
 453        vslice.SetExtent(0, nx, 0, ny, k, k)
 454        vslice.Update()
 455        m = Mesh(vslice.GetOutput())
 456        m.pipeline = utils.OperationNode(f"zslice {k}", parents=[self], c="#4cc9f0:#e9c46a")
 457        return m
 458
 459    def slice_plane(self, origin: List[float], normal: List[float], autocrop=False, border=0.5, mode="linear") -> Mesh:
 460        """
 461        Extract the slice along a given plane position and normal.
 462
 463        Two metadata arrays are added to the output Mesh:
 464            - "shape" : contains the shape of the slice
 465            - "original_bounds" : contains the original bounds of the slice
 466        One can access them with e.g. `myslice.metadata["shape"]`.
 467
 468        Arguments:
 469            origin : (list)
 470                position of the plane
 471            normal : (list)
 472                normal to the plane
 473            autocrop : (bool)
 474                crop the output to the minimal possible size
 475            border : (float)
 476                add a border to the output slice
 477            mode : (str)
 478                interpolation mode, one of the following: "linear", "nearest", "cubic"
 479
 480        Example:
 481            - [slice_plane1.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/slice_plane1.py)
 482
 483                ![](https://vedo.embl.es/images/volumetric/slicePlane1.gif)
 484            
 485            - [slice_plane2.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/slice_plane2.py)
 486                
 487                ![](https://vedo.embl.es/images/volumetric/slicePlane2.png)
 488
 489            - [slice_plane3.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/slice_plane3.py)
 490
 491                ![](https://vedo.embl.es/images/volumetric/slicePlane3.jpg)
 492        """
 493        newaxis = utils.versor(normal)
 494        pos = np.array(origin)
 495        initaxis = (0, 0, 1)
 496        crossvec = np.cross(initaxis, newaxis)
 497        angle = np.arccos(np.dot(initaxis, newaxis))
 498        T = vtki.vtkTransform()
 499        T.PostMultiply()
 500        T.RotateWXYZ(np.rad2deg(angle), crossvec.tolist())
 501        T.Translate(pos.tolist())
 502
 503        reslice = vtki.new("ImageReslice")
 504        reslice.SetResliceAxes(T.GetMatrix())
 505        reslice.SetInputData(self.dataset)
 506        reslice.SetOutputDimensionality(2)
 507        reslice.SetTransformInputSampling(True)
 508        reslice.SetGenerateStencilOutput(False)
 509        if border:
 510            reslice.SetBorder(True)
 511            reslice.SetBorderThickness(border)
 512        else:
 513            reslice.SetBorder(False)
 514        if mode == "linear":
 515            reslice.SetInterpolationModeToLinear()
 516        elif mode == "nearest":
 517            reslice.SetInterpolationModeToNearestNeighbor()
 518        elif mode == "cubic":
 519            reslice.SetInterpolationModeToCubic()
 520        else:
 521            vedo.logger.error(f"in slice_plane(): unknown interpolation mode {mode}")
 522            raise ValueError()
 523        reslice.SetAutoCropOutput(not autocrop)
 524        reslice.Update()
 525        img = reslice.GetOutput()
 526
 527        vslice = vtki.new("ImageDataGeometryFilter")
 528        vslice.SetInputData(img)
 529        vslice.Update()
 530
 531        msh = Mesh(vslice.GetOutput()).apply_transform(T)
 532        msh.properties.LightingOff()
 533
 534        d0, d1, _ = img.GetDimensions()
 535        varr1 = utils.numpy2vtk([d1, d0], name="shape")
 536        varr2 = utils.numpy2vtk(img.GetBounds(), name="original_bounds")
 537        msh.dataset.GetFieldData().AddArray(varr1)
 538        msh.dataset.GetFieldData().AddArray(varr2)
 539        msh.pipeline = utils.OperationNode("slice_plane", parents=[self], c="#4cc9f0:#e9c46a")
 540        return msh
 541    
 542    def slab(self, slice_range=(), axis='z', operation="mean") -> Mesh:
 543        """
 544        Extract a slab from a `Volume` by combining 
 545        all of the slices of an image to create a single slice.
 546
 547        Returns a `Mesh` containing metadata which
 548        can be accessed with e.g. `mesh.metadata["slab_range"]`.
 549
 550        Metadata:
 551            slab_range : (list)
 552                contains the range of slices extracted
 553            slab_axis : (str)
 554                contains the axis along which the slab was extracted
 555            slab_operation : (str)
 556                contains the operation performed on the slab
 557            slab_bounding_box : (list)
 558                contains the bounding box of the slab
 559
 560        Arguments:
 561            slice_range : (list)
 562                range of slices to extract
 563            axis : (str)
 564                axis along which to extract the slab
 565            operation : (str)
 566                operation to perform on the slab,
 567                allowed values are: "sum", "min", "max", "mean".
 568        
 569        Example:
 570            - [slab.py](https://github.com/marcomusy/vedo/blob/master/examples/volumetric/slab_vol.py)
 571
 572            ![](https://vedo.embl.es/images/volumetric/slab_vol.jpg)
 573        """
 574        if len(slice_range) != 2:
 575            vedo.logger.error("in slab(): slice_range is empty or invalid")
 576            raise ValueError()
 577        
 578        islab = vtki.new("ImageSlab")
 579        islab.SetInputData(self.dataset)
 580
 581        if operation in ["+", "add", "sum"]:
 582            islab.SetOperationToSum()
 583        elif "min" in operation:
 584            islab.SetOperationToMin()
 585        elif "max" in operation:
 586            islab.SetOperationToMax()
 587        elif "mean" in operation:
 588            islab.SetOperationToMean()
 589        else:
 590            vedo.logger.error(f"in slab(): unknown operation {operation}")
 591            raise ValueError()
 592
 593        dims = self.dimensions()
 594        if axis == 'x':
 595            islab.SetOrientationToX()
 596            if slice_range[0]  > dims[0]-1:
 597                slice_range[0] = dims[0]-1
 598            if slice_range[1]  > dims[0]-1:
 599                slice_range[1] = dims[0]-1
 600        elif axis == 'y':
 601            islab.SetOrientationToY()
 602            if slice_range[0]  > dims[1]-1:
 603                slice_range[0] = dims[1]-1
 604            if slice_range[1]  > dims[1]-1:
 605                slice_range[1] = dims[1]-1
 606        elif axis == 'z':
 607            islab.SetOrientationToZ()
 608            if slice_range[0]  > dims[2]-1:
 609                slice_range[0] = dims[2]-1
 610            if slice_range[1]  > dims[2]-1:
 611                slice_range[1] = dims[2]-1
 612        else:
 613            vedo.logger.error(f"Error in slab(): unknown axis {axis}")
 614            raise RuntimeError()
 615        
 616        islab.SetSliceRange(slice_range)
 617        islab.Update()
 618
 619        msh = Mesh(islab.GetOutput()).lighting('off')
 620        msh.mapper.SetLookupTable(utils.ctf2lut(self, msh))
 621        msh.mapper.SetScalarRange(self.scalar_range())
 622
 623        msh.metadata["slab_range"] = slice_range
 624        msh.metadata["slab_axis"]  = axis
 625        msh.metadata["slab_operation"] = operation
 626
 627        # compute bounds of slab
 628        origin = list(self.origin())
 629        spacing = list(self.spacing())
 630        if axis == 'x':
 631            msh.metadata["slab_bounding_box"] = [
 632                origin[0] + slice_range[0]*spacing[0],
 633                origin[0] + slice_range[1]*spacing[0],
 634                origin[1],
 635                origin[1] + dims[1]*spacing[1],
 636                origin[2],
 637                origin[2] + dims[2]*spacing[2],
 638            ]
 639        elif axis == 'y':
 640            msh.metadata["slab_bounding_box"] = [
 641                origin[0],
 642                origin[0] + dims[0]*spacing[0],
 643                origin[1] + slice_range[0]*spacing[1],
 644                origin[1] + slice_range[1]*spacing[1],
 645                origin[2],
 646                origin[2] + dims[2]*spacing[2],
 647            ]
 648        elif axis == 'z':
 649            msh.metadata["slab_bounding_box"] = [
 650                origin[0],
 651                origin[0] + dims[0]*spacing[0],
 652                origin[1],
 653                origin[1] + dims[1]*spacing[1],
 654                origin[2] + slice_range[0]*spacing[2],
 655                origin[2] + slice_range[1]*spacing[2],
 656            ]
 657
 658        msh.pipeline = utils.OperationNode(
 659            f"slab{slice_range}", 
 660            comment=f"axis={axis}, operation={operation}",
 661            parents=[self],
 662            c="#4cc9f0:#e9c46a",
 663        )
 664        msh.name = "SlabMesh"
 665        return msh
 666
 667
 668    def warp(
 669            self,
 670            source: Union["vedo.Points", List],
 671            target: Union["vedo.Points", List],
 672            sigma=1, mode="3d", fit=True,
 673        ) -> Self:
 674        """
 675        Warp volume scalars within a Volume by specifying
 676        source and target sets of points.
 677
 678        Arguments:
 679            source : (Points, list)
 680                the list of source points
 681            target : (Points, list)
 682                the list of target points
 683            fit : (bool)
 684                fit/adapt the old bounding box to the warped geometry
 685        """
 686        if isinstance(source, vedo.Points):
 687            source = source.vertices
 688        if isinstance(target, vedo.Points):
 689            target = target.vertices
 690
 691        NLT = transformations.NonLinearTransform()
 692        NLT.source_points = source
 693        NLT.target_points = target
 694        NLT.sigma = sigma
 695        NLT.mode = mode
 696
 697        self.apply_transform(NLT, fit=fit)
 698        self.pipeline = utils.OperationNode("warp", parents=[self], c="#4cc9f0")
 699        return self
 700
 701    def apply_transform(
 702            self, 
 703            T: Union[transformations.LinearTransform, transformations.NonLinearTransform],
 704            fit=True, interpolation="cubic",
 705        ) -> Self:
 706        """
 707        Apply a transform to the scalars in the volume.
 708
 709        Arguments:
 710            T : (LinearTransform, NonLinearTransform)
 711                The transformation to be applied
 712            fit : (bool)
 713                fit/adapt the old bounding box to the modified geometry
 714            interpolation : (str)
 715                one of the following: "nearest", "linear", "cubic"
 716        """
 717        if utils.is_sequence(T):
 718            T = transformations.LinearTransform(T)
 719
 720        TI = T.compute_inverse()
 721
 722        reslice = vtki.new("ImageReslice")
 723        reslice.SetInputData(self.dataset)
 724        reslice.SetResliceTransform(TI.T)
 725        reslice.SetOutputDimensionality(3)
 726        if "lin" in interpolation.lower():
 727            reslice.SetInterpolationModeToLinear()
 728        elif "near" in interpolation.lower():
 729            reslice.SetInterpolationModeToNearestNeighbor()
 730        elif "cubic" in interpolation.lower():
 731            reslice.SetInterpolationModeToCubic()
 732        else:
 733            vedo.logger.error(
 734                f"in apply_transform: unknown interpolation mode {interpolation}")
 735            raise ValueError()
 736        reslice.SetAutoCropOutput(fit)
 737        reslice.Update()
 738        self._update(reslice.GetOutput())
 739        self.transform = T
 740        self.pipeline = utils.OperationNode(
 741            "apply_transform", parents=[self], c="#4cc9f0")
 742        return self
 743
 744    def imagedata(self) -> vtki.vtkImageData:
 745        """
 746        DEPRECATED:
 747        Use `Volume.dataset` instead.
 748
 749        Return the underlying `vtkImagaData` object.
 750        """
 751        print("Volume.imagedata() is deprecated, use Volume.dataset instead")
 752        return self.dataset
 753    
 754    def modified(self) -> Self:
 755        """
 756        Mark the object as modified.
 757
 758        Example:
 759
 760        - [numpy2volume0.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/numpy2volume0.py)
 761        """
 762        scals = self.dataset.GetPointData().GetScalars()
 763        if scals:
 764            scals.Modified()
 765        return self
 766
 767    def tonumpy(self) -> np.ndarray:
 768        """
 769        Get read-write access to voxels of a Volume object as a numpy array.
 770
 771        When you set values in the output image, you don't want numpy to reallocate the array
 772        but instead set values in the existing array, so use the [:] operator.
 773
 774        Example:
 775            `arr[:] = arr*2 + 15`
 776
 777        If the array is modified add a call to:
 778        `volume.modified()`
 779        when all your modifications are completed.
 780        """
 781        narray_shape = tuple(reversed(self.dataset.GetDimensions()))
 782
 783        scals = self.dataset.GetPointData().GetScalars()
 784        comps = scals.GetNumberOfComponents()
 785        if comps == 1:
 786            narray = utils.vtk2numpy(scals).reshape(narray_shape)
 787            narray = np.transpose(narray, axes=[2, 1, 0])
 788        else:
 789            narray = utils.vtk2numpy(scals).reshape(*narray_shape, comps)
 790            narray = np.transpose(narray, axes=[2, 1, 0, 3])
 791
 792        # narray = utils.vtk2numpy(self.dataset.GetPointData().GetScalars()).reshape(narray_shape)
 793        # narray = np.transpose(narray, axes=[2, 1, 0])
 794        return narray
 795
 796    @property
 797    def shape(self) -> np.ndarray:
 798        """Return the nr. of voxels in the 3 dimensions."""
 799        return np.array(self.dataset.GetDimensions())
 800
 801    def dimensions(self) -> np.ndarray:
 802        """Return the nr. of voxels in the 3 dimensions."""
 803        return np.array(self.dataset.GetDimensions())
 804
 805    def scalar_range(self) -> np.ndarray:
 806        """Return the range of the scalar values."""
 807        return np.array(self.dataset.GetScalarRange())
 808
 809    def spacing(self, s=None) -> Union[Self, Iterable[float]]:
 810        """Set/get the voxels size in the 3 dimensions."""
 811        if s is not None:
 812            self.dataset.SetSpacing(s)
 813            return self
 814        return np.array(self.dataset.GetSpacing())
 815
 816    def origin(self, s=None) -> Union[Self, Iterable[float]]:
 817        """
 818        Set/get the origin of the volumetric dataset.
 819
 820        The origin is the position in world coordinates of the point index (0,0,0).
 821        This point does not have to be part of the dataset, in other words,
 822        the dataset extent does not have to start at (0,0,0) and the origin 
 823        can be outside of the dataset bounding box. 
 824        The origin plus spacing determine the position in space of the points.
 825        """
 826        if s is not None:
 827            self.dataset.SetOrigin(s)
 828            return self
 829        return np.array(self.dataset.GetOrigin())
 830    
 831    def pos(self, p=None) -> Union[Self, Iterable[float]]:
 832        """Set/get the position of the volumetric dataset."""
 833        if p is not None:
 834            self.origin(p)
 835            return self
 836        return self.origin()
 837
 838    def center(self) -> np.ndarray:
 839        """Get the center of the volumetric dataset."""
 840        # note that this does not have the set method like origin and spacing
 841        return np.array(self.dataset.GetCenter())
 842    
 843    def shift(self, s: list) -> Self:
 844        """Shift the volumetric dataset by a vector."""
 845        self.origin(self.origin() + np.array(s))
 846        return self
 847
 848    def rotate_x(self, angle: float, rad=False, around=None) -> Self:
 849        """
 850        Rotate around x-axis. If angle is in radians set `rad=True`.
 851
 852        Use `around` to define a pivoting point.
 853        """
 854        if angle == 0:
 855            return self
 856        LT = transformations.LinearTransform().rotate_x(angle, rad, around)
 857        return self.apply_transform(LT, fit=True, interpolation="linear")
 858
 859    def rotate_y(self, angle: float, rad=False, around=None) -> Self:
 860        """
 861        Rotate around y-axis. If angle is in radians set `rad=True`.
 862
 863        Use `around` to define a pivoting point.
 864        """
 865        if angle == 0:
 866            return self
 867        LT = transformations.LinearTransform().rotate_y(angle, rad, around)
 868        return self.apply_transform(LT, fit=True, interpolation="linear")
 869
 870    def rotate_z(self, angle: float, rad=False, around=None) -> Self:
 871        """
 872        Rotate around z-axis. If angle is in radians set `rad=True`.
 873
 874        Use `around` to define a pivoting point.
 875        """
 876        if angle == 0:
 877            return self
 878        LT = transformations.LinearTransform().rotate_z(angle, rad, around)
 879        return self.apply_transform(LT, fit=True, interpolation="linear")
 880
 881    def get_cell_from_ijk(self, ijk: list) -> int:
 882        """
 883        Get the voxel id number at the given ijk coordinates.
 884
 885        Arguments:
 886            ijk : (list)
 887                the ijk coordinates of the voxel
 888        """
 889        return self.dataset.ComputeCellId(ijk)
 890    
 891    def get_point_from_ijk(self, ijk: list) -> int:
 892        """
 893        Get the point id number at the given ijk coordinates.
 894
 895        Arguments:
 896            ijk : (list)
 897                the ijk coordinates of the voxel
 898        """
 899        return self.dataset.ComputePointId(ijk)
 900
 901    def permute_axes(self, x: int, y: int, z: int) -> Self:
 902        """
 903        Reorder the axes of the Volume by specifying
 904        the input axes which are supposed to become the new X, Y, and Z.
 905        """
 906        imp = vtki.new("ImagePermute")
 907        imp.SetFilteredAxes(x, y, z)
 908        imp.SetInputData(self.dataset)
 909        imp.Update()
 910        self._update(imp.GetOutput())
 911        self.pipeline = utils.OperationNode(
 912            f"permute_axes({(x,y,z)})", parents=[self], c="#4cc9f0"
 913        )
 914        return self
 915
 916    def resample(self, new_spacing: List[float], interpolation=1) -> Self:
 917        """
 918        Resamples a `Volume` to be larger or smaller.
 919
 920        This method modifies the spacing of the input.
 921        Linear interpolation is used to resample the data.
 922
 923        Arguments:
 924            new_spacing : (list)
 925                a list of 3 new spacings for the 3 axes
 926            interpolation : (int)
 927                0=nearest_neighbor, 1=linear, 2=cubic
 928        """
 929        rsp = vtki.new("ImageResample")
 930        oldsp = self.spacing()
 931        for i in range(3):
 932            if oldsp[i] != new_spacing[i]:
 933                rsp.SetAxisOutputSpacing(i, new_spacing[i])
 934        rsp.InterpolateOn()
 935        rsp.SetInterpolationMode(interpolation)
 936        rsp.OptimizationOn()
 937        rsp.Update()
 938        self._update(rsp.GetOutput())
 939        self.pipeline = utils.OperationNode(
 940            "resample", comment=f"spacing: {tuple(new_spacing)}", parents=[self], c="#4cc9f0"
 941        )
 942        return self
 943
 944
 945    def threshold(self, above=None, below=None, replace=None, replace_value=None) -> Self:
 946        """
 947        Binary or continuous volume thresholding.
 948        Find the voxels that contain a value above/below the input values
 949        and replace them with a new value (default is 0).
 950        """
 951        th = vtki.new("ImageThreshold")
 952        th.SetInputData(self.dataset)
 953
 954        # sanity checks
 955        if above is not None and below is not None:
 956            if above == below:
 957                return self
 958            if above > below:
 959                vedo.logger.warning("in volume.threshold(), above > below, skip.")
 960                return self
 961
 962        ## cases
 963        if below is not None and above is not None:
 964            th.ThresholdBetween(above, below)
 965
 966        elif above is not None:
 967            th.ThresholdByUpper(above)
 968
 969        elif below is not None:
 970            th.ThresholdByLower(below)
 971
 972        ##
 973        if replace is not None:
 974            th.SetReplaceIn(True)
 975            th.SetInValue(replace)
 976        else:
 977            th.SetReplaceIn(False)
 978
 979        if replace_value is not None:
 980            th.SetReplaceOut(True)
 981            th.SetOutValue(replace_value)
 982        else:
 983            th.SetReplaceOut(False)
 984
 985        th.Update()
 986        self._update(th.GetOutput())
 987        self.pipeline = utils.OperationNode("threshold", parents=[self], c="#4cc9f0")
 988        return self
 989
 990    def crop(self, left=None, right=None, back=None, front=None, bottom=None, top=None, VOI=()) -> Self:
 991        """
 992        Crop a `Volume` object.
 993
 994        Arguments:
 995            left : (float)
 996                fraction to crop from the left plane (negative x)
 997            right : (float)
 998                fraction to crop from the right plane (positive x)
 999            back : (float)
1000                fraction to crop from the back plane (negative y)
1001            front : (float)
1002                fraction to crop from the front plane (positive y)
1003            bottom : (float)
1004                fraction to crop from the bottom plane (negative z)
1005            top : (float)
1006                fraction to crop from the top plane (positive z)
1007            VOI : (list)
1008                extract Volume Of Interest expressed in voxel numbers
1009
1010        Example:
1011            `vol.crop(VOI=(xmin, xmax, ymin, ymax, zmin, zmax)) # all integers nrs`
1012        """
1013        extractVOI = vtki.new("ExtractVOI")
1014        extractVOI.SetInputData(self.dataset)
1015
1016        if VOI:
1017            extractVOI.SetVOI(VOI)
1018        else:
1019            d = self.dataset.GetDimensions()
1020            bx0, bx1, by0, by1, bz0, bz1 = 0, d[0]-1, 0, d[1]-1, 0, d[2]-1
1021            if left is not None:   bx0 = int((d[0]-1)*left)
1022            if right is not None:  bx1 = int((d[0]-1)*(1-right))
1023            if back is not None:   by0 = int((d[1]-1)*back)
1024            if front is not None:  by1 = int((d[1]-1)*(1-front))
1025            if bottom is not None: bz0 = int((d[2]-1)*bottom)
1026            if top is not None:    bz1 = int((d[2]-1)*(1-top))
1027            extractVOI.SetVOI(bx0, bx1, by0, by1, bz0, bz1)
1028        extractVOI.Update()
1029        self._update(extractVOI.GetOutput())
1030
1031        self.pipeline = utils.OperationNode(
1032            "crop", parents=[self], c="#4cc9f0", comment=f"dims={tuple(self.dimensions())}"
1033        )
1034        return self
1035
1036    def append(self, *volumes, axis="z", preserve_extents=False) -> Self:
1037        """
1038        Take the components from multiple inputs and merges them into one output.
1039        Except for the append axis, all inputs must have the same extent.
1040        All inputs must have the same number of scalar components.
1041        The output has the same origin and spacing as the first input.
1042        The origin and spacing of all other inputs are ignored.
1043        All inputs must have the same scalar type.
1044
1045        Arguments:
1046            axis : (int, str)
1047                axis expanded to hold the multiple images
1048            preserve_extents : (bool)
1049                if True, the extent of the inputs is used to place
1050                the image in the output. The whole extent of the output is the union of the input
1051                whole extents. Any portion of the output not covered by the inputs is set to zero.
1052                The origin and spacing is taken from the first input.
1053
1054        Example:
1055            ```python
1056            from vedo import Volume, dataurl
1057            vol = Volume(dataurl+'embryo.tif')
1058            vol.append(vol, axis='x').show().close()
1059            ```
1060            ![](https://vedo.embl.es/images/feats/volume_append.png)
1061        """
1062        ima = vtki.new("ImageAppend")
1063        ima.SetInputData(self.dataset)
1064        # if not utils.is_sequence(volumes):
1065        #     volumes = [volumes]
1066        for volume in volumes:
1067            if isinstance(volume, vtki.vtkImageData):
1068                ima.AddInputData(volume)
1069            else:
1070                ima.AddInputData(volume.dataset)
1071        ima.SetPreserveExtents(preserve_extents)
1072        if axis == "x":
1073            axis = 0
1074        elif axis == "y":
1075            axis = 1
1076        elif axis == "z":
1077            axis = 2
1078        ima.SetAppendAxis(axis)
1079        ima.Update()
1080        self._update(ima.GetOutput())
1081
1082        self.pipeline = utils.OperationNode(
1083            "append",
1084            parents=[self, *volumes],
1085            c="#4cc9f0",
1086            comment=f"dims={tuple(self.dimensions())}",
1087        )
1088        return self
1089
1090    def pad(self, voxels=10, value=0) -> Self:
1091        """
1092        Add the specified number of voxels at the `Volume` borders.
1093        Voxels can be a list formatted as `[nx0, nx1, ny0, ny1, nz0, nz1]`.
1094
1095        Arguments:
1096            voxels : (int, list)
1097                number of voxels to be added (or a list of length 4)
1098            value : (int)
1099                intensity value (gray-scale color) of the padding
1100
1101        Example:
1102            ```python
1103            from vedo import Volume, dataurl, show
1104            iso = Volume(dataurl+'embryo.tif').isosurface()
1105            vol = iso.binarize(spacing=(100, 100, 100)).pad(10)
1106            vol.dilate([15,15,15])
1107            show(iso, vol.isosurface(), N=2, axes=1)
1108            ```
1109            ![](https://vedo.embl.es/images/volumetric/volume_pad.png)
1110        """
1111        x0, x1, y0, y1, z0, z1 = self.dataset.GetExtent()
1112        pf = vtki.new("ImageConstantPad")
1113        pf.SetInputData(self.dataset)
1114        pf.SetConstant(value)
1115        if utils.is_sequence(voxels):
1116            pf.SetOutputWholeExtent(
1117                x0 - voxels[0], x1 + voxels[1],
1118                y0 - voxels[2], y1 + voxels[3],
1119                z0 - voxels[4], z1 + voxels[5],
1120            )
1121        else:
1122            pf.SetOutputWholeExtent(
1123                x0 - voxels, x1 + voxels,
1124                y0 - voxels, y1 + voxels,
1125                z0 - voxels, z1 + voxels,
1126            )
1127        pf.Update()
1128        self._update(pf.GetOutput())
1129        self.pipeline = utils.OperationNode(
1130            "pad", comment=f"{voxels} voxels", parents=[self], c="#f28482"
1131        )
1132        return self
1133
1134    def resize(self, newdims: List[int]) -> Self:
1135        """Increase or reduce the number of voxels of a Volume with interpolation."""
1136        rsz = vtki.new("ImageResize")
1137        rsz.SetResizeMethodToOutputDimensions()
1138        rsz.SetInputData(self.dataset)
1139        rsz.SetOutputDimensions(newdims)
1140        rsz.Update()
1141        self.dataset = rsz.GetOutput()
1142        self._update(self.dataset)
1143        self.pipeline = utils.OperationNode(
1144            "resize", parents=[self], c="#4cc9f0", comment=f"dims={tuple(self.dimensions())}"
1145        )
1146        return self
1147
1148    def normalize(self) -> Self:
1149        """Normalize that scalar components for each point."""
1150        norm = vtki.new("ImageNormalize")
1151        norm.SetInputData(self.dataset)
1152        norm.Update()
1153        self._update(norm.GetOutput())
1154        self.pipeline = utils.OperationNode("normalize", parents=[self], c="#4cc9f0")
1155        return self
1156
1157    def mirror(self, axis="x") -> Self:
1158        """
1159        Mirror flip along one of the cartesian axes.
1160        """
1161        img = self.dataset
1162
1163        ff = vtki.new("ImageFlip")
1164        ff.SetInputData(img)
1165        if axis.lower() == "x":
1166            ff.SetFilteredAxis(0)
1167        elif axis.lower() == "y":
1168            ff.SetFilteredAxis(1)
1169        elif axis.lower() == "z":
1170            ff.SetFilteredAxis(2)
1171        else:
1172            vedo.logger.error("mirror must be set to either x, y, z or n")
1173            raise RuntimeError()
1174        ff.Update()
1175        self._update(ff.GetOutput())
1176        self.pipeline = utils.OperationNode(f"mirror {axis}", parents=[self], c="#4cc9f0")
1177        return self
1178
1179    def operation(self, operation: str, volume2=None) -> "Volume":
1180        """
1181        Perform operations with `Volume` objects.
1182        Keyword `volume2` can be a constant `float`.
1183
1184        Possible operations are:
1185        ```
1186        and, or, xor, nand, nor, not,
1187        +, -, /, 1/x, sin, cos, exp, log,
1188        abs, **2, sqrt, min, max, atan, atan2, median,
1189        mag, dot, gradient, divergence, laplacian.
1190        ```
1191
1192        Example:
1193        ```py
1194        from vedo import Box, show
1195        vol1 = Box(size=(35,10, 5)).binarize()
1196        vol2 = Box(size=( 5,10,35)).binarize()
1197        vol = vol1.operation("xor", vol2)
1198        show([[vol1, vol2], 
1199            ["vol1 xor vol2", vol]],
1200            N=2, axes=1, viewup="z",
1201        ).close()
1202        ```
1203
1204        Note:
1205            For logic operations, the two volumes must have the same bounds.
1206            If they do not, a larger image is created to contain both and the
1207            volumes are resampled onto the larger image before the operation is
1208            performed. This can be slow and memory intensive.
1209
1210        See also:
1211            - [volume_operations.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/volume_operations.py)
1212        """
1213        op = operation.lower()
1214        image1 = self.dataset
1215
1216        if op in ["and", "or", "xor", "nand", "nor"]:
1217
1218            if not np.allclose(image1.GetBounds(), volume2.dataset.GetBounds()):
1219                # create a larger image to contain both
1220                b1 = image1.GetBounds()
1221                b2 = volume2.dataset.GetBounds()
1222                b = [
1223                    min(b1[0], b2[0]),
1224                    max(b1[1], b2[1]),
1225                    min(b1[2], b2[2]),
1226                    max(b1[3], b2[3]),
1227                    min(b1[4], b2[4]),
1228                    max(b1[5], b2[5]),
1229                ]
1230                dims1 = image1.GetDimensions()
1231                dims2 = volume2.dataset.GetDimensions()
1232                dims = [max(dims1[0], dims2[0]), max(dims1[1], dims2[1]), max(dims1[2], dims2[2])]
1233
1234                image = vtki.vtkImageData()
1235                image.SetDimensions(dims)
1236                spacing = (
1237                    (b[1] - b[0]) / dims[0],
1238                    (b[3] - b[2]) / dims[1],
1239                    (b[5] - b[4]) / dims[2],
1240                )
1241                image.SetSpacing(spacing)
1242                image.SetOrigin((b[0], b[2], b[4]))
1243                image.AllocateScalars(vtki.VTK_UNSIGNED_CHAR, 1)
1244                image.GetPointData().GetScalars().FillComponent(0, 0)
1245
1246                interp1 = vtki.new("ImageReslice")
1247                interp1.SetInputData(image1)
1248                interp1.SetOutputExtent(image.GetExtent())
1249                interp1.SetOutputOrigin(image.GetOrigin())
1250                interp1.SetOutputSpacing(image.GetSpacing())
1251                interp1.SetInterpolationModeToNearestNeighbor()
1252                interp1.Update()
1253                imageA = interp1.GetOutput()
1254
1255                interp2 = vtki.new("ImageReslice")
1256                interp2.SetInputData(volume2.dataset)
1257                interp2.SetOutputExtent(image.GetExtent())
1258                interp2.SetOutputOrigin(image.GetOrigin())
1259                interp2.SetOutputSpacing(image.GetSpacing())
1260                interp2.SetInterpolationModeToNearestNeighbor()
1261                interp2.Update()
1262                imageB = interp2.GetOutput()
1263
1264            else:
1265                imageA = image1
1266                imageB = volume2.dataset
1267
1268            img_logic = vtki.new("ImageLogic")
1269            img_logic.SetInput1Data(imageA)
1270            img_logic.SetInput2Data(imageB)
1271            img_logic.SetOperation(["and", "or", "xor", "nand", "nor"].index(op))
1272            img_logic.Update()
1273
1274            out_vol = Volume(img_logic.GetOutput())
1275            out_vol.pipeline = utils.OperationNode(
1276                "operation", comment=f"{op}", parents=[self, volume2], c="#4cc9f0", shape="cylinder"
1277            )
1278            return out_vol  ######################################################
1279
1280        if volume2 and isinstance(volume2, Volume):
1281            # assert image1.GetScalarType() == volume2.dataset.GetScalarType(), "volumes have different scalar types"
1282            # make sure they have the same bounds:
1283            assert np.allclose(image1.GetBounds(), volume2.dataset.GetBounds()), "volumes have different bounds"
1284            # make sure they have the same spacing:
1285            assert np.allclose(image1.GetSpacing(), volume2.dataset.GetSpacing()), "volumes have different spacing"
1286            # make sure they have the same origin:
1287            assert np.allclose(image1.GetOrigin(), volume2.dataset.GetOrigin()), "volumes have different origin"
1288
1289        mf = None
1290        if op in ["median"]:
1291            mf = vtki.new("ImageMedian3D")
1292            mf.SetInputData(image1)
1293        elif op in ["mag"]:
1294            mf = vtki.new("ImageMagnitude")
1295            mf.SetInputData(image1)
1296        elif op in ["dot"]:
1297            mf = vtki.new("ImageDotProduct")
1298            mf.SetInput1Data(image1)
1299            mf.SetInput2Data(volume2.dataset)
1300        elif op in ["grad", "gradient"]:
1301            mf = vtki.new("ImageGradient")
1302            mf.SetDimensionality(3)
1303            mf.SetInputData(image1)
1304        elif op in ["div", "divergence"]:
1305            mf = vtki.new("ImageDivergence")
1306            mf.SetInputData(image1)
1307        elif op in ["laplacian"]:
1308            mf = vtki.new("ImageLaplacian")
1309            mf.SetDimensionality(3)
1310            mf.SetInputData(image1)
1311        elif op in ["not"]:
1312            mf = vtki.new("ImageLogic")
1313            mf.SetInput1Data(image1)
1314            mf.SetOperation(4)
1315
1316        if mf is not None:
1317            mf.Update()
1318            vol = Volume(mf.GetOutput())
1319            vol.pipeline = utils.OperationNode(
1320                "operation", comment=f"{op}", parents=[self], c="#4cc9f0", shape="cylinder"
1321            )
1322            return vol  ######################################################
1323
1324        mat = vtki.new("ImageMathematics")
1325        mat.SetInput1Data(image1)
1326
1327        K = None
1328
1329        if utils.is_number(volume2):
1330            K = volume2
1331            mat.SetConstantK(K)
1332            mat.SetConstantC(K)
1333
1334        elif volume2 is not None:  # assume image2 is a constant value
1335            mat.SetInput2Data(volume2.dataset)
1336
1337        # ###########################
1338        if op in ["+", "add", "plus"]:
1339            if K:
1340                mat.SetOperationToAddConstant()
1341            else:
1342                mat.SetOperationToAdd()
1343
1344        elif op in ["-", "subtract", "minus"]:
1345            if K:
1346                mat.SetConstantC(-float(K))
1347                mat.SetOperationToAddConstant()
1348            else:
1349                mat.SetOperationToSubtract()
1350
1351        elif op in ["*", "multiply", "times"]:
1352            if K:
1353                mat.SetOperationToMultiplyByK()
1354            else:
1355                mat.SetOperationToMultiply()
1356
1357        elif op in ["/", "divide"]:
1358            if K:
1359                mat.SetConstantK(1.0 / K)
1360                mat.SetOperationToMultiplyByK()
1361            else:
1362                mat.SetOperationToDivide()
1363
1364        elif op in ["1/x", "invert"]:
1365            mat.SetOperationToInvert()
1366        elif op in ["sin"]:
1367            mat.SetOperationToSin()
1368        elif op in ["cos"]:
1369            mat.SetOperationToCos()
1370        elif op in ["exp"]:
1371            mat.SetOperationToExp()
1372        elif op in ["log"]:
1373            mat.SetOperationToLog()
1374        elif op in ["abs"]:
1375            mat.SetOperationToAbsoluteValue()
1376        elif op in ["**2", "square"]:
1377            mat.SetOperationToSquare()
1378        elif op in ["sqrt", "sqr"]:
1379            mat.SetOperationToSquareRoot()
1380        elif op in ["min"]:
1381            mat.SetOperationToMin()
1382        elif op in ["max"]:
1383            mat.SetOperationToMax()
1384        elif op in ["atan"]:
1385            mat.SetOperationToATAN()
1386        elif op in ["atan2"]:
1387            mat.SetOperationToATAN2()
1388        else:
1389            vedo.logger.error(f"unknown operation {operation}")
1390            raise RuntimeError()
1391        mat.Update()
1392
1393        self._update(mat.GetOutput())
1394
1395        self.pipeline = utils.OperationNode(
1396            "operation", comment=f"{op}", parents=[self, volume2], shape="cylinder", c="#4cc9f0"
1397        )
1398        return self
1399
1400    def frequency_pass_filter(self, low_cutoff=None, high_cutoff=None, order=1) -> Self:
1401        """
1402        Low-pass and high-pass filtering become trivial in the frequency domain.
1403        A portion of the pixels/voxels are simply masked or attenuated.
1404        This function applies a high pass Butterworth filter that attenuates the
1405        frequency domain image.
1406
1407        The gradual attenuation of the filter is important.
1408        A simple high-pass filter would simply mask a set of pixels in the frequency domain,
1409        but the abrupt transition would cause a ringing effect in the spatial domain.
1410
1411        Arguments:
1412            low_cutoff : (list)
1413                the cutoff frequencies for x, y and z
1414            high_cutoff : (list)
1415                the cutoff frequencies for x, y and z
1416            order : (int)
1417                order determines sharpness of the cutoff curve
1418        """
1419        # https://lorensen.github.io/VTKExamples/site/Cxx/ImageProcessing/IdealHighPass
1420        fft = vtki.new("ImageFFT")
1421        fft.SetInputData(self.dataset)
1422        fft.Update()
1423        out = fft.GetOutput()
1424
1425        if high_cutoff:
1426            blp = vtki.new("ImageButterworthLowPass")
1427            blp.SetInputData(out)
1428            blp.SetCutOff(high_cutoff)
1429            blp.SetOrder(order)
1430            blp.Update()
1431            out = blp.GetOutput()
1432
1433        if low_cutoff:
1434            bhp = vtki.new("ImageButterworthHighPass")
1435            bhp.SetInputData(out)
1436            bhp.SetCutOff(low_cutoff)
1437            bhp.SetOrder(order)
1438            bhp.Update()
1439            out = bhp.GetOutput()
1440
1441        rfft = vtki.new("ImageRFFT")
1442        rfft.SetInputData(out)
1443        rfft.Update()
1444
1445        ecomp = vtki.new("ImageExtractComponents")
1446        ecomp.SetInputData(rfft.GetOutput())
1447        ecomp.SetComponents(0)
1448        ecomp.Update()
1449        self._update(ecomp.GetOutput())
1450        self.pipeline = utils.OperationNode("frequency_pass_filter", parents=[self], c="#4cc9f0")
1451        return self
1452
1453    def smooth_gaussian(self, sigma=(2, 2, 2), radius=None) -> Self:
1454        """
1455        Performs a convolution of the input Volume with a gaussian.
1456
1457        Arguments:
1458            sigma : (float, list)
1459                standard deviation(s) in voxel units.
1460                A list can be given to smooth in the three direction differently.
1461            radius : (float, list)
1462                radius factor(s) determine how far out the gaussian
1463                kernel will go before being clamped to zero. A list can be given too.
1464        """
1465        gsf = vtki.new("ImageGaussianSmooth")
1466        gsf.SetDimensionality(3)
1467        gsf.SetInputData(self.dataset)
1468        if utils.is_sequence(sigma):
1469            gsf.SetStandardDeviations(sigma)
1470        else:
1471            gsf.SetStandardDeviation(sigma)
1472        if radius is not None:
1473            if utils.is_sequence(radius):
1474                gsf.SetRadiusFactors(radius)
1475            else:
1476                gsf.SetRadiusFactor(radius)
1477        gsf.Update()
1478        self._update(gsf.GetOutput())
1479        self.pipeline = utils.OperationNode("smooth_gaussian", parents=[self], c="#4cc9f0")
1480        return self
1481
1482    def smooth_median(self, neighbours=(2, 2, 2)) -> Self:
1483        """
1484        Median filter that replaces each pixel with the median value
1485        from a rectangular neighborhood around that pixel.
1486        """
1487        imgm = vtki.new("ImageMedian3D")
1488        imgm.SetInputData(self.dataset)
1489        if utils.is_sequence(neighbours):
1490            imgm.SetKernelSize(neighbours[0], neighbours[1], neighbours[2])
1491        else:
1492            imgm.SetKernelSize(neighbours, neighbours, neighbours)
1493        imgm.Update()
1494        self._update(imgm.GetOutput())
1495        self.pipeline = utils.OperationNode("smooth_median", parents=[self], c="#4cc9f0")
1496        return self
1497
1498    def erode(self, neighbours=(2, 2, 2)) -> Self:
1499        """
1500        Replace a voxel with the minimum over an ellipsoidal neighborhood of voxels.
1501        If `neighbours` of an axis is 1, no processing is done on that axis.
1502
1503        Examples:
1504            - [erode_dilate.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/erode_dilate.py)
1505
1506                ![](https://vedo.embl.es/images/volumetric/erode_dilate.png)
1507        """
1508        ver = vtki.new("ImageContinuousErode3D")
1509        ver.SetInputData(self.dataset)
1510        ver.SetKernelSize(neighbours[0], neighbours[1], neighbours[2])
1511        ver.Update()
1512        self._update(ver.GetOutput())
1513        self.pipeline = utils.OperationNode("erode", parents=[self], c="#4cc9f0")
1514        return self
1515
1516    def dilate(self, neighbours=(2, 2, 2)) -> Self:
1517        """
1518        Replace a voxel with the maximum over an ellipsoidal neighborhood of voxels.
1519        If `neighbours` of an axis is 1, no processing is done on that axis.
1520
1521        Check also `erode()` and `pad()`.
1522
1523        Examples:
1524            - [erode_dilate.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/erode_dilate.py)
1525        """
1526        ver = vtki.new("ImageContinuousDilate3D")
1527        ver.SetInputData(self.dataset)
1528        ver.SetKernelSize(neighbours[0], neighbours[1], neighbours[2])
1529        ver.Update()
1530        self._update(ver.GetOutput())
1531        self.pipeline = utils.OperationNode("dilate", parents=[self], c="#4cc9f0")
1532        return self
1533
1534    def magnitude(self) -> Self:
1535        """Colapses components with magnitude function."""
1536        imgm = vtki.new("ImageMagnitude")
1537        imgm.SetInputData(self.dataset)
1538        imgm.Update()
1539        self._update(imgm.GetOutput())
1540        self.pipeline = utils.OperationNode("magnitude", parents=[self], c="#4cc9f0")
1541        return self
1542
1543    def topoints(self) -> "vedo.Points":
1544        """
1545        Extract all image voxels as points.
1546        This function takes an input `Volume` and creates an `Mesh`
1547        that contains the points and the point attributes.
1548
1549        Examples:
1550            - [vol2points.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/vol2points.py)
1551        """
1552        v2p = vtki.new("ImageToPoints")
1553        v2p.SetInputData(self.dataset)
1554        v2p.Update()
1555        mpts = vedo.Points(v2p.GetOutput())
1556        mpts.pipeline = utils.OperationNode("topoints", parents=[self], c="#4cc9f0:#e9c46a")
1557        return mpts
1558
1559    def euclidean_distance(self, anisotropy=False, max_distance=None) -> "Volume":
1560        """
1561        Implementation of the Euclidean DT (Distance Transform) using Saito's algorithm.
1562        The distance map produced contains the square of the Euclidean distance values.
1563        The algorithm has a O(n^(D+1)) complexity over n x n x...x n images in D dimensions.
1564
1565        Check out also: https://en.wikipedia.org/wiki/Distance_transform
1566
1567        Arguments:
1568            anisotropy : bool
1569                used to define whether Spacing should be used in the
1570                computation of the distances.
1571            max_distance : bool
1572                any distance bigger than max_distance will not be
1573                computed but set to this specified value instead.
1574
1575        Examples:
1576            - [euclidian_dist.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/euclidian_dist.py)
1577        """
1578        euv = vtki.new("ImageEuclideanDistance")
1579        euv.SetInputData(self.dataset)
1580        euv.SetConsiderAnisotropy(anisotropy)
1581        if max_distance is not None:
1582            euv.InitializeOn()
1583            euv.SetMaximumDistance(max_distance)
1584        euv.SetAlgorithmToSaito()
1585        euv.Update()
1586        vol = Volume(euv.GetOutput())
1587        vol.pipeline = utils.OperationNode("euclidean_distance", parents=[self], c="#4cc9f0")
1588        return vol
1589
1590    def correlation_with(self, vol2: "Volume", dim=2) -> "Volume":
1591        """
1592        Find the correlation between two volumetric data sets.
1593        Keyword `dim` determines whether the correlation will be 3D, 2D or 1D.
1594        The default is a 2D Correlation.
1595
1596        The output size will match the size of the first input.
1597        The second input is considered the correlation kernel.
1598        """
1599        imc = vtki.new("ImageCorrelation")
1600        imc.SetInput1Data(self.dataset)
1601        imc.SetInput2Data(vol2.dataset)
1602        imc.SetDimensionality(dim)
1603        imc.Update()
1604        vol = Volume(imc.GetOutput())
1605
1606        vol.pipeline = utils.OperationNode("correlation_with", parents=[self, vol2], c="#4cc9f0")
1607        return vol
1608
1609    def scale_voxels(self, scale=1) -> Self:
1610        """Scale the voxel content by factor `scale`."""
1611        rsl = vtki.new("ImageReslice")
1612        rsl.SetInputData(self.dataset)
1613        rsl.SetScalarScale(scale)
1614        rsl.Update()
1615        self._update(rsl.GetOutput())
1616        self.pipeline = utils.OperationNode(
1617            "scale_voxels", comment=f"scale={scale}", parents=[self], c="#4cc9f0"
1618        )
1619        return self
  33class Volume(VolumeAlgorithms, VolumeVisual):
  34    """
  35    Class to describe dataset that are defined on "voxels",
  36    the 3D equivalent of 2D pixels.
  37    """
  38    def __init__(
  39        self,
  40        inputobj=None,
  41        dims=None,
  42        origin=None,
  43        spacing=None,
  44    ) -> None:
  45        """
  46        This class can be initialized with a numpy object,
  47        a `vtkImageData` or a list of 2D bmp files.
  48
  49        Arguments:
  50            origin : (list)
  51                set volume origin coordinates
  52            spacing : (list)
  53                voxel dimensions in x, y and z.
  54            dims : (list)
  55                specify the dimensions of the volume.
  56
  57        Example:
  58            ```python
  59            from vedo import Volume
  60            vol = Volume("path/to/mydata/rec*.bmp")
  61            vol.show()
  62            ```
  63
  64        Examples:
  65            - [numpy2volume1.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/numpy2volume1.py)
  66
  67                ![](https://vedo.embl.es/images/volumetric/numpy2volume1.png)
  68
  69            - [read_volume2.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/read_volume2.py)
  70
  71                ![](https://vedo.embl.es/images/volumetric/read_volume2.png)
  72
  73        .. note::
  74            if a `list` of values is used for `alphas` this is interpreted
  75            as a transfer function along the range of the scalar.
  76        """
  77        super().__init__()
  78
  79        self.name = "Volume"
  80        self.filename = ""
  81        self.file_size = ""
  82
  83        self.info = {}
  84        self.time =  time.time()
  85
  86        self.actor = vtki.vtkVolume()
  87        self.actor.retrieve_object = weak_ref_to(self)
  88        self.properties = self.actor.GetProperty()
  89
  90        self.transform = None
  91        self.point_locator = None
  92        self.cell_locator = None
  93        self.line_locator = None
  94
  95        ###################
  96        if isinstance(inputobj, str):
  97            if "https://" in inputobj:
  98                inputobj = vedo.file_io.download(inputobj, verbose=False)  # fpath
  99            elif os.path.isfile(inputobj):
 100                self.filename = inputobj
 101            else:
 102                inputobj = sorted(glob.glob(inputobj))
 103
 104        ###################
 105        inputtype = str(type(inputobj))
 106
 107        # print('Volume inputtype', inputtype, c='b')
 108
 109        if inputobj is None:
 110            img = vtki.vtkImageData()
 111
 112        elif utils.is_sequence(inputobj):
 113
 114            if isinstance(inputobj[0], str) and ".bmp" in inputobj[0].lower():
 115                # scan sequence of BMP files
 116                ima = vtki.new("ImageAppend")
 117                ima.SetAppendAxis(2)
 118                pb = utils.ProgressBar(0, len(inputobj))
 119                for i in pb.range():
 120                    f = inputobj[i]
 121                    if "_rec_spr" in f: # OPT specific
 122                        continue
 123                    picr = vtki.new("BMPReader")
 124                    picr.SetFileName(f)
 125                    picr.Update()
 126                    mgf = vtki.new("ImageMagnitude")
 127                    mgf.SetInputData(picr.GetOutput())
 128                    mgf.Update()
 129                    ima.AddInputData(mgf.GetOutput())
 130                    pb.print("loading...")
 131                ima.Update()
 132                img = ima.GetOutput()
 133
 134            else:
 135
 136                if len(inputobj.shape) == 1:
 137                    varr = utils.numpy2vtk(inputobj)
 138                else:
 139                    varr = utils.numpy2vtk(inputobj.ravel(order="F"))
 140                varr.SetName("input_scalars")
 141
 142                img = vtki.vtkImageData()
 143                if dims is not None:
 144                    img.SetDimensions(dims[2], dims[1], dims[0])
 145                else:
 146                    if len(inputobj.shape) == 1:
 147                        vedo.logger.error("must set dimensions (dims keyword) in Volume")
 148                        raise RuntimeError()
 149                    img.SetDimensions(inputobj.shape)
 150                img.GetPointData().AddArray(varr)
 151                img.GetPointData().SetActiveScalars(varr.GetName())
 152
 153        elif isinstance(inputobj, vtki.vtkImageData):
 154            img = inputobj
 155
 156        elif isinstance(inputobj, str):
 157            if "https://" in inputobj:
 158                inputobj = vedo.file_io.download(inputobj, verbose=False)
 159            img = vedo.file_io.loadImageData(inputobj)
 160            self.filename = inputobj
 161
 162        else:
 163            vedo.logger.error(f"cannot understand input type {inputtype}")
 164            return
 165
 166        if dims is not None:
 167            img.SetDimensions(dims)
 168
 169        if origin is not None:
 170            img.SetOrigin(origin)
 171
 172        if spacing is not None:
 173            img.SetSpacing(spacing)
 174
 175        self.dataset = img
 176        self.transform = None
 177
 178        #####################################
 179        mapper = vtki.new("SmartVolumeMapper")
 180        mapper.SetInputData(img)
 181        self.actor.SetMapper(mapper)
 182
 183        if img.GetPointData().GetScalars():
 184            if img.GetPointData().GetScalars().GetNumberOfComponents() == 1:
 185                self.properties.SetShade(True)
 186                self.properties.SetInterpolationType(1)
 187                self.cmap("RdBu_r")
 188                self.alpha([0.0, 0.0, 0.2, 0.4, 0.8, 1.0])
 189                self.alpha_gradient(None)
 190                self.properties.SetScalarOpacityUnitDistance(1.0)
 191
 192        self.pipeline = utils.OperationNode(
 193            "Volume", comment=f"dims={tuple(self.dimensions())}", c="#4cc9f0"
 194        )
 195        #######################################################################
 196
 197    @property
 198    def mapper(self):
 199        """Return the underlying `vtkMapper` object."""
 200        return self.actor.GetMapper()
 201    
 202    @mapper.setter
 203    def mapper(self, mapper):
 204        """
 205        Set the underlying `vtkMapper` object.
 206        
 207        Arguments:
 208            mapper : (str, vtkMapper)
 209                either 'gpu', 'opengl_gpu', 'fixed' or 'smart'
 210        """
 211        if isinstance(mapper, 
 212            (vtki.get_class("Mapper"), vtki.get_class("ImageResliceMapper"))
 213        ):
 214            pass
 215        elif mapper is None:
 216            mapper = vtki.new("SmartVolumeMapper")
 217        elif "gpu" in mapper:
 218            mapper = vtki.new("GPUVolumeRayCastMapper")
 219        elif "opengl_gpu" in mapper:
 220            mapper = vtki.new("OpenGLGPUVolumeRayCastMapper")
 221        elif "smart" in mapper:
 222            mapper = vtki.new("SmartVolumeMapper")
 223        elif "fixed" in mapper:
 224            mapper = vtki.new("FixedPointVolumeRayCastMapper")
 225        else:
 226            print("Error unknown mapper type", [mapper])
 227            raise RuntimeError()
 228        self.actor.SetMapper(mapper)
 229
 230    def c(self, *args, **kwargs) -> Self:
 231        """Deprecated. Use `Volume.cmap()` instead."""
 232        vedo.logger.warning("Volume.c() is deprecated, use Volume.cmap() instead")
 233        return self.cmap(*args, **kwargs)
 234
 235    def _update(self, data, reset_locators=False):
 236        # reset_locators here is dummy
 237        self.dataset = data
 238        self.mapper.SetInputData(data)
 239        self.dataset.GetPointData().Modified()
 240        self.mapper.Modified()
 241        self.mapper.Update()
 242        return self
 243
 244    def __str__(self):
 245        """Print a summary for the `Volume` object."""
 246        module = self.__class__.__module__
 247        name = self.__class__.__name__
 248        out = vedo.printc(
 249            f"{module}.{name} at ({hex(self.memory_address())})".ljust(75),
 250            c="c", bold=True, invert=True, return_string=True,
 251        )
 252        out += "\x1b[0m\x1b[36;1m"
 253
 254        out+= "name".ljust(14) + ": " + str(self.name) + "\n"
 255        if self.filename:
 256            out+= "filename".ljust(14) + ": " + str(self.filename) + "\n"
 257
 258        out+= "dimensions".ljust(14) + ": " + str(self.shape) + "\n"
 259
 260        out+= "origin".ljust(14) + ": "
 261        out+= utils.precision(self.origin(), 6) + "\n"
 262
 263        out+= "center".ljust(14) + ": "
 264        out+= utils.precision(self.center(), 6) + "\n"
 265
 266        out+= "spacing".ljust(14)    + ": "
 267        out+= utils.precision(self.spacing(), 6) + "\n"
 268
 269        bnds = self.bounds()
 270        bx1, bx2 = utils.precision(bnds[0], 3), utils.precision(bnds[1], 3)
 271        by1, by2 = utils.precision(bnds[2], 3), utils.precision(bnds[3], 3)
 272        bz1, bz2 = utils.precision(bnds[4], 3), utils.precision(bnds[5], 3)
 273        out+= "bounds".ljust(14) + ":"
 274        out+= " x=(" + bx1 + ", " + bx2 + "),"
 275        out+= " y=(" + by1 + ", " + by2 + "),"
 276        out+= " z=(" + bz1 + ", " + bz2 + ")\n"
 277
 278        out+= "memory size".ljust(14) + ": "
 279        out+= str(int(self.dataset.GetActualMemorySize()/1024+0.5))+" MB\n"
 280
 281        st = self.dataset.GetScalarTypeAsString()
 282        out+= "scalar size".ljust(14) + ": "
 283        out+= str(self.dataset.GetScalarSize()) + f" bytes ({st})\n"
 284        out+= "scalar range".ljust(14) + ": "
 285        out+= str(self.dataset.GetScalarRange()) + "\n"
 286
 287        #utils.print_histogram(self, logscale=True, bins=8, height=15, c="b", bold=True)
 288        return out.rstrip() + "\x1b[0m"
 289
 290    def _repr_html_(self):
 291        """
 292        HTML representation of the Volume object for Jupyter Notebooks.
 293
 294        Returns:
 295            HTML text with the image and some properties.
 296        """
 297        import io
 298        import base64
 299        from PIL import Image
 300
 301        library_name = "vedo.volume.Volume"
 302        help_url = "https://vedo.embl.es/docs/vedo/volume.html"
 303
 304        arr = self.thumbnail(azimuth=0, elevation=-60, zoom=1.4, axes=True)
 305
 306        im = Image.fromarray(arr)
 307        buffered = io.BytesIO()
 308        im.save(buffered, format="PNG", quality=100)
 309        encoded = base64.b64encode(buffered.getvalue()).decode("utf-8")
 310        url = "data:image/png;base64," + encoded
 311        image = f"<img src='{url}'></img>"
 312
 313        # statisitics
 314        bounds = "<br/>".join(
 315            [
 316                utils.precision(min_x, 4) + " ... " + utils.precision(max_x, 4)
 317                for min_x, max_x in zip(self.bounds()[::2], self.bounds()[1::2])
 318            ]
 319        )
 320
 321        help_text = ""
 322        if self.name:
 323            help_text += f"<b> {self.name}: &nbsp&nbsp</b>"
 324        help_text += '<b><a href="' + help_url + '" target="_blank">' + library_name + "</a></b>"
 325        if self.filename:
 326            dots = ""
 327            if len(self.filename) > 30:
 328                dots = "..."
 329            help_text += f"<br/><code><i>({dots}{self.filename[-30:]})</i></code>"
 330
 331        pdata = ""
 332        if self.dataset.GetPointData().GetScalars():
 333            if self.dataset.GetPointData().GetScalars().GetName():
 334                name = self.dataset.GetPointData().GetScalars().GetName()
 335                pdata = "<tr><td><b> point data array </b></td><td>" + name + "</td></tr>"
 336
 337        cdata = ""
 338        if self.dataset.GetCellData().GetScalars():
 339            if self.dataset.GetCellData().GetScalars().GetName():
 340                name = self.dataset.GetCellData().GetScalars().GetName()
 341                cdata = "<tr><td><b> voxel data array </b></td><td>" + name + "</td></tr>"
 342
 343        img = self.dataset
 344
 345        allt = [
 346            "<table>",
 347            "<tr>",
 348            "<td>",
 349            image,
 350            "</td>",
 351            "<td style='text-align: center; vertical-align: center;'><br/>",
 352            help_text,
 353            "<table>",
 354            "<tr><td><b> bounds </b> <br/> (x/y/z) </td><td>" + str(bounds) + "</td></tr>",
 355            "<tr><td><b> dimensions </b></td><td>" + str(img.GetDimensions()) + "</td></tr>",
 356            "<tr><td><b> voxel spacing </b></td><td>"
 357            + utils.precision(img.GetSpacing(), 3)
 358            + "</td></tr>",
 359            "<tr><td><b> in memory size </b></td><td>"
 360            + str(int(img.GetActualMemorySize() / 1024))
 361            + "MB</td></tr>",
 362            pdata,
 363            cdata,
 364            "<tr><td><b> scalar range </b></td><td>"
 365            + utils.precision(img.GetScalarRange(), 4)
 366            + "</td></tr>",
 367            "</table>",
 368            "</table>",
 369        ]
 370        return "\n".join(allt)
 371
 372    def copy(self, deep=True) -> "Volume":
 373        """Return a copy of the Volume. Alias of `clone()`."""
 374        return self.clone(deep=deep)
 375
 376    def clone(self, deep=True) -> "Volume":
 377        """Return a clone copy of the Volume. Alias of `copy()`."""
 378        if deep:
 379            newimg = vtki.vtkImageData()
 380            newimg.CopyStructure(self.dataset)
 381            newimg.CopyAttributes(self.dataset)
 382            newvol = Volume(newimg)
 383        else:
 384            newvol = Volume(self.dataset)
 385
 386        prop = vtki.vtkVolumeProperty()
 387        prop.DeepCopy(self.properties)
 388        newvol.actor.SetProperty(prop)
 389        newvol.properties = prop
 390
 391        newvol.pipeline = utils.OperationNode("clone", parents=[self], c="#bbd0ff", shape="diamond")
 392        return newvol
 393
 394    def astype(self, dtype: Union[str, int]) -> Self:
 395        """
 396        Reset the type of the scalars array.
 397
 398        Arguments:
 399            dtype : (str)
 400                the type of the scalars array in
 401                ["int8", "uint8", "int16", "uint16", "int32", "uint32", "float32", "float64"]
 402        """
 403        if dtype in ["int8", "uint8", "int16", "uint16", "int32", "uint32", "float32", "float64"]:
 404            caster = vtki.new("ImageCast")
 405            caster.SetInputData(self.dataset)
 406            caster.SetOutputScalarType(int(vtki.array_types[dtype]))
 407            caster.ClampOverflowOn()
 408            caster.Update()
 409            self._update(caster.GetOutput())
 410            self.pipeline = utils.OperationNode(f"astype({dtype})", parents=[self], c="#4cc9f0")
 411        else:
 412            vedo.logger.error(f"astype(): unknown type {dtype}")
 413            raise ValueError()
 414        return self
 415
 416    def component_weight(self, i: int, weight: float) -> Self:
 417        """Set the scalar component weight in range [0,1]."""
 418        self.properties.SetComponentWeight(i, weight)
 419        return self
 420
 421    def xslice(self, i: int) -> Mesh:
 422        """Extract the slice at index `i` of volume along x-axis."""
 423        vslice = vtki.new("ImageDataGeometryFilter")
 424        vslice.SetInputData(self.dataset)
 425        nx, ny, nz = self.dataset.GetDimensions()
 426        if i > nx - 1:
 427            i = nx - 1
 428        vslice.SetExtent(i, i, 0, ny, 0, nz)
 429        vslice.Update()
 430        m = Mesh(vslice.GetOutput())
 431        m.pipeline = utils.OperationNode(f"xslice {i}", parents=[self], c="#4cc9f0:#e9c46a")
 432        return m
 433
 434    def yslice(self, j: int) -> Mesh:
 435        """Extract the slice at index `j` of volume along y-axis."""
 436        vslice = vtki.new("ImageDataGeometryFilter")
 437        vslice.SetInputData(self.dataset)
 438        nx, ny, nz = self.dataset.GetDimensions()
 439        if j > ny - 1:
 440            j = ny - 1
 441        vslice.SetExtent(0, nx, j, j, 0, nz)
 442        vslice.Update()
 443        m = Mesh(vslice.GetOutput())
 444        m.pipeline = utils.OperationNode(f"yslice {j}", parents=[self], c="#4cc9f0:#e9c46a")
 445        return m
 446
 447    def zslice(self, k: int) -> Mesh:
 448        """Extract the slice at index `i` of volume along z-axis."""
 449        vslice = vtki.new("ImageDataGeometryFilter")
 450        vslice.SetInputData(self.dataset)
 451        nx, ny, nz = self.dataset.GetDimensions()
 452        if k > nz - 1:
 453            k = nz - 1
 454        vslice.SetExtent(0, nx, 0, ny, k, k)
 455        vslice.Update()
 456        m = Mesh(vslice.GetOutput())
 457        m.pipeline = utils.OperationNode(f"zslice {k}", parents=[self], c="#4cc9f0:#e9c46a")
 458        return m
 459
 460    def slice_plane(self, origin: List[float], normal: List[float], autocrop=False, border=0.5, mode="linear") -> Mesh:
 461        """
 462        Extract the slice along a given plane position and normal.
 463
 464        Two metadata arrays are added to the output Mesh:
 465            - "shape" : contains the shape of the slice
 466            - "original_bounds" : contains the original bounds of the slice
 467        One can access them with e.g. `myslice.metadata["shape"]`.
 468
 469        Arguments:
 470            origin : (list)
 471                position of the plane
 472            normal : (list)
 473                normal to the plane
 474            autocrop : (bool)
 475                crop the output to the minimal possible size
 476            border : (float)
 477                add a border to the output slice
 478            mode : (str)
 479                interpolation mode, one of the following: "linear", "nearest", "cubic"
 480
 481        Example:
 482            - [slice_plane1.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/slice_plane1.py)
 483
 484                ![](https://vedo.embl.es/images/volumetric/slicePlane1.gif)
 485            
 486            - [slice_plane2.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/slice_plane2.py)
 487                
 488                ![](https://vedo.embl.es/images/volumetric/slicePlane2.png)
 489
 490            - [slice_plane3.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/slice_plane3.py)
 491
 492                ![](https://vedo.embl.es/images/volumetric/slicePlane3.jpg)
 493        """
 494        newaxis = utils.versor(normal)
 495        pos = np.array(origin)
 496        initaxis = (0, 0, 1)
 497        crossvec = np.cross(initaxis, newaxis)
 498        angle = np.arccos(np.dot(initaxis, newaxis))
 499        T = vtki.vtkTransform()
 500        T.PostMultiply()
 501        T.RotateWXYZ(np.rad2deg(angle), crossvec.tolist())
 502        T.Translate(pos.tolist())
 503
 504        reslice = vtki.new("ImageReslice")
 505        reslice.SetResliceAxes(T.GetMatrix())
 506        reslice.SetInputData(self.dataset)
 507        reslice.SetOutputDimensionality(2)
 508        reslice.SetTransformInputSampling(True)
 509        reslice.SetGenerateStencilOutput(False)
 510        if border:
 511            reslice.SetBorder(True)
 512            reslice.SetBorderThickness(border)
 513        else:
 514            reslice.SetBorder(False)
 515        if mode == "linear":
 516            reslice.SetInterpolationModeToLinear()
 517        elif mode == "nearest":
 518            reslice.SetInterpolationModeToNearestNeighbor()
 519        elif mode == "cubic":
 520            reslice.SetInterpolationModeToCubic()
 521        else:
 522            vedo.logger.error(f"in slice_plane(): unknown interpolation mode {mode}")
 523            raise ValueError()
 524        reslice.SetAutoCropOutput(not autocrop)
 525        reslice.Update()
 526        img = reslice.GetOutput()
 527
 528        vslice = vtki.new("ImageDataGeometryFilter")
 529        vslice.SetInputData(img)
 530        vslice.Update()
 531
 532        msh = Mesh(vslice.GetOutput()).apply_transform(T)
 533        msh.properties.LightingOff()
 534
 535        d0, d1, _ = img.GetDimensions()
 536        varr1 = utils.numpy2vtk([d1, d0], name="shape")
 537        varr2 = utils.numpy2vtk(img.GetBounds(), name="original_bounds")
 538        msh.dataset.GetFieldData().AddArray(varr1)
 539        msh.dataset.GetFieldData().AddArray(varr2)
 540        msh.pipeline = utils.OperationNode("slice_plane", parents=[self], c="#4cc9f0:#e9c46a")
 541        return msh
 542    
 543    def slab(self, slice_range=(), axis='z', operation="mean") -> Mesh:
 544        """
 545        Extract a slab from a `Volume` by combining 
 546        all of the slices of an image to create a single slice.
 547
 548        Returns a `Mesh` containing metadata which
 549        can be accessed with e.g. `mesh.metadata["slab_range"]`.
 550
 551        Metadata:
 552            slab_range : (list)
 553                contains the range of slices extracted
 554            slab_axis : (str)
 555                contains the axis along which the slab was extracted
 556            slab_operation : (str)
 557                contains the operation performed on the slab
 558            slab_bounding_box : (list)
 559                contains the bounding box of the slab
 560
 561        Arguments:
 562            slice_range : (list)
 563                range of slices to extract
 564            axis : (str)
 565                axis along which to extract the slab
 566            operation : (str)
 567                operation to perform on the slab,
 568                allowed values are: "sum", "min", "max", "mean".
 569        
 570        Example:
 571            - [slab.py](https://github.com/marcomusy/vedo/blob/master/examples/volumetric/slab_vol.py)
 572
 573            ![](https://vedo.embl.es/images/volumetric/slab_vol.jpg)
 574        """
 575        if len(slice_range) != 2:
 576            vedo.logger.error("in slab(): slice_range is empty or invalid")
 577            raise ValueError()
 578        
 579        islab = vtki.new("ImageSlab")
 580        islab.SetInputData(self.dataset)
 581
 582        if operation in ["+", "add", "sum"]:
 583            islab.SetOperationToSum()
 584        elif "min" in operation:
 585            islab.SetOperationToMin()
 586        elif "max" in operation:
 587            islab.SetOperationToMax()
 588        elif "mean" in operation:
 589            islab.SetOperationToMean()
 590        else:
 591            vedo.logger.error(f"in slab(): unknown operation {operation}")
 592            raise ValueError()
 593
 594        dims = self.dimensions()
 595        if axis == 'x':
 596            islab.SetOrientationToX()
 597            if slice_range[0]  > dims[0]-1:
 598                slice_range[0] = dims[0]-1
 599            if slice_range[1]  > dims[0]-1:
 600                slice_range[1] = dims[0]-1
 601        elif axis == 'y':
 602            islab.SetOrientationToY()
 603            if slice_range[0]  > dims[1]-1:
 604                slice_range[0] = dims[1]-1
 605            if slice_range[1]  > dims[1]-1:
 606                slice_range[1] = dims[1]-1
 607        elif axis == 'z':
 608            islab.SetOrientationToZ()
 609            if slice_range[0]  > dims[2]-1:
 610                slice_range[0] = dims[2]-1
 611            if slice_range[1]  > dims[2]-1:
 612                slice_range[1] = dims[2]-1
 613        else:
 614            vedo.logger.error(f"Error in slab(): unknown axis {axis}")
 615            raise RuntimeError()
 616        
 617        islab.SetSliceRange(slice_range)
 618        islab.Update()
 619
 620        msh = Mesh(islab.GetOutput()).lighting('off')
 621        msh.mapper.SetLookupTable(utils.ctf2lut(self, msh))
 622        msh.mapper.SetScalarRange(self.scalar_range())
 623
 624        msh.metadata["slab_range"] = slice_range
 625        msh.metadata["slab_axis"]  = axis
 626        msh.metadata["slab_operation"] = operation
 627
 628        # compute bounds of slab
 629        origin = list(self.origin())
 630        spacing = list(self.spacing())
 631        if axis == 'x':
 632            msh.metadata["slab_bounding_box"] = [
 633                origin[0] + slice_range[0]*spacing[0],
 634                origin[0] + slice_range[1]*spacing[0],
 635                origin[1],
 636                origin[1] + dims[1]*spacing[1],
 637                origin[2],
 638                origin[2] + dims[2]*spacing[2],
 639            ]
 640        elif axis == 'y':
 641            msh.metadata["slab_bounding_box"] = [
 642                origin[0],
 643                origin[0] + dims[0]*spacing[0],
 644                origin[1] + slice_range[0]*spacing[1],
 645                origin[1] + slice_range[1]*spacing[1],
 646                origin[2],
 647                origin[2] + dims[2]*spacing[2],
 648            ]
 649        elif axis == 'z':
 650            msh.metadata["slab_bounding_box"] = [
 651                origin[0],
 652                origin[0] + dims[0]*spacing[0],
 653                origin[1],
 654                origin[1] + dims[1]*spacing[1],
 655                origin[2] + slice_range[0]*spacing[2],
 656                origin[2] + slice_range[1]*spacing[2],
 657            ]
 658
 659        msh.pipeline = utils.OperationNode(
 660            f"slab{slice_range}", 
 661            comment=f"axis={axis}, operation={operation}",
 662            parents=[self],
 663            c="#4cc9f0:#e9c46a",
 664        )
 665        msh.name = "SlabMesh"
 666        return msh
 667
 668
 669    def warp(
 670            self,
 671            source: Union["vedo.Points", List],
 672            target: Union["vedo.Points", List],
 673            sigma=1, mode="3d", fit=True,
 674        ) -> Self:
 675        """
 676        Warp volume scalars within a Volume by specifying
 677        source and target sets of points.
 678
 679        Arguments:
 680            source : (Points, list)
 681                the list of source points
 682            target : (Points, list)
 683                the list of target points
 684            fit : (bool)
 685                fit/adapt the old bounding box to the warped geometry
 686        """
 687        if isinstance(source, vedo.Points):
 688            source = source.vertices
 689        if isinstance(target, vedo.Points):
 690            target = target.vertices
 691
 692        NLT = transformations.NonLinearTransform()
 693        NLT.source_points = source
 694        NLT.target_points = target
 695        NLT.sigma = sigma
 696        NLT.mode = mode
 697
 698        self.apply_transform(NLT, fit=fit)
 699        self.pipeline = utils.OperationNode("warp", parents=[self], c="#4cc9f0")
 700        return self
 701
 702    def apply_transform(
 703            self, 
 704            T: Union[transformations.LinearTransform, transformations.NonLinearTransform],
 705            fit=True, interpolation="cubic",
 706        ) -> Self:
 707        """
 708        Apply a transform to the scalars in the volume.
 709
 710        Arguments:
 711            T : (LinearTransform, NonLinearTransform)
 712                The transformation to be applied
 713            fit : (bool)
 714                fit/adapt the old bounding box to the modified geometry
 715            interpolation : (str)
 716                one of the following: "nearest", "linear", "cubic"
 717        """
 718        if utils.is_sequence(T):
 719            T = transformations.LinearTransform(T)
 720
 721        TI = T.compute_inverse()
 722
 723        reslice = vtki.new("ImageReslice")
 724        reslice.SetInputData(self.dataset)
 725        reslice.SetResliceTransform(TI.T)
 726        reslice.SetOutputDimensionality(3)
 727        if "lin" in interpolation.lower():
 728            reslice.SetInterpolationModeToLinear()
 729        elif "near" in interpolation.lower():
 730            reslice.SetInterpolationModeToNearestNeighbor()
 731        elif "cubic" in interpolation.lower():
 732            reslice.SetInterpolationModeToCubic()
 733        else:
 734            vedo.logger.error(
 735                f"in apply_transform: unknown interpolation mode {interpolation}")
 736            raise ValueError()
 737        reslice.SetAutoCropOutput(fit)
 738        reslice.Update()
 739        self._update(reslice.GetOutput())
 740        self.transform = T
 741        self.pipeline = utils.OperationNode(
 742            "apply_transform", parents=[self], c="#4cc9f0")
 743        return self
 744
 745    def imagedata(self) -> vtki.vtkImageData:
 746        """
 747        DEPRECATED:
 748        Use `Volume.dataset` instead.
 749
 750        Return the underlying `vtkImagaData` object.
 751        """
 752        print("Volume.imagedata() is deprecated, use Volume.dataset instead")
 753        return self.dataset
 754    
 755    def modified(self) -> Self:
 756        """
 757        Mark the object as modified.
 758
 759        Example:
 760
 761        - [numpy2volume0.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/numpy2volume0.py)
 762        """
 763        scals = self.dataset.GetPointData().GetScalars()
 764        if scals:
 765            scals.Modified()
 766        return self
 767
 768    def tonumpy(self) -> np.ndarray:
 769        """
 770        Get read-write access to voxels of a Volume object as a numpy array.
 771
 772        When you set values in the output image, you don't want numpy to reallocate the array
 773        but instead set values in the existing array, so use the [:] operator.
 774
 775        Example:
 776            `arr[:] = arr*2 + 15`
 777
 778        If the array is modified add a call to:
 779        `volume.modified()`
 780        when all your modifications are completed.
 781        """
 782        narray_shape = tuple(reversed(self.dataset.GetDimensions()))
 783
 784        scals = self.dataset.GetPointData().GetScalars()
 785        comps = scals.GetNumberOfComponents()
 786        if comps == 1:
 787            narray = utils.vtk2numpy(scals).reshape(narray_shape)
 788            narray = np.transpose(narray, axes=[2, 1, 0])
 789        else:
 790            narray = utils.vtk2numpy(scals).reshape(*narray_shape, comps)
 791            narray = np.transpose(narray, axes=[2, 1, 0, 3])
 792
 793        # narray = utils.vtk2numpy(self.dataset.GetPointData().GetScalars()).reshape(narray_shape)
 794        # narray = np.transpose(narray, axes=[2, 1, 0])
 795        return narray
 796
 797    @property
 798    def shape(self) -> np.ndarray:
 799        """Return the nr. of voxels in the 3 dimensions."""
 800        return np.array(self.dataset.GetDimensions())
 801
 802    def dimensions(self) -> np.ndarray:
 803        """Return the nr. of voxels in the 3 dimensions."""
 804        return np.array(self.dataset.GetDimensions())
 805
 806    def scalar_range(self) -> np.ndarray:
 807        """Return the range of the scalar values."""
 808        return np.array(self.dataset.GetScalarRange())
 809
 810    def spacing(self, s=None) -> Union[Self, Iterable[float]]:
 811        """Set/get the voxels size in the 3 dimensions."""
 812        if s is not None:
 813            self.dataset.SetSpacing(s)
 814            return self
 815        return np.array(self.dataset.GetSpacing())
 816
 817    def origin(self, s=None) -> Union[Self, Iterable[float]]:
 818        """
 819        Set/get the origin of the volumetric dataset.
 820
 821        The origin is the position in world coordinates of the point index (0,0,0).
 822        This point does not have to be part of the dataset, in other words,
 823        the dataset extent does not have to start at (0,0,0) and the origin 
 824        can be outside of the dataset bounding box. 
 825        The origin plus spacing determine the position in space of the points.
 826        """
 827        if s is not None:
 828            self.dataset.SetOrigin(s)
 829            return self
 830        return np.array(self.dataset.GetOrigin())
 831    
 832    def pos(self, p=None) -> Union[Self, Iterable[float]]:
 833        """Set/get the position of the volumetric dataset."""
 834        if p is not None:
 835            self.origin(p)
 836            return self
 837        return self.origin()
 838
 839    def center(self) -> np.ndarray:
 840        """Get the center of the volumetric dataset."""
 841        # note that this does not have the set method like origin and spacing
 842        return np.array(self.dataset.GetCenter())
 843    
 844    def shift(self, s: list) -> Self:
 845        """Shift the volumetric dataset by a vector."""
 846        self.origin(self.origin() + np.array(s))
 847        return self
 848
 849    def rotate_x(self, angle: float, rad=False, around=None) -> Self:
 850        """
 851        Rotate around x-axis. If angle is in radians set `rad=True`.
 852
 853        Use `around` to define a pivoting point.
 854        """
 855        if angle == 0:
 856            return self
 857        LT = transformations.LinearTransform().rotate_x(angle, rad, around)
 858        return self.apply_transform(LT, fit=True, interpolation="linear")
 859
 860    def rotate_y(self, angle: float, rad=False, around=None) -> Self:
 861        """
 862        Rotate around y-axis. If angle is in radians set `rad=True`.
 863
 864        Use `around` to define a pivoting point.
 865        """
 866        if angle == 0:
 867            return self
 868        LT = transformations.LinearTransform().rotate_y(angle, rad, around)
 869        return self.apply_transform(LT, fit=True, interpolation="linear")
 870
 871    def rotate_z(self, angle: float, rad=False, around=None) -> Self:
 872        """
 873        Rotate around z-axis. If angle is in radians set `rad=True`.
 874
 875        Use `around` to define a pivoting point.
 876        """
 877        if angle == 0:
 878            return self
 879        LT = transformations.LinearTransform().rotate_z(angle, rad, around)
 880        return self.apply_transform(LT, fit=True, interpolation="linear")
 881
 882    def get_cell_from_ijk(self, ijk: list) -> int:
 883        """
 884        Get the voxel id number at the given ijk coordinates.
 885
 886        Arguments:
 887            ijk : (list)
 888                the ijk coordinates of the voxel
 889        """
 890        return self.dataset.ComputeCellId(ijk)
 891    
 892    def get_point_from_ijk(self, ijk: list) -> int:
 893        """
 894        Get the point id number at the given ijk coordinates.
 895
 896        Arguments:
 897            ijk : (list)
 898                the ijk coordinates of the voxel
 899        """
 900        return self.dataset.ComputePointId(ijk)
 901
 902    def permute_axes(self, x: int, y: int, z: int) -> Self:
 903        """
 904        Reorder the axes of the Volume by specifying
 905        the input axes which are supposed to become the new X, Y, and Z.
 906        """
 907        imp = vtki.new("ImagePermute")
 908        imp.SetFilteredAxes(x, y, z)
 909        imp.SetInputData(self.dataset)
 910        imp.Update()
 911        self._update(imp.GetOutput())
 912        self.pipeline = utils.OperationNode(
 913            f"permute_axes({(x,y,z)})", parents=[self], c="#4cc9f0"
 914        )
 915        return self
 916
 917    def resample(self, new_spacing: List[float], interpolation=1) -> Self:
 918        """
 919        Resamples a `Volume` to be larger or smaller.
 920
 921        This method modifies the spacing of the input.
 922        Linear interpolation is used to resample the data.
 923
 924        Arguments:
 925            new_spacing : (list)
 926                a list of 3 new spacings for the 3 axes
 927            interpolation : (int)
 928                0=nearest_neighbor, 1=linear, 2=cubic
 929        """
 930        rsp = vtki.new("ImageResample")
 931        oldsp = self.spacing()
 932        for i in range(3):
 933            if oldsp[i] != new_spacing[i]:
 934                rsp.SetAxisOutputSpacing(i, new_spacing[i])
 935        rsp.InterpolateOn()
 936        rsp.SetInterpolationMode(interpolation)
 937        rsp.OptimizationOn()
 938        rsp.Update()
 939        self._update(rsp.GetOutput())
 940        self.pipeline = utils.OperationNode(
 941            "resample", comment=f"spacing: {tuple(new_spacing)}", parents=[self], c="#4cc9f0"
 942        )
 943        return self
 944
 945
 946    def threshold(self, above=None, below=None, replace=None, replace_value=None) -> Self:
 947        """
 948        Binary or continuous volume thresholding.
 949        Find the voxels that contain a value above/below the input values
 950        and replace them with a new value (default is 0).
 951        """
 952        th = vtki.new("ImageThreshold")
 953        th.SetInputData(self.dataset)
 954
 955        # sanity checks
 956        if above is not None and below is not None:
 957            if above == below:
 958                return self
 959            if above > below:
 960                vedo.logger.warning("in volume.threshold(), above > below, skip.")
 961                return self
 962
 963        ## cases
 964        if below is not None and above is not None:
 965            th.ThresholdBetween(above, below)
 966
 967        elif above is not None:
 968            th.ThresholdByUpper(above)
 969
 970        elif below is not None:
 971            th.ThresholdByLower(below)
 972
 973        ##
 974        if replace is not None:
 975            th.SetReplaceIn(True)
 976            th.SetInValue(replace)
 977        else:
 978            th.SetReplaceIn(False)
 979
 980        if replace_value is not None:
 981            th.SetReplaceOut(True)
 982            th.SetOutValue(replace_value)
 983        else:
 984            th.SetReplaceOut(False)
 985
 986        th.Update()
 987        self._update(th.GetOutput())
 988        self.pipeline = utils.OperationNode("threshold", parents=[self], c="#4cc9f0")
 989        return self
 990
 991    def crop(self, left=None, right=None, back=None, front=None, bottom=None, top=None, VOI=()) -> Self:
 992        """
 993        Crop a `Volume` object.
 994
 995        Arguments:
 996            left : (float)
 997                fraction to crop from the left plane (negative x)
 998            right : (float)
 999                fraction to crop from the right plane (positive x)
1000            back : (float)
1001                fraction to crop from the back plane (negative y)
1002            front : (float)
1003                fraction to crop from the front plane (positive y)
1004            bottom : (float)
1005                fraction to crop from the bottom plane (negative z)
1006            top : (float)
1007                fraction to crop from the top plane (positive z)
1008            VOI : (list)
1009                extract Volume Of Interest expressed in voxel numbers
1010
1011        Example:
1012            `vol.crop(VOI=(xmin, xmax, ymin, ymax, zmin, zmax)) # all integers nrs`
1013        """
1014        extractVOI = vtki.new("ExtractVOI")
1015        extractVOI.SetInputData(self.dataset)
1016
1017        if VOI:
1018            extractVOI.SetVOI(VOI)
1019        else:
1020            d = self.dataset.GetDimensions()
1021            bx0, bx1, by0, by1, bz0, bz1 = 0, d[0]-1, 0, d[1]-1, 0, d[2]-1
1022            if left is not None:   bx0 = int((d[0]-1)*left)
1023            if right is not None:  bx1 = int((d[0]-1)*(1-right))
1024            if back is not None:   by0 = int((d[1]-1)*back)
1025            if front is not None:  by1 = int((d[1]-1)*(1-front))
1026            if bottom is not None: bz0 = int((d[2]-1)*bottom)
1027            if top is not None:    bz1 = int((d[2]-1)*(1-top))
1028            extractVOI.SetVOI(bx0, bx1, by0, by1, bz0, bz1)
1029        extractVOI.Update()
1030        self._update(extractVOI.GetOutput())
1031
1032        self.pipeline = utils.OperationNode(
1033            "crop", parents=[self], c="#4cc9f0", comment=f"dims={tuple(self.dimensions())}"
1034        )
1035        return self
1036
1037    def append(self, *volumes, axis="z", preserve_extents=False) -> Self:
1038        """
1039        Take the components from multiple inputs and merges them into one output.
1040        Except for the append axis, all inputs must have the same extent.
1041        All inputs must have the same number of scalar components.
1042        The output has the same origin and spacing as the first input.
1043        The origin and spacing of all other inputs are ignored.
1044        All inputs must have the same scalar type.
1045
1046        Arguments:
1047            axis : (int, str)
1048                axis expanded to hold the multiple images
1049            preserve_extents : (bool)
1050                if True, the extent of the inputs is used to place
1051                the image in the output. The whole extent of the output is the union of the input
1052                whole extents. Any portion of the output not covered by the inputs is set to zero.
1053                The origin and spacing is taken from the first input.
1054
1055        Example:
1056            ```python
1057            from vedo import Volume, dataurl
1058            vol = Volume(dataurl+'embryo.tif')
1059            vol.append(vol, axis='x').show().close()
1060            ```
1061            ![](https://vedo.embl.es/images/feats/volume_append.png)
1062        """
1063        ima = vtki.new("ImageAppend")
1064        ima.SetInputData(self.dataset)
1065        # if not utils.is_sequence(volumes):
1066        #     volumes = [volumes]
1067        for volume in volumes:
1068            if isinstance(volume, vtki.vtkImageData):
1069                ima.AddInputData(volume)
1070            else:
1071                ima.AddInputData(volume.dataset)
1072        ima.SetPreserveExtents(preserve_extents)
1073        if axis == "x":
1074            axis = 0
1075        elif axis == "y":
1076            axis = 1
1077        elif axis == "z":
1078            axis = 2
1079        ima.SetAppendAxis(axis)
1080        ima.Update()
1081        self._update(ima.GetOutput())
1082
1083        self.pipeline = utils.OperationNode(
1084            "append",
1085            parents=[self, *volumes],
1086            c="#4cc9f0",
1087            comment=f"dims={tuple(self.dimensions())}",
1088        )
1089        return self
1090
1091    def pad(self, voxels=10, value=0) -> Self:
1092        """
1093        Add the specified number of voxels at the `Volume` borders.
1094        Voxels can be a list formatted as `[nx0, nx1, ny0, ny1, nz0, nz1]`.
1095
1096        Arguments:
1097            voxels : (int, list)
1098                number of voxels to be added (or a list of length 4)
1099            value : (int)
1100                intensity value (gray-scale color) of the padding
1101
1102        Example:
1103            ```python
1104            from vedo import Volume, dataurl, show
1105            iso = Volume(dataurl+'embryo.tif').isosurface()
1106            vol = iso.binarize(spacing=(100, 100, 100)).pad(10)
1107            vol.dilate([15,15,15])
1108            show(iso, vol.isosurface(), N=2, axes=1)
1109            ```
1110            ![](https://vedo.embl.es/images/volumetric/volume_pad.png)
1111        """
1112        x0, x1, y0, y1, z0, z1 = self.dataset.GetExtent()
1113        pf = vtki.new("ImageConstantPad")
1114        pf.SetInputData(self.dataset)
1115        pf.SetConstant(value)
1116        if utils.is_sequence(voxels):
1117            pf.SetOutputWholeExtent(
1118                x0 - voxels[0], x1 + voxels[1],
1119                y0 - voxels[2], y1 + voxels[3],
1120                z0 - voxels[4], z1 + voxels[5],
1121            )
1122        else:
1123            pf.SetOutputWholeExtent(
1124                x0 - voxels, x1 + voxels,
1125                y0 - voxels, y1 + voxels,
1126                z0 - voxels, z1 + voxels,
1127            )
1128        pf.Update()
1129        self._update(pf.GetOutput())
1130        self.pipeline = utils.OperationNode(
1131            "pad", comment=f"{voxels} voxels", parents=[self], c="#f28482"
1132        )
1133        return self
1134
1135    def resize(self, newdims: List[int]) -> Self:
1136        """Increase or reduce the number of voxels of a Volume with interpolation."""
1137        rsz = vtki.new("ImageResize")
1138        rsz.SetResizeMethodToOutputDimensions()
1139        rsz.SetInputData(self.dataset)
1140        rsz.SetOutputDimensions(newdims)
1141        rsz.Update()
1142        self.dataset = rsz.GetOutput()
1143        self._update(self.dataset)
1144        self.pipeline = utils.OperationNode(
1145            "resize", parents=[self], c="#4cc9f0", comment=f"dims={tuple(self.dimensions())}"
1146        )
1147        return self
1148
1149    def normalize(self) -> Self:
1150        """Normalize that scalar components for each point."""
1151        norm = vtki.new("ImageNormalize")
1152        norm.SetInputData(self.dataset)
1153        norm.Update()
1154        self._update(norm.GetOutput())
1155        self.pipeline = utils.OperationNode("normalize", parents=[self], c="#4cc9f0")
1156        return self
1157
1158    def mirror(self, axis="x") -> Self:
1159        """
1160        Mirror flip along one of the cartesian axes.
1161        """
1162        img = self.dataset
1163
1164        ff = vtki.new("ImageFlip")
1165        ff.SetInputData(img)
1166        if axis.lower() == "x":
1167            ff.SetFilteredAxis(0)
1168        elif axis.lower() == "y":
1169            ff.SetFilteredAxis(1)
1170        elif axis.lower() == "z":
1171            ff.SetFilteredAxis(2)
1172        else:
1173            vedo.logger.error("mirror must be set to either x, y, z or n")
1174            raise RuntimeError()
1175        ff.Update()
1176        self._update(ff.GetOutput())
1177        self.pipeline = utils.OperationNode(f"mirror {axis}", parents=[self], c="#4cc9f0")
1178        return self
1179
1180    def operation(self, operation: str, volume2=None) -> "Volume":
1181        """
1182        Perform operations with `Volume` objects.
1183        Keyword `volume2` can be a constant `float`.
1184
1185        Possible operations are:
1186        ```
1187        and, or, xor, nand, nor, not,
1188        +, -, /, 1/x, sin, cos, exp, log,
1189        abs, **2, sqrt, min, max, atan, atan2, median,
1190        mag, dot, gradient, divergence, laplacian.
1191        ```
1192
1193        Example:
1194        ```py
1195        from vedo import Box, show
1196        vol1 = Box(size=(35,10, 5)).binarize()
1197        vol2 = Box(size=( 5,10,35)).binarize()
1198        vol = vol1.operation("xor", vol2)
1199        show([[vol1, vol2], 
1200            ["vol1 xor vol2", vol]],
1201            N=2, axes=1, viewup="z",
1202        ).close()
1203        ```
1204
1205        Note:
1206            For logic operations, the two volumes must have the same bounds.
1207            If they do not, a larger image is created to contain both and the
1208            volumes are resampled onto the larger image before the operation is
1209            performed. This can be slow and memory intensive.
1210
1211        See also:
1212            - [volume_operations.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/volume_operations.py)
1213        """
1214        op = operation.lower()
1215        image1 = self.dataset
1216
1217        if op in ["and", "or", "xor", "nand", "nor"]:
1218
1219            if not np.allclose(image1.GetBounds(), volume2.dataset.GetBounds()):
1220                # create a larger image to contain both
1221                b1 = image1.GetBounds()
1222                b2 = volume2.dataset.GetBounds()
1223                b = [
1224                    min(b1[0], b2[0]),
1225                    max(b1[1], b2[1]),
1226                    min(b1[2], b2[2]),
1227                    max(b1[3], b2[3]),
1228                    min(b1[4], b2[4]),
1229                    max(b1[5], b2[5]),
1230                ]
1231                dims1 = image1.GetDimensions()
1232                dims2 = volume2.dataset.GetDimensions()
1233                dims = [max(dims1[0], dims2[0]), max(dims1[1], dims2[1]), max(dims1[2], dims2[2])]
1234
1235                image = vtki.vtkImageData()
1236                image.SetDimensions(dims)
1237                spacing = (
1238                    (b[1] - b[0]) / dims[0],
1239                    (b[3] - b[2]) / dims[1],
1240                    (b[5] - b[4]) / dims[2],
1241                )
1242                image.SetSpacing(spacing)
1243                image.SetOrigin((b[0], b[2], b[4]))
1244                image.AllocateScalars(vtki.VTK_UNSIGNED_CHAR, 1)
1245                image.GetPointData().GetScalars().FillComponent(0, 0)
1246
1247                interp1 = vtki.new("ImageReslice")
1248                interp1.SetInputData(image1)
1249                interp1.SetOutputExtent(image.GetExtent())
1250                interp1.SetOutputOrigin(image.GetOrigin())
1251                interp1.SetOutputSpacing(image.GetSpacing())
1252                interp1.SetInterpolationModeToNearestNeighbor()
1253                interp1.Update()
1254                imageA = interp1.GetOutput()
1255
1256                interp2 = vtki.new("ImageReslice")
1257                interp2.SetInputData(volume2.dataset)
1258                interp2.SetOutputExtent(image.GetExtent())
1259                interp2.SetOutputOrigin(image.GetOrigin())
1260                interp2.SetOutputSpacing(image.GetSpacing())
1261                interp2.SetInterpolationModeToNearestNeighbor()
1262                interp2.Update()
1263                imageB = interp2.GetOutput()
1264
1265            else:
1266                imageA = image1
1267                imageB = volume2.dataset
1268
1269            img_logic = vtki.new("ImageLogic")
1270            img_logic.SetInput1Data(imageA)
1271            img_logic.SetInput2Data(imageB)
1272            img_logic.SetOperation(["and", "or", "xor", "nand", "nor"].index(op))
1273            img_logic.Update()
1274
1275            out_vol = Volume(img_logic.GetOutput())
1276            out_vol.pipeline = utils.OperationNode(
1277                "operation", comment=f"{op}", parents=[self, volume2], c="#4cc9f0", shape="cylinder"
1278            )
1279            return out_vol  ######################################################
1280
1281        if volume2 and isinstance(volume2, Volume):
1282            # assert image1.GetScalarType() == volume2.dataset.GetScalarType(), "volumes have different scalar types"
1283            # make sure they have the same bounds:
1284            assert np.allclose(image1.GetBounds(), volume2.dataset.GetBounds()), "volumes have different bounds"
1285            # make sure they have the same spacing:
1286            assert np.allclose(image1.GetSpacing(), volume2.dataset.GetSpacing()), "volumes have different spacing"
1287            # make sure they have the same origin:
1288            assert np.allclose(image1.GetOrigin(), volume2.dataset.GetOrigin()), "volumes have different origin"
1289
1290        mf = None
1291        if op in ["median"]:
1292            mf = vtki.new("ImageMedian3D")
1293            mf.SetInputData(image1)
1294        elif op in ["mag"]:
1295            mf = vtki.new("ImageMagnitude")
1296            mf.SetInputData(image1)
1297        elif op in ["dot"]:
1298            mf = vtki.new("ImageDotProduct")
1299            mf.SetInput1Data(image1)
1300            mf.SetInput2Data(volume2.dataset)
1301        elif op in ["grad", "gradient"]:
1302            mf = vtki.new("ImageGradient")
1303            mf.SetDimensionality(3)
1304            mf.SetInputData(image1)
1305        elif op in ["div", "divergence"]:
1306            mf = vtki.new("ImageDivergence")
1307            mf.SetInputData(image1)
1308        elif op in ["laplacian"]:
1309            mf = vtki.new("ImageLaplacian")
1310            mf.SetDimensionality(3)
1311            mf.SetInputData(image1)
1312        elif op in ["not"]:
1313            mf = vtki.new("ImageLogic")
1314            mf.SetInput1Data(image1)
1315            mf.SetOperation(4)
1316
1317        if mf is not None:
1318            mf.Update()
1319            vol = Volume(mf.GetOutput())
1320            vol.pipeline = utils.OperationNode(
1321                "operation", comment=f"{op}", parents=[self], c="#4cc9f0", shape="cylinder"
1322            )
1323            return vol  ######################################################
1324
1325        mat = vtki.new("ImageMathematics")
1326        mat.SetInput1Data(image1)
1327
1328        K = None
1329
1330        if utils.is_number(volume2):
1331            K = volume2
1332            mat.SetConstantK(K)
1333            mat.SetConstantC(K)
1334
1335        elif volume2 is not None:  # assume image2 is a constant value
1336            mat.SetInput2Data(volume2.dataset)
1337
1338        # ###########################
1339        if op in ["+", "add", "plus"]:
1340            if K:
1341                mat.SetOperationToAddConstant()
1342            else:
1343                mat.SetOperationToAdd()
1344
1345        elif op in ["-", "subtract", "minus"]:
1346            if K:
1347                mat.SetConstantC(-float(K))
1348                mat.SetOperationToAddConstant()
1349            else:
1350                mat.SetOperationToSubtract()
1351
1352        elif op in ["*", "multiply", "times"]:
1353            if K:
1354                mat.SetOperationToMultiplyByK()
1355            else:
1356                mat.SetOperationToMultiply()
1357
1358        elif op in ["/", "divide"]:
1359            if K:
1360                mat.SetConstantK(1.0 / K)
1361                mat.SetOperationToMultiplyByK()
1362            else:
1363                mat.SetOperationToDivide()
1364
1365        elif op in ["1/x", "invert"]:
1366            mat.SetOperationToInvert()
1367        elif op in ["sin"]:
1368            mat.SetOperationToSin()
1369        elif op in ["cos"]:
1370            mat.SetOperationToCos()
1371        elif op in ["exp"]:
1372            mat.SetOperationToExp()
1373        elif op in ["log"]:
1374            mat.SetOperationToLog()
1375        elif op in ["abs"]:
1376            mat.SetOperationToAbsoluteValue()
1377        elif op in ["**2", "square"]:
1378            mat.SetOperationToSquare()
1379        elif op in ["sqrt", "sqr"]:
1380            mat.SetOperationToSquareRoot()
1381        elif op in ["min"]:
1382            mat.SetOperationToMin()
1383        elif op in ["max"]:
1384            mat.SetOperationToMax()
1385        elif op in ["atan"]:
1386            mat.SetOperationToATAN()
1387        elif op in ["atan2"]:
1388            mat.SetOperationToATAN2()
1389        else:
1390            vedo.logger.error(f"unknown operation {operation}")
1391            raise RuntimeError()
1392        mat.Update()
1393
1394        self._update(mat.GetOutput())
1395
1396        self.pipeline = utils.OperationNode(
1397            "operation", comment=f"{op}", parents=[self, volume2], shape="cylinder", c="#4cc9f0"
1398        )
1399        return self
1400
1401    def frequency_pass_filter(self, low_cutoff=None, high_cutoff=None, order=1) -> Self:
1402        """
1403        Low-pass and high-pass filtering become trivial in the frequency domain.
1404        A portion of the pixels/voxels are simply masked or attenuated.
1405        This function applies a high pass Butterworth filter that attenuates the
1406        frequency domain image.
1407
1408        The gradual attenuation of the filter is important.
1409        A simple high-pass filter would simply mask a set of pixels in the frequency domain,
1410        but the abrupt transition would cause a ringing effect in the spatial domain.
1411
1412        Arguments:
1413            low_cutoff : (list)
1414                the cutoff frequencies for x, y and z
1415            high_cutoff : (list)
1416                the cutoff frequencies for x, y and z
1417            order : (int)
1418                order determines sharpness of the cutoff curve
1419        """
1420        # https://lorensen.github.io/VTKExamples/site/Cxx/ImageProcessing/IdealHighPass
1421        fft = vtki.new("ImageFFT")
1422        fft.SetInputData(self.dataset)
1423        fft.Update()
1424        out = fft.GetOutput()
1425
1426        if high_cutoff:
1427            blp = vtki.new("ImageButterworthLowPass")
1428            blp.SetInputData(out)
1429            blp.SetCutOff(high_cutoff)
1430            blp.SetOrder(order)
1431            blp.Update()
1432            out = blp.GetOutput()
1433
1434        if low_cutoff:
1435            bhp = vtki.new("ImageButterworthHighPass")
1436            bhp.SetInputData(out)
1437            bhp.SetCutOff(low_cutoff)
1438            bhp.SetOrder(order)
1439            bhp.Update()
1440            out = bhp.GetOutput()
1441
1442        rfft = vtki.new("ImageRFFT")
1443        rfft.SetInputData(out)
1444        rfft.Update()
1445
1446        ecomp = vtki.new("ImageExtractComponents")
1447        ecomp.SetInputData(rfft.GetOutput())
1448        ecomp.SetComponents(0)
1449        ecomp.Update()
1450        self._update(ecomp.GetOutput())
1451        self.pipeline = utils.OperationNode("frequency_pass_filter", parents=[self], c="#4cc9f0")
1452        return self
1453
1454    def smooth_gaussian(self, sigma=(2, 2, 2), radius=None) -> Self:
1455        """
1456        Performs a convolution of the input Volume with a gaussian.
1457
1458        Arguments:
1459            sigma : (float, list)
1460                standard deviation(s) in voxel units.
1461                A list can be given to smooth in the three direction differently.
1462            radius : (float, list)
1463                radius factor(s) determine how far out the gaussian
1464                kernel will go before being clamped to zero. A list can be given too.
1465        """
1466        gsf = vtki.new("ImageGaussianSmooth")
1467        gsf.SetDimensionality(3)
1468        gsf.SetInputData(self.dataset)
1469        if utils.is_sequence(sigma):
1470            gsf.SetStandardDeviations(sigma)
1471        else:
1472            gsf.SetStandardDeviation(sigma)
1473        if radius is not None:
1474            if utils.is_sequence(radius):
1475                gsf.SetRadiusFactors(radius)
1476            else:
1477                gsf.SetRadiusFactor(radius)
1478        gsf.Update()
1479        self._update(gsf.GetOutput())
1480        self.pipeline = utils.OperationNode("smooth_gaussian", parents=[self], c="#4cc9f0")
1481        return self
1482
1483    def smooth_median(self, neighbours=(2, 2, 2)) -> Self:
1484        """
1485        Median filter that replaces each pixel with the median value
1486        from a rectangular neighborhood around that pixel.
1487        """
1488        imgm = vtki.new("ImageMedian3D")
1489        imgm.SetInputData(self.dataset)
1490        if utils.is_sequence(neighbours):
1491            imgm.SetKernelSize(neighbours[0], neighbours[1], neighbours[2])
1492        else:
1493            imgm.SetKernelSize(neighbours, neighbours, neighbours)
1494        imgm.Update()
1495        self._update(imgm.GetOutput())
1496        self.pipeline = utils.OperationNode("smooth_median", parents=[self], c="#4cc9f0")
1497        return self
1498
1499    def erode(self, neighbours=(2, 2, 2)) -> Self:
1500        """
1501        Replace a voxel with the minimum over an ellipsoidal neighborhood of voxels.
1502        If `neighbours` of an axis is 1, no processing is done on that axis.
1503
1504        Examples:
1505            - [erode_dilate.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/erode_dilate.py)
1506
1507                ![](https://vedo.embl.es/images/volumetric/erode_dilate.png)
1508        """
1509        ver = vtki.new("ImageContinuousErode3D")
1510        ver.SetInputData(self.dataset)
1511        ver.SetKernelSize(neighbours[0], neighbours[1], neighbours[2])
1512        ver.Update()
1513        self._update(ver.GetOutput())
1514        self.pipeline = utils.OperationNode("erode", parents=[self], c="#4cc9f0")
1515        return self
1516
1517    def dilate(self, neighbours=(2, 2, 2)) -> Self:
1518        """
1519        Replace a voxel with the maximum over an ellipsoidal neighborhood of voxels.
1520        If `neighbours` of an axis is 1, no processing is done on that axis.
1521
1522        Check also `erode()` and `pad()`.
1523
1524        Examples:
1525            - [erode_dilate.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/erode_dilate.py)
1526        """
1527        ver = vtki.new("ImageContinuousDilate3D")
1528        ver.SetInputData(self.dataset)
1529        ver.SetKernelSize(neighbours[0], neighbours[1], neighbours[2])
1530        ver.Update()
1531        self._update(ver.GetOutput())
1532        self.pipeline = utils.OperationNode("dilate", parents=[self], c="#4cc9f0")
1533        return self
1534
1535    def magnitude(self) -> Self:
1536        """Colapses components with magnitude function."""
1537        imgm = vtki.new("ImageMagnitude")
1538        imgm.SetInputData(self.dataset)
1539        imgm.Update()
1540        self._update(imgm.GetOutput())
1541        self.pipeline = utils.OperationNode("magnitude", parents=[self], c="#4cc9f0")
1542        return self
1543
1544    def topoints(self) -> "vedo.Points":
1545        """
1546        Extract all image voxels as points.
1547        This function takes an input `Volume` and creates an `Mesh`
1548        that contains the points and the point attributes.
1549
1550        Examples:
1551            - [vol2points.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/vol2points.py)
1552        """
1553        v2p = vtki.new("ImageToPoints")
1554        v2p.SetInputData(self.dataset)
1555        v2p.Update()
1556        mpts = vedo.Points(v2p.GetOutput())
1557        mpts.pipeline = utils.OperationNode("topoints", parents=[self], c="#4cc9f0:#e9c46a")
1558        return mpts
1559
1560    def euclidean_distance(self, anisotropy=False, max_distance=None) -> "Volume":
1561        """
1562        Implementation of the Euclidean DT (Distance Transform) using Saito's algorithm.
1563        The distance map produced contains the square of the Euclidean distance values.
1564        The algorithm has a O(n^(D+1)) complexity over n x n x...x n images in D dimensions.
1565
1566        Check out also: https://en.wikipedia.org/wiki/Distance_transform
1567
1568        Arguments:
1569            anisotropy : bool
1570                used to define whether Spacing should be used in the
1571                computation of the distances.
1572            max_distance : bool
1573                any distance bigger than max_distance will not be
1574                computed but set to this specified value instead.
1575
1576        Examples:
1577            - [euclidian_dist.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/euclidian_dist.py)
1578        """
1579        euv = vtki.new("ImageEuclideanDistance")
1580        euv.SetInputData(self.dataset)
1581        euv.SetConsiderAnisotropy(anisotropy)
1582        if max_distance is not None:
1583            euv.InitializeOn()
1584            euv.SetMaximumDistance(max_distance)
1585        euv.SetAlgorithmToSaito()
1586        euv.Update()
1587        vol = Volume(euv.GetOutput())
1588        vol.pipeline = utils.OperationNode("euclidean_distance", parents=[self], c="#4cc9f0")
1589        return vol
1590
1591    def correlation_with(self, vol2: "Volume", dim=2) -> "Volume":
1592        """
1593        Find the correlation between two volumetric data sets.
1594        Keyword `dim` determines whether the correlation will be 3D, 2D or 1D.
1595        The default is a 2D Correlation.
1596
1597        The output size will match the size of the first input.
1598        The second input is considered the correlation kernel.
1599        """
1600        imc = vtki.new("ImageCorrelation")
1601        imc.SetInput1Data(self.dataset)
1602        imc.SetInput2Data(vol2.dataset)
1603        imc.SetDimensionality(dim)
1604        imc.Update()
1605        vol = Volume(imc.GetOutput())
1606
1607        vol.pipeline = utils.OperationNode("correlation_with", parents=[self, vol2], c="#4cc9f0")
1608        return vol
1609
1610    def scale_voxels(self, scale=1) -> Self:
1611        """Scale the voxel content by factor `scale`."""
1612        rsl = vtki.new("ImageReslice")
1613        rsl.SetInputData(self.dataset)
1614        rsl.SetScalarScale(scale)
1615        rsl.Update()
1616        self._update(rsl.GetOutput())
1617        self.pipeline = utils.OperationNode(
1618            "scale_voxels", comment=f"scale={scale}", parents=[self], c="#4cc9f0"
1619        )
1620        return self

Class to describe dataset that are defined on "voxels", the 3D equivalent of 2D pixels.

Volume(inputobj=None, dims=None, origin=None, spacing=None)
 38    def __init__(
 39        self,
 40        inputobj=None,
 41        dims=None,
 42        origin=None,
 43        spacing=None,
 44    ) -> None:
 45        """
 46        This class can be initialized with a numpy object,
 47        a `vtkImageData` or a list of 2D bmp files.
 48
 49        Arguments:
 50            origin : (list)
 51                set volume origin coordinates
 52            spacing : (list)
 53                voxel dimensions in x, y and z.
 54            dims : (list)
 55                specify the dimensions of the volume.
 56
 57        Example:
 58            ```python
 59            from vedo import Volume
 60            vol = Volume("path/to/mydata/rec*.bmp")
 61            vol.show()
 62            ```
 63
 64        Examples:
 65            - [numpy2volume1.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/numpy2volume1.py)
 66
 67                ![](https://vedo.embl.es/images/volumetric/numpy2volume1.png)
 68
 69            - [read_volume2.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/read_volume2.py)
 70
 71                ![](https://vedo.embl.es/images/volumetric/read_volume2.png)
 72
 73        .. note::
 74            if a `list` of values is used for `alphas` this is interpreted
 75            as a transfer function along the range of the scalar.
 76        """
 77        super().__init__()
 78
 79        self.name = "Volume"
 80        self.filename = ""
 81        self.file_size = ""
 82
 83        self.info = {}
 84        self.time =  time.time()
 85
 86        self.actor = vtki.vtkVolume()
 87        self.actor.retrieve_object = weak_ref_to(self)
 88        self.properties = self.actor.GetProperty()
 89
 90        self.transform = None
 91        self.point_locator = None
 92        self.cell_locator = None
 93        self.line_locator = None
 94
 95        ###################
 96        if isinstance(inputobj, str):
 97            if "https://" in inputobj:
 98                inputobj = vedo.file_io.download(inputobj, verbose=False)  # fpath
 99            elif os.path.isfile(inputobj):
100                self.filename = inputobj
101            else:
102                inputobj = sorted(glob.glob(inputobj))
103
104        ###################
105        inputtype = str(type(inputobj))
106
107        # print('Volume inputtype', inputtype, c='b')
108
109        if inputobj is None:
110            img = vtki.vtkImageData()
111
112        elif utils.is_sequence(inputobj):
113
114            if isinstance(inputobj[0], str) and ".bmp" in inputobj[0].lower():
115                # scan sequence of BMP files
116                ima = vtki.new("ImageAppend")
117                ima.SetAppendAxis(2)
118                pb = utils.ProgressBar(0, len(inputobj))
119                for i in pb.range():
120                    f = inputobj[i]
121                    if "_rec_spr" in f: # OPT specific
122                        continue
123                    picr = vtki.new("BMPReader")
124                    picr.SetFileName(f)
125                    picr.Update()
126                    mgf = vtki.new("ImageMagnitude")
127                    mgf.SetInputData(picr.GetOutput())
128                    mgf.Update()
129                    ima.AddInputData(mgf.GetOutput())
130                    pb.print("loading...")
131                ima.Update()
132                img = ima.GetOutput()
133
134            else:
135
136                if len(inputobj.shape) == 1:
137                    varr = utils.numpy2vtk(inputobj)
138                else:
139                    varr = utils.numpy2vtk(inputobj.ravel(order="F"))
140                varr.SetName("input_scalars")
141
142                img = vtki.vtkImageData()
143                if dims is not None:
144                    img.SetDimensions(dims[2], dims[1], dims[0])
145                else:
146                    if len(inputobj.shape) == 1:
147                        vedo.logger.error("must set dimensions (dims keyword) in Volume")
148                        raise RuntimeError()
149                    img.SetDimensions(inputobj.shape)
150                img.GetPointData().AddArray(varr)
151                img.GetPointData().SetActiveScalars(varr.GetName())
152
153        elif isinstance(inputobj, vtki.vtkImageData):
154            img = inputobj
155
156        elif isinstance(inputobj, str):
157            if "https://" in inputobj:
158                inputobj = vedo.file_io.download(inputobj, verbose=False)
159            img = vedo.file_io.loadImageData(inputobj)
160            self.filename = inputobj
161
162        else:
163            vedo.logger.error(f"cannot understand input type {inputtype}")
164            return
165
166        if dims is not None:
167            img.SetDimensions(dims)
168
169        if origin is not None:
170            img.SetOrigin(origin)
171
172        if spacing is not None:
173            img.SetSpacing(spacing)
174
175        self.dataset = img
176        self.transform = None
177
178        #####################################
179        mapper = vtki.new("SmartVolumeMapper")
180        mapper.SetInputData(img)
181        self.actor.SetMapper(mapper)
182
183        if img.GetPointData().GetScalars():
184            if img.GetPointData().GetScalars().GetNumberOfComponents() == 1:
185                self.properties.SetShade(True)
186                self.properties.SetInterpolationType(1)
187                self.cmap("RdBu_r")
188                self.alpha([0.0, 0.0, 0.2, 0.4, 0.8, 1.0])
189                self.alpha_gradient(None)
190                self.properties.SetScalarOpacityUnitDistance(1.0)
191
192        self.pipeline = utils.OperationNode(
193            "Volume", comment=f"dims={tuple(self.dimensions())}", c="#4cc9f0"
194        )
195        #######################################################################

This class can be initialized with a numpy object, a vtkImageData or a list of 2D bmp files.

Arguments:
  • origin : (list) set volume origin coordinates
  • spacing : (list) voxel dimensions in x, y and z.
  • dims : (list) specify the dimensions of the volume.
Example:
from vedo import Volume
vol = Volume("path/to/mydata/rec*.bmp")
vol.show()
Examples:

if a list of values is used for alphas this is interpreted as a transfer function along the range of the scalar.

mapper
197    @property
198    def mapper(self):
199        """Return the underlying `vtkMapper` object."""
200        return self.actor.GetMapper()

Return the underlying vtkMapper object.

def c(self, *args, **kwargs) -> Self:
230    def c(self, *args, **kwargs) -> Self:
231        """Deprecated. Use `Volume.cmap()` instead."""
232        vedo.logger.warning("Volume.c() is deprecated, use Volume.cmap() instead")
233        return self.cmap(*args, **kwargs)

Deprecated. Use Volume.cmap() instead.

def copy(self, deep=True) -> Volume:
372    def copy(self, deep=True) -> "Volume":
373        """Return a copy of the Volume. Alias of `clone()`."""
374        return self.clone(deep=deep)

Return a copy of the Volume. Alias of clone().

def clone(self, deep=True) -> Volume:
376    def clone(self, deep=True) -> "Volume":
377        """Return a clone copy of the Volume. Alias of `copy()`."""
378        if deep:
379            newimg = vtki.vtkImageData()
380            newimg.CopyStructure(self.dataset)
381            newimg.CopyAttributes(self.dataset)
382            newvol = Volume(newimg)
383        else:
384            newvol = Volume(self.dataset)
385
386        prop = vtki.vtkVolumeProperty()
387        prop.DeepCopy(self.properties)
388        newvol.actor.SetProperty(prop)
389        newvol.properties = prop
390
391        newvol.pipeline = utils.OperationNode("clone", parents=[self], c="#bbd0ff", shape="diamond")
392        return newvol

Return a clone copy of the Volume. Alias of copy().

def astype(self, dtype: Union[str, int]) -> Self:
394    def astype(self, dtype: Union[str, int]) -> Self:
395        """
396        Reset the type of the scalars array.
397
398        Arguments:
399            dtype : (str)
400                the type of the scalars array in
401                ["int8", "uint8", "int16", "uint16", "int32", "uint32", "float32", "float64"]
402        """
403        if dtype in ["int8", "uint8", "int16", "uint16", "int32", "uint32", "float32", "float64"]:
404            caster = vtki.new("ImageCast")
405            caster.SetInputData(self.dataset)
406            caster.SetOutputScalarType(int(vtki.array_types[dtype]))
407            caster.ClampOverflowOn()
408            caster.Update()
409            self._update(caster.GetOutput())
410            self.pipeline = utils.OperationNode(f"astype({dtype})", parents=[self], c="#4cc9f0")
411        else:
412            vedo.logger.error(f"astype(): unknown type {dtype}")
413            raise ValueError()
414        return self

Reset the type of the scalars array.

Arguments:
  • dtype : (str) the type of the scalars array in ["int8", "uint8", "int16", "uint16", "int32", "uint32", "float32", "float64"]
def component_weight(self, i: int, weight: float) -> Self:
416    def component_weight(self, i: int, weight: float) -> Self:
417        """Set the scalar component weight in range [0,1]."""
418        self.properties.SetComponentWeight(i, weight)
419        return self

Set the scalar component weight in range [0,1].

def xslice(self, i: int) -> vedo.mesh.Mesh:
421    def xslice(self, i: int) -> Mesh:
422        """Extract the slice at index `i` of volume along x-axis."""
423        vslice = vtki.new("ImageDataGeometryFilter")
424        vslice.SetInputData(self.dataset)
425        nx, ny, nz = self.dataset.GetDimensions()
426        if i > nx - 1:
427            i = nx - 1
428        vslice.SetExtent(i, i, 0, ny, 0, nz)
429        vslice.Update()
430        m = Mesh(vslice.GetOutput())
431        m.pipeline = utils.OperationNode(f"xslice {i}", parents=[self], c="#4cc9f0:#e9c46a")
432        return m

Extract the slice at index i of volume along x-axis.

def yslice(self, j: int) -> vedo.mesh.Mesh:
434    def yslice(self, j: int) -> Mesh:
435        """Extract the slice at index `j` of volume along y-axis."""
436        vslice = vtki.new("ImageDataGeometryFilter")
437        vslice.SetInputData(self.dataset)
438        nx, ny, nz = self.dataset.GetDimensions()
439        if j > ny - 1:
440            j = ny - 1
441        vslice.SetExtent(0, nx, j, j, 0, nz)
442        vslice.Update()
443        m = Mesh(vslice.GetOutput())
444        m.pipeline = utils.OperationNode(f"yslice {j}", parents=[self], c="#4cc9f0:#e9c46a")
445        return m

Extract the slice at index j of volume along y-axis.

def zslice(self, k: int) -> vedo.mesh.Mesh:
447    def zslice(self, k: int) -> Mesh:
448        """Extract the slice at index `i` of volume along z-axis."""
449        vslice = vtki.new("ImageDataGeometryFilter")
450        vslice.SetInputData(self.dataset)
451        nx, ny, nz = self.dataset.GetDimensions()
452        if k > nz - 1:
453            k = nz - 1
454        vslice.SetExtent(0, nx, 0, ny, k, k)
455        vslice.Update()
456        m = Mesh(vslice.GetOutput())
457        m.pipeline = utils.OperationNode(f"zslice {k}", parents=[self], c="#4cc9f0:#e9c46a")
458        return m

Extract the slice at index i of volume along z-axis.

def slice_plane( self, origin: List[float], normal: List[float], autocrop=False, border=0.5, mode='linear') -> vedo.mesh.Mesh:
460    def slice_plane(self, origin: List[float], normal: List[float], autocrop=False, border=0.5, mode="linear") -> Mesh:
461        """
462        Extract the slice along a given plane position and normal.
463
464        Two metadata arrays are added to the output Mesh:
465            - "shape" : contains the shape of the slice
466            - "original_bounds" : contains the original bounds of the slice
467        One can access them with e.g. `myslice.metadata["shape"]`.
468
469        Arguments:
470            origin : (list)
471                position of the plane
472            normal : (list)
473                normal to the plane
474            autocrop : (bool)
475                crop the output to the minimal possible size
476            border : (float)
477                add a border to the output slice
478            mode : (str)
479                interpolation mode, one of the following: "linear", "nearest", "cubic"
480
481        Example:
482            - [slice_plane1.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/slice_plane1.py)
483
484                ![](https://vedo.embl.es/images/volumetric/slicePlane1.gif)
485            
486            - [slice_plane2.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/slice_plane2.py)
487                
488                ![](https://vedo.embl.es/images/volumetric/slicePlane2.png)
489
490            - [slice_plane3.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/slice_plane3.py)
491
492                ![](https://vedo.embl.es/images/volumetric/slicePlane3.jpg)
493        """
494        newaxis = utils.versor(normal)
495        pos = np.array(origin)
496        initaxis = (0, 0, 1)
497        crossvec = np.cross(initaxis, newaxis)
498        angle = np.arccos(np.dot(initaxis, newaxis))
499        T = vtki.vtkTransform()
500        T.PostMultiply()
501        T.RotateWXYZ(np.rad2deg(angle), crossvec.tolist())
502        T.Translate(pos.tolist())
503
504        reslice = vtki.new("ImageReslice")
505        reslice.SetResliceAxes(T.GetMatrix())
506        reslice.SetInputData(self.dataset)
507        reslice.SetOutputDimensionality(2)
508        reslice.SetTransformInputSampling(True)
509        reslice.SetGenerateStencilOutput(False)
510        if border:
511            reslice.SetBorder(True)
512            reslice.SetBorderThickness(border)
513        else:
514            reslice.SetBorder(False)
515        if mode == "linear":
516            reslice.SetInterpolationModeToLinear()
517        elif mode == "nearest":
518            reslice.SetInterpolationModeToNearestNeighbor()
519        elif mode == "cubic":
520            reslice.SetInterpolationModeToCubic()
521        else:
522            vedo.logger.error(f"in slice_plane(): unknown interpolation mode {mode}")
523            raise ValueError()
524        reslice.SetAutoCropOutput(not autocrop)
525        reslice.Update()
526        img = reslice.GetOutput()
527
528        vslice = vtki.new("ImageDataGeometryFilter")
529        vslice.SetInputData(img)
530        vslice.Update()
531
532        msh = Mesh(vslice.GetOutput()).apply_transform(T)
533        msh.properties.LightingOff()
534
535        d0, d1, _ = img.GetDimensions()
536        varr1 = utils.numpy2vtk([d1, d0], name="shape")
537        varr2 = utils.numpy2vtk(img.GetBounds(), name="original_bounds")
538        msh.dataset.GetFieldData().AddArray(varr1)
539        msh.dataset.GetFieldData().AddArray(varr2)
540        msh.pipeline = utils.OperationNode("slice_plane", parents=[self], c="#4cc9f0:#e9c46a")
541        return msh

Extract the slice along a given plane position and normal.

Two metadata arrays are added to the output Mesh:
  • "shape" : contains the shape of the slice
  • "original_bounds" : contains the original bounds of the slice

One can access them with e.g. myslice.metadata["shape"].

Arguments:
  • origin : (list) position of the plane
  • normal : (list) normal to the plane
  • autocrop : (bool) crop the output to the minimal possible size
  • border : (float) add a border to the output slice
  • mode : (str) interpolation mode, one of the following: "linear", "nearest", "cubic"
Example:
def slab(self, slice_range=(), axis='z', operation='mean') -> vedo.mesh.Mesh:
543    def slab(self, slice_range=(), axis='z', operation="mean") -> Mesh:
544        """
545        Extract a slab from a `Volume` by combining 
546        all of the slices of an image to create a single slice.
547
548        Returns a `Mesh` containing metadata which
549        can be accessed with e.g. `mesh.metadata["slab_range"]`.
550
551        Metadata:
552            slab_range : (list)
553                contains the range of slices extracted
554            slab_axis : (str)
555                contains the axis along which the slab was extracted
556            slab_operation : (str)
557                contains the operation performed on the slab
558            slab_bounding_box : (list)
559                contains the bounding box of the slab
560
561        Arguments:
562            slice_range : (list)
563                range of slices to extract
564            axis : (str)
565                axis along which to extract the slab
566            operation : (str)
567                operation to perform on the slab,
568                allowed values are: "sum", "min", "max", "mean".
569        
570        Example:
571            - [slab.py](https://github.com/marcomusy/vedo/blob/master/examples/volumetric/slab_vol.py)
572
573            ![](https://vedo.embl.es/images/volumetric/slab_vol.jpg)
574        """
575        if len(slice_range) != 2:
576            vedo.logger.error("in slab(): slice_range is empty or invalid")
577            raise ValueError()
578        
579        islab = vtki.new("ImageSlab")
580        islab.SetInputData(self.dataset)
581
582        if operation in ["+", "add", "sum"]:
583            islab.SetOperationToSum()
584        elif "min" in operation:
585            islab.SetOperationToMin()
586        elif "max" in operation:
587            islab.SetOperationToMax()
588        elif "mean" in operation:
589            islab.SetOperationToMean()
590        else:
591            vedo.logger.error(f"in slab(): unknown operation {operation}")
592            raise ValueError()
593
594        dims = self.dimensions()
595        if axis == 'x':
596            islab.SetOrientationToX()
597            if slice_range[0]  > dims[0]-1:
598                slice_range[0] = dims[0]-1
599            if slice_range[1]  > dims[0]-1:
600                slice_range[1] = dims[0]-1
601        elif axis == 'y':
602            islab.SetOrientationToY()
603            if slice_range[0]  > dims[1]-1:
604                slice_range[0] = dims[1]-1
605            if slice_range[1]  > dims[1]-1:
606                slice_range[1] = dims[1]-1
607        elif axis == 'z':
608            islab.SetOrientationToZ()
609            if slice_range[0]  > dims[2]-1:
610                slice_range[0] = dims[2]-1
611            if slice_range[1]  > dims[2]-1:
612                slice_range[1] = dims[2]-1
613        else:
614            vedo.logger.error(f"Error in slab(): unknown axis {axis}")
615            raise RuntimeError()
616        
617        islab.SetSliceRange(slice_range)
618        islab.Update()
619
620        msh = Mesh(islab.GetOutput()).lighting('off')
621        msh.mapper.SetLookupTable(utils.ctf2lut(self, msh))
622        msh.mapper.SetScalarRange(self.scalar_range())
623
624        msh.metadata["slab_range"] = slice_range
625        msh.metadata["slab_axis"]  = axis
626        msh.metadata["slab_operation"] = operation
627
628        # compute bounds of slab
629        origin = list(self.origin())
630        spacing = list(self.spacing())
631        if axis == 'x':
632            msh.metadata["slab_bounding_box"] = [
633                origin[0] + slice_range[0]*spacing[0],
634                origin[0] + slice_range[1]*spacing[0],
635                origin[1],
636                origin[1] + dims[1]*spacing[1],
637                origin[2],
638                origin[2] + dims[2]*spacing[2],
639            ]
640        elif axis == 'y':
641            msh.metadata["slab_bounding_box"] = [
642                origin[0],
643                origin[0] + dims[0]*spacing[0],
644                origin[1] + slice_range[0]*spacing[1],
645                origin[1] + slice_range[1]*spacing[1],
646                origin[2],
647                origin[2] + dims[2]*spacing[2],
648            ]
649        elif axis == 'z':
650            msh.metadata["slab_bounding_box"] = [
651                origin[0],
652                origin[0] + dims[0]*spacing[0],
653                origin[1],
654                origin[1] + dims[1]*spacing[1],
655                origin[2] + slice_range[0]*spacing[2],
656                origin[2] + slice_range[1]*spacing[2],
657            ]
658
659        msh.pipeline = utils.OperationNode(
660            f"slab{slice_range}", 
661            comment=f"axis={axis}, operation={operation}",
662            parents=[self],
663            c="#4cc9f0:#e9c46a",
664        )
665        msh.name = "SlabMesh"
666        return msh

Extract a slab from a Volume by combining all of the slices of an image to create a single slice.

Returns a Mesh containing metadata which can be accessed with e.g. mesh.metadata["slab_range"].

Metadata:

slab_range : (list) contains the range of slices extracted slab_axis : (str) contains the axis along which the slab was extracted slab_operation : (str) contains the operation performed on the slab slab_bounding_box : (list) contains the bounding box of the slab

Arguments:
  • slice_range : (list) range of slices to extract
  • axis : (str) axis along which to extract the slab
  • operation : (str) operation to perform on the slab, allowed values are: "sum", "min", "max", "mean".
Example:

def warp( self, source: Union[vedo.pointcloud.Points, List], target: Union[vedo.pointcloud.Points, List], sigma=1, mode='3d', fit=True) -> Self:
669    def warp(
670            self,
671            source: Union["vedo.Points", List],
672            target: Union["vedo.Points", List],
673            sigma=1, mode="3d", fit=True,
674        ) -> Self:
675        """
676        Warp volume scalars within a Volume by specifying
677        source and target sets of points.
678
679        Arguments:
680            source : (Points, list)
681                the list of source points
682            target : (Points, list)
683                the list of target points
684            fit : (bool)
685                fit/adapt the old bounding box to the warped geometry
686        """
687        if isinstance(source, vedo.Points):
688            source = source.vertices
689        if isinstance(target, vedo.Points):
690            target = target.vertices
691
692        NLT = transformations.NonLinearTransform()
693        NLT.source_points = source
694        NLT.target_points = target
695        NLT.sigma = sigma
696        NLT.mode = mode
697
698        self.apply_transform(NLT, fit=fit)
699        self.pipeline = utils.OperationNode("warp", parents=[self], c="#4cc9f0")
700        return self

Warp volume scalars within a Volume by specifying source and target sets of points.

Arguments:
  • source : (Points, list) the list of source points
  • target : (Points, list) the list of target points
  • fit : (bool) fit/adapt the old bounding box to the warped geometry
def apply_transform( self, T: Union[vedo.transformations.LinearTransform, vedo.transformations.NonLinearTransform], fit=True, interpolation='cubic') -> Self:
702    def apply_transform(
703            self, 
704            T: Union[transformations.LinearTransform, transformations.NonLinearTransform],
705            fit=True, interpolation="cubic",
706        ) -> Self:
707        """
708        Apply a transform to the scalars in the volume.
709
710        Arguments:
711            T : (LinearTransform, NonLinearTransform)
712                The transformation to be applied
713            fit : (bool)
714                fit/adapt the old bounding box to the modified geometry
715            interpolation : (str)
716                one of the following: "nearest", "linear", "cubic"
717        """
718        if utils.is_sequence(T):
719            T = transformations.LinearTransform(T)
720
721        TI = T.compute_inverse()
722
723        reslice = vtki.new("ImageReslice")
724        reslice.SetInputData(self.dataset)
725        reslice.SetResliceTransform(TI.T)
726        reslice.SetOutputDimensionality(3)
727        if "lin" in interpolation.lower():
728            reslice.SetInterpolationModeToLinear()
729        elif "near" in interpolation.lower():
730            reslice.SetInterpolationModeToNearestNeighbor()
731        elif "cubic" in interpolation.lower():
732            reslice.SetInterpolationModeToCubic()
733        else:
734            vedo.logger.error(
735                f"in apply_transform: unknown interpolation mode {interpolation}")
736            raise ValueError()
737        reslice.SetAutoCropOutput(fit)
738        reslice.Update()
739        self._update(reslice.GetOutput())
740        self.transform = T
741        self.pipeline = utils.OperationNode(
742            "apply_transform", parents=[self], c="#4cc9f0")
743        return self

Apply a transform to the scalars in the volume.

Arguments:
  • T : (LinearTransform, NonLinearTransform) The transformation to be applied
  • fit : (bool) fit/adapt the old bounding box to the modified geometry
  • interpolation : (str) one of the following: "nearest", "linear", "cubic"
def imagedata(self) -> vtkmodules.vtkCommonDataModel.vtkImageData:
745    def imagedata(self) -> vtki.vtkImageData:
746        """
747        DEPRECATED:
748        Use `Volume.dataset` instead.
749
750        Return the underlying `vtkImagaData` object.
751        """
752        print("Volume.imagedata() is deprecated, use Volume.dataset instead")
753        return self.dataset

DEPRECATED: Use Volume.dataset instead.

Return the underlying vtkImagaData object.

def modified(self) -> Self:
755    def modified(self) -> Self:
756        """
757        Mark the object as modified.
758
759        Example:
760
761        - [numpy2volume0.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/numpy2volume0.py)
762        """
763        scals = self.dataset.GetPointData().GetScalars()
764        if scals:
765            scals.Modified()
766        return self

Mark the object as modified.

Example:

def tonumpy(self) -> numpy.ndarray:
768    def tonumpy(self) -> np.ndarray:
769        """
770        Get read-write access to voxels of a Volume object as a numpy array.
771
772        When you set values in the output image, you don't want numpy to reallocate the array
773        but instead set values in the existing array, so use the [:] operator.
774
775        Example:
776            `arr[:] = arr*2 + 15`
777
778        If the array is modified add a call to:
779        `volume.modified()`
780        when all your modifications are completed.
781        """
782        narray_shape = tuple(reversed(self.dataset.GetDimensions()))
783
784        scals = self.dataset.GetPointData().GetScalars()
785        comps = scals.GetNumberOfComponents()
786        if comps == 1:
787            narray = utils.vtk2numpy(scals).reshape(narray_shape)
788            narray = np.transpose(narray, axes=[2, 1, 0])
789        else:
790            narray = utils.vtk2numpy(scals).reshape(*narray_shape, comps)
791            narray = np.transpose(narray, axes=[2, 1, 0, 3])
792
793        # narray = utils.vtk2numpy(self.dataset.GetPointData().GetScalars()).reshape(narray_shape)
794        # narray = np.transpose(narray, axes=[2, 1, 0])
795        return narray

Get read-write access to voxels of a Volume object as a numpy array.

When you set values in the output image, you don't want numpy to reallocate the array but instead set values in the existing array, so use the [:] operator.

Example:

arr[:] = arr*2 + 15

If the array is modified add a call to: volume.modified() when all your modifications are completed.

shape: numpy.ndarray
797    @property
798    def shape(self) -> np.ndarray:
799        """Return the nr. of voxels in the 3 dimensions."""
800        return np.array(self.dataset.GetDimensions())

Return the nr. of voxels in the 3 dimensions.

def dimensions(self) -> numpy.ndarray:
802    def dimensions(self) -> np.ndarray:
803        """Return the nr. of voxels in the 3 dimensions."""
804        return np.array(self.dataset.GetDimensions())

Return the nr. of voxels in the 3 dimensions.

def scalar_range(self) -> numpy.ndarray:
806    def scalar_range(self) -> np.ndarray:
807        """Return the range of the scalar values."""
808        return np.array(self.dataset.GetScalarRange())

Return the range of the scalar values.

def spacing(self, s=None) -> Union[Self, Iterable[float]]:
810    def spacing(self, s=None) -> Union[Self, Iterable[float]]:
811        """Set/get the voxels size in the 3 dimensions."""
812        if s is not None:
813            self.dataset.SetSpacing(s)
814            return self
815        return np.array(self.dataset.GetSpacing())

Set/get the voxels size in the 3 dimensions.

def origin(self, s=None) -> Union[Self, Iterable[float]]:
817    def origin(self, s=None) -> Union[Self, Iterable[float]]:
818        """
819        Set/get the origin of the volumetric dataset.
820
821        The origin is the position in world coordinates of the point index (0,0,0).
822        This point does not have to be part of the dataset, in other words,
823        the dataset extent does not have to start at (0,0,0) and the origin 
824        can be outside of the dataset bounding box. 
825        The origin plus spacing determine the position in space of the points.
826        """
827        if s is not None:
828            self.dataset.SetOrigin(s)
829            return self
830        return np.array(self.dataset.GetOrigin())

Set/get the origin of the volumetric dataset.

The origin is the position in world coordinates of the point index (0,0,0). This point does not have to be part of the dataset, in other words, the dataset extent does not have to start at (0,0,0) and the origin can be outside of the dataset bounding box. The origin plus spacing determine the position in space of the points.

def pos(self, p=None) -> Union[Self, Iterable[float]]:
832    def pos(self, p=None) -> Union[Self, Iterable[float]]:
833        """Set/get the position of the volumetric dataset."""
834        if p is not None:
835            self.origin(p)
836            return self
837        return self.origin()

Set/get the position of the volumetric dataset.

def center(self) -> numpy.ndarray:
839    def center(self) -> np.ndarray:
840        """Get the center of the volumetric dataset."""
841        # note that this does not have the set method like origin and spacing
842        return np.array(self.dataset.GetCenter())

Get the center of the volumetric dataset.

def shift(self, s: list) -> Self:
844    def shift(self, s: list) -> Self:
845        """Shift the volumetric dataset by a vector."""
846        self.origin(self.origin() + np.array(s))
847        return self

Shift the volumetric dataset by a vector.

def rotate_x(self, angle: float, rad=False, around=None) -> Self:
849    def rotate_x(self, angle: float, rad=False, around=None) -> Self:
850        """
851        Rotate around x-axis. If angle is in radians set `rad=True`.
852
853        Use `around` to define a pivoting point.
854        """
855        if angle == 0:
856            return self
857        LT = transformations.LinearTransform().rotate_x(angle, rad, around)
858        return self.apply_transform(LT, fit=True, interpolation="linear")

Rotate around x-axis. If angle is in radians set rad=True.

Use around to define a pivoting point.

def rotate_y(self, angle: float, rad=False, around=None) -> Self:
860    def rotate_y(self, angle: float, rad=False, around=None) -> Self:
861        """
862        Rotate around y-axis. If angle is in radians set `rad=True`.
863
864        Use `around` to define a pivoting point.
865        """
866        if angle == 0:
867            return self
868        LT = transformations.LinearTransform().rotate_y(angle, rad, around)
869        return self.apply_transform(LT, fit=True, interpolation="linear")

Rotate around y-axis. If angle is in radians set rad=True.

Use around to define a pivoting point.

def rotate_z(self, angle: float, rad=False, around=None) -> Self:
871    def rotate_z(self, angle: float, rad=False, around=None) -> Self:
872        """
873        Rotate around z-axis. If angle is in radians set `rad=True`.
874
875        Use `around` to define a pivoting point.
876        """
877        if angle == 0:
878            return self
879        LT = transformations.LinearTransform().rotate_z(angle, rad, around)
880        return self.apply_transform(LT, fit=True, interpolation="linear")

Rotate around z-axis. If angle is in radians set rad=True.

Use around to define a pivoting point.

def get_cell_from_ijk(self, ijk: list) -> int:
882    def get_cell_from_ijk(self, ijk: list) -> int:
883        """
884        Get the voxel id number at the given ijk coordinates.
885
886        Arguments:
887            ijk : (list)
888                the ijk coordinates of the voxel
889        """
890        return self.dataset.ComputeCellId(ijk)

Get the voxel id number at the given ijk coordinates.

Arguments:
  • ijk : (list) the ijk coordinates of the voxel
def get_point_from_ijk(self, ijk: list) -> int:
892    def get_point_from_ijk(self, ijk: list) -> int:
893        """
894        Get the point id number at the given ijk coordinates.
895
896        Arguments:
897            ijk : (list)
898                the ijk coordinates of the voxel
899        """
900        return self.dataset.ComputePointId(ijk)

Get the point id number at the given ijk coordinates.

Arguments:
  • ijk : (list) the ijk coordinates of the voxel
def permute_axes(self, x: int, y: int, z: int) -> Self:
902    def permute_axes(self, x: int, y: int, z: int) -> Self:
903        """
904        Reorder the axes of the Volume by specifying
905        the input axes which are supposed to become the new X, Y, and Z.
906        """
907        imp = vtki.new("ImagePermute")
908        imp.SetFilteredAxes(x, y, z)
909        imp.SetInputData(self.dataset)
910        imp.Update()
911        self._update(imp.GetOutput())
912        self.pipeline = utils.OperationNode(
913            f"permute_axes({(x,y,z)})", parents=[self], c="#4cc9f0"
914        )
915        return self

Reorder the axes of the Volume by specifying the input axes which are supposed to become the new X, Y, and Z.

def resample(self, new_spacing: List[float], interpolation=1) -> Self:
917    def resample(self, new_spacing: List[float], interpolation=1) -> Self:
918        """
919        Resamples a `Volume` to be larger or smaller.
920
921        This method modifies the spacing of the input.
922        Linear interpolation is used to resample the data.
923
924        Arguments:
925            new_spacing : (list)
926                a list of 3 new spacings for the 3 axes
927            interpolation : (int)
928                0=nearest_neighbor, 1=linear, 2=cubic
929        """
930        rsp = vtki.new("ImageResample")
931        oldsp = self.spacing()
932        for i in range(3):
933            if oldsp[i] != new_spacing[i]:
934                rsp.SetAxisOutputSpacing(i, new_spacing[i])
935        rsp.InterpolateOn()
936        rsp.SetInterpolationMode(interpolation)
937        rsp.OptimizationOn()
938        rsp.Update()
939        self._update(rsp.GetOutput())
940        self.pipeline = utils.OperationNode(
941            "resample", comment=f"spacing: {tuple(new_spacing)}", parents=[self], c</