vedo.volume

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