vedo.volume

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