vedo.volume

   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