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