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