vedo.volume

Work with volumetric datasets (voxel data).

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

Class to describe dataset that are defined on "voxels", the 3D equivalent of 2D pixels.

Volume(inputobj=None, dims=None, origin=None, spacing=None)
 35    def __init__(
 36        self,
 37        inputobj=None,
 38        dims=None,
 39        origin=None,
 40        spacing=None,
 41    ):
 42        """
 43        This class can be initialized with a numpy object,
 44        a `vtkImageData` or a list of 2D bmp files.
 45
 46        Arguments:
 47            origin : (list)
 48                set volume origin coordinates
 49            spacing : (list)
 50                voxel dimensions in x, y and z.
 51            dims : (list)
 52                specify the dimensions of the volume.
 53
 54        Example:
 55            ```python
 56            from vedo import Volume
 57            vol = Volume("path/to/mydata/rec*.bmp")
 58            vol.show()
 59            ```
 60
 61        Examples:
 62            - [numpy2volume1.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/numpy2volume1.py)
 63
 64                ![](https://vedo.embl.es/images/volumetric/numpy2volume1.png)
 65
 66            - [read_volume2.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/read_volume2.py)
 67
 68                ![](https://vedo.embl.es/images/volumetric/read_volume2.png)
 69
 70        .. note::
 71            if a `list` of values is used for `alphas` this is interpreted
 72            as a transfer function along the range of the scalar.
 73        """
 74        self.name = "Volume"
 75        self.filename = ""
 76        self.info = {}
 77        self.time = 0
 78        self.rendered_at = set()
 79
 80        self.actor = vtk.vtkVolume()
 81        self.actor.retrieve_object = weak_ref_to(self)
 82        self.properties = self.actor.GetProperty()
 83
 84        ###################
 85        if isinstance(inputobj, str):
 86            if "https://" in inputobj:
 87                inputobj = vedo.file_io.download(inputobj, verbose=False)  # fpath
 88            elif os.path.isfile(inputobj):
 89                self.filename = inputobj
 90            else:
 91                inputobj = sorted(glob.glob(inputobj))
 92
 93        ###################
 94        inputtype = str(type(inputobj))
 95
 96        # print('Volume inputtype', inputtype, c='b')
 97
 98        if inputobj is None:
 99            img = vtk.vtkImageData()
100
101        elif utils.is_sequence(inputobj):
102
103            if isinstance(inputobj[0], str) and ".bmp" in inputobj[0].lower():
104                # scan sequence of BMP files
105                ima = vtk.new("ImageAppend")
106                ima.SetAppendAxis(2)
107                pb = utils.ProgressBar(0, len(inputobj))
108                for i in pb.range():
109                    f = inputobj[i]
110                    if "_rec_spr" in f: # OPT specific
111                        continue
112                    picr = vtk.new("BMPReader")
113                    picr.SetFileName(f)
114                    picr.Update()
115                    mgf = vtk.new("ImageMagnitude")
116                    mgf.SetInputData(picr.GetOutput())
117                    mgf.Update()
118                    ima.AddInputData(mgf.GetOutput())
119                    pb.print("loading...")
120                ima.Update()
121                img = ima.GetOutput()
122
123            else:
124
125                if len(inputobj.shape) == 1:
126                    varr = utils.numpy2vtk(inputobj)
127                else:
128                    varr = utils.numpy2vtk(inputobj.ravel(order="F"))
129                varr.SetName("input_scalars")
130
131                img = vtk.vtkImageData()
132                if dims is not None:
133                    img.SetDimensions(dims[2], dims[1], dims[0])
134                else:
135                    if len(inputobj.shape) == 1:
136                        vedo.logger.error("must set dimensions (dims keyword) in Volume")
137                        raise RuntimeError()
138                    img.SetDimensions(inputobj.shape)
139                img.GetPointData().AddArray(varr)
140                img.GetPointData().SetActiveScalars(varr.GetName())
141
142        elif isinstance(inputobj, vtk.vtkImageData):
143            img = inputobj
144
145        elif isinstance(inputobj, str):
146            if "https://" in inputobj:
147                inputobj = vedo.file_io.download(inputobj, verbose=False)
148            img = vedo.file_io.loadImageData(inputobj)
149            self.filename = inputobj
150
151        else:
152            vedo.logger.error(f"cannot understand input type {inputtype}")
153            return
154
155        if dims is not None:
156            img.SetDimensions(dims)
157
158        if origin is not None:
159            img.SetOrigin(origin)
160
161        if spacing is not None:
162            img.SetSpacing(spacing)
163
164        self.dataset = img
165        self.transform = None
166
167        #####################################
168        mapper = vtk.new("SmartVolumeMapper")
169        mapper.SetInputData(img)
170        self.actor.SetMapper(mapper)
171
172        if img.GetPointData().GetScalars():
173            if img.GetPointData().GetScalars().GetNumberOfComponents() == 1:
174                self.properties.SetShade(True)
175                self.properties.SetInterpolationType(1)
176                self.cmap("RdBu_r")
177                self.alpha([0.0, 0.0, 0.2, 0.4, 0.8, 1.0])
178                self.alpha_gradient(None)
179                self.properties.SetScalarOpacityUnitDistance(1.0)
180
181        self.pipeline = utils.OperationNode(
182            "Volume", comment=f"dims={tuple(self.dimensions())}", c="#4cc9f0"
183        )
184        #######################################################################

This class can be initialized with a numpy object, a vtkImageData or a list of 2D bmp files.

Arguments:
  • origin : (list) set volume origin coordinates
  • spacing : (list) voxel dimensions in x, y and z.
  • dims : (list) specify the dimensions of the volume.
Example:
from vedo import Volume
vol = Volume("path/to/mydata/rec*.bmp")
vol.show()
Examples:

if a list of values is used for alphas this is interpreted as a transfer function along the range of the scalar.

mapper

Return the underlying vtkMapper object.

def c(self, *args, **kwargs):
218    def c(self, *args, **kwargs):
219        """Deprecated. Use `Volume.cmap()` instead."""
220        vedo.logger.warning("Volume.c() is deprecated, use Volume.cmap() instead")
221        return self.cmap(*args, **kwargs)

Deprecated. Use Volume.cmap() instead.

def copy(self, deep=True):
359    def copy(self, deep=True):
360        """Return a copy of the Volume. Alias of `clone()`."""
361        return self.clone(deep=deep)

Return a copy of the Volume. Alias of clone().

def clone(self, deep=True):
363    def clone(self, deep=True):
364        """Return a clone copy of the Volume. Alias of `copy()`."""
365        if deep:
366            newimg = vtk.vtkImageData()
367            newimg.CopyStructure(self.dataset)
368            newimg.CopyAttributes(self.dataset)
369            newvol = Volume(newimg)
370        else:
371            newvol = Volume(self.dataset)
372
373        prop = vtk.vtkVolumeProperty()
374        prop.DeepCopy(self.properties)
375        newvol.actor.SetProperty(prop)
376        newvol.properties = prop
377
378        newvol.pipeline = utils.OperationNode("clone", parents=[self], c="#bbd0ff", shape="diamond")
379        return newvol

Return a clone copy of the Volume. Alias of copy().

def component_weight(self, i, weight):
381    def component_weight(self, i, weight):
382        """Set the scalar component weight in range [0,1]."""
383        self.properties.SetComponentWeight(i, weight)
384        return self

Set the scalar component weight in range [0,1].

def xslice(self, i):
386    def xslice(self, i):
387        """Extract the slice at index `i` of volume along x-axis."""
388        vslice = vtk.new("ImageDataGeometryFilter")
389        vslice.SetInputData(self.dataset)
390        nx, ny, nz = self.dataset.GetDimensions()
391        if i > nx - 1:
392            i = nx - 1
393        vslice.SetExtent(i, i, 0, ny, 0, nz)
394        vslice.Update()
395        m = Mesh(vslice.GetOutput())
396        m.pipeline = utils.OperationNode(f"xslice {i}", parents=[self], c="#4cc9f0:#e9c46a")
397        return m

Extract the slice at index i of volume along x-axis.

def yslice(self, j):
399    def yslice(self, j):
400        """Extract the slice at index `j` of volume along y-axis."""
401        vslice = vtk.new("ImageDataGeometryFilter")
402        vslice.SetInputData(self.dataset)
403        nx, ny, nz = self.dataset.GetDimensions()
404        if j > ny - 1:
405            j = ny - 1
406        vslice.SetExtent(0, nx, j, j, 0, nz)
407        vslice.Update()
408        m = Mesh(vslice.GetOutput())
409        m.pipeline = utils.OperationNode(f"yslice {j}", parents=[self], c="#4cc9f0:#e9c46a")
410        return m

Extract the slice at index j of volume along y-axis.

def zslice(self, k):
412    def zslice(self, k):
413        """Extract the slice at index `i` of volume along z-axis."""
414        vslice = vtk.new("ImageDataGeometryFilter")
415        vslice.SetInputData(self.dataset)
416        nx, ny, nz = self.dataset.GetDimensions()
417        if k > nz - 1:
418            k = nz - 1
419        vslice.SetExtent(0, nx, 0, ny, k, k)
420        vslice.Update()
421        m = Mesh(vslice.GetOutput())
422        m.pipeline = utils.OperationNode(f"zslice {k}", parents=[self], c="#4cc9f0:#e9c46a")
423        return m

Extract the slice at index i of volume along z-axis.

def slice_plane(self, origin=(0, 0, 0), normal=(1, 1, 1), autocrop=False):
425    def slice_plane(self, origin=(0, 0, 0), normal=(1, 1, 1), autocrop=False):
426        """
427        Extract the slice along a given plane position and normal.
428
429        Example:
430            - [slice_plane1.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/slice_plane1.py)
431
432                ![](https://vedo.embl.es/images/volumetric/slicePlane1.gif)
433        """
434        reslice = vtk.new("ImageReslice")
435        reslice.SetInputData(self.dataset)
436        reslice.SetOutputDimensionality(2)
437        newaxis = utils.versor(normal)
438        pos = np.array(origin)
439        initaxis = (0, 0, 1)
440        crossvec = np.cross(initaxis, newaxis)
441        angle = np.arccos(np.dot(initaxis, newaxis))
442        T = vtk.vtkTransform()
443        T.PostMultiply()
444        T.RotateWXYZ(np.rad2deg(angle), crossvec)
445        T.Translate(pos)
446        M = T.GetMatrix()
447        reslice.SetResliceAxes(M)
448        reslice.SetInterpolationModeToLinear()
449        reslice.SetAutoCropOutput(not autocrop)
450        reslice.Update()
451        vslice = vtk.new("ImageDataGeometryFilter")
452        vslice.SetInputData(reslice.GetOutput())
453        vslice.Update()
454        msh = Mesh(vslice.GetOutput())
455        msh.apply_transform(T)
456        msh.pipeline = utils.OperationNode(
457            "slice_plane", parents=[self], c="#4cc9f0:#e9c46a")
458        return msh

Extract the slice along a given plane position and normal.

Example:
def slab(self, slice_range=(), axis='z', operation='mean'):
460    def slab(self, slice_range=(), axis='z', operation="mean"):
461        """
462        Extract a slab from a `Volume` by combining 
463        all of the slices of an image to create a single slice.
464
465        Returns a `Mesh` containing metadata which
466        can be accessed with e.g. `mesh.metadata["slab_range"]`.
467
468        Metadata:
469            slab_range : (list)
470                contains the range of slices extracted
471            slab_axis : (str)
472                contains the axis along which the slab was extracted
473            slab_operation : (str)
474                contains the operation performed on the slab
475            slab_bounding_box : (list)
476                contains the bounding box of the slab
477
478        Arguments:
479            slice_range : (list)
480                range of slices to extract
481            axis : (str)
482                axis along which to extract the slab
483            operation : (str)
484                operation to perform on the slab,
485                allowed values are: "sum", "min", "max", "mean".
486        """
487        if len(slice_range) != 2:
488            vedo.logger.error("in slab(): slice_range is empty or invalid")
489            raise ValueError()
490        
491        islab = vtk.new("ImageSlab")
492        islab.SetInputData(self.dataset)
493
494        if operation in ["+", "add", "sum"]:
495            islab.SetOperationToSum()
496        elif "min" in operation:
497            islab.SetOperationToMin()
498        elif "max" in operation:
499            islab.SetOperationToMax()
500        elif "mean" in operation:
501            islab.SetOperationToMean()
502        else:
503            vedo.logger.error(f"in slab(): unknown operation {operation}")
504            raise ValueError()
505
506        dims = self.dimensions()
507        if axis == 'x':
508            islab.SetOrientationToX()
509            if slice_range[0]  > dims[0]-1:
510                slice_range[0] = dims[0]-1
511            if slice_range[1]  > dims[0]-1:
512                slice_range[1] = dims[0]-1
513        elif axis == 'y':
514            islab.SetOrientationToY()
515            if slice_range[0]  > dims[1]-1:
516                slice_range[0] = dims[1]-1
517            if slice_range[1]  > dims[1]-1:
518                slice_range[1] = dims[1]-1
519        elif axis == 'z':
520            islab.SetOrientationToZ()
521            if slice_range[0]  > dims[2]-1:
522                slice_range[0] = dims[2]-1
523            if slice_range[1]  > dims[2]-1:
524                slice_range[1] = dims[2]-1
525        else:
526            vedo.logger.error(f"Error in slab(): unknown axis {axis}")
527            raise RuntimeError()
528        
529        islab.SetSliceRange(slice_range)
530        islab.Update()
531
532        msh = Mesh(islab.GetOutput()).lighting('off')
533        msh.mapper.SetLookupTable(utils.ctf2lut(self, msh))
534        msh.mapper.SetScalarRange(self.scalar_range())
535
536        msh.metadata["slab_range"] = slice_range
537        msh.metadata["slab_axis"]  = axis
538        msh.metadata["slab_operation"] = operation
539
540        # compute bounds of slab
541        origin = self.origin()
542        spacing = self.spacing()
543        if axis == 'x':
544            msh.metadata["slab_bounding_box"] = [
545                origin[0] + slice_range[0]*spacing[0],
546                origin[0] + slice_range[1]*spacing[0],
547                origin[1],
548                origin[1] + dims[1]*spacing[1],
549                origin[2],
550                origin[2] + dims[2]*spacing[2],
551            ]
552        elif axis == 'y':
553            msh.metadata["slab_bounding_box"] = [
554                origin[0],
555                origin[0] + dims[0]*spacing[0],
556                origin[1] + slice_range[0]*spacing[1],
557                origin[1] + slice_range[1]*spacing[1],
558                origin[2],
559                origin[2] + dims[2]*spacing[2],
560            ]
561        elif axis == 'z':
562            msh.metadata["slab_bounding_box"] = [
563                origin[0],
564                origin[0] + dims[0]*spacing[0],
565                origin[1],
566                origin[1] + dims[1]*spacing[1],
567                origin[2] + slice_range[0]*spacing[2],
568                origin[2] + slice_range[1]*spacing[2],
569            ]
570
571        msh.pipeline = utils.OperationNode(
572            f"slab{slice_range}", 
573            comment=f"axis={axis}, operation={operation}",
574            parents=[self],
575            c="#4cc9f0:#e9c46a",
576        )
577        msh.name = "SlabMesh"
578        return msh

Extract a slab from a Volume by combining all of the slices of an image to create a single slice.

Returns a Mesh containing metadata which can be accessed with e.g. mesh.metadata["slab_range"].

Metadata:

slab_range : (list) contains the range of slices extracted slab_axis : (str) contains the axis along which the slab was extracted slab_operation : (str) contains the operation performed on the slab slab_bounding_box : (list) contains the bounding box of the slab

Arguments:
  • slice_range : (list) range of slices to extract
  • axis : (str) axis along which to extract the slab
  • operation : (str) operation to perform on the slab, allowed values are: "sum", "min", "max", "mean".
def warp(self, source, target, sigma=1, mode='3d', fit=True):
581    def warp(self, source, target, sigma=1, mode="3d", fit=True):
582        """
583        Warp volume scalars within a Volume by specifying
584        source and target sets of points.
585
586        Arguments:
587            source : (Points, list)
588                the list of source points
589            target : (Points, list)
590                the list of target points
591            fit : (bool)
592                fit/adapt the old bounding box to the warped geometry
593        """
594        if isinstance(source, vedo.Points):
595            source = source.vertices
596        if isinstance(target, vedo.Points):
597            target = target.vertices
598
599        NLT = transformations.NonLinearTransform()
600        NLT.source_points = source
601        NLT.target_points = target
602        NLT.sigma = sigma
603        NLT.mode = mode
604
605        self.apply_transform(NLT, fit=fit)
606        self.pipeline = utils.OperationNode("warp", parents=[self], c="#4cc9f0")
607        return self

Warp volume scalars within a Volume by specifying source and target sets of points.

Arguments:
  • source : (Points, list) the list of source points
  • target : (Points, list) the list of target points
  • fit : (bool) fit/adapt the old bounding box to the warped geometry
def apply_transform(self, T, fit=True, interpolation='cubic'):
609    def apply_transform(self, T, fit=True, interpolation="cubic"):
610        """
611        Apply a transform to the scalars in the volume.
612
613        Arguments:
614            T : (LinearTransform, NonLinearTransform)
615                The transformation to be applied
616            fit : (bool)
617                fit/adapt the old bounding box to the modified geometry
618            interpolation : (str)
619                one of the following: "linear", "nearest", "cubic"
620        """
621        TI = T.compute_inverse()
622        reslice = vtk.new("ImageReslice")
623        reslice.SetInputData(self.dataset)
624        reslice.SetResliceTransform(TI.T)
625        reslice.SetOutputDimensionality(3)
626        if "lin" in interpolation.lower():
627            reslice.SetInterpolationModeToLinear()
628        elif "near" in interpolation.lower():
629            reslice.SetInterpolationModeToNearestNeighbor()
630        elif "cubic" in interpolation.lower():
631            reslice.SetInterpolationModeToCubic()
632        else:
633            vedo.logger.error(
634                f"in apply_transform: unknown interpolation mode {interpolation}")
635            raise ValueError()
636        reslice.SetAutoCropOutput(fit)
637        reslice.Update()
638        self._update(reslice.GetOutput())
639        self.transform = T
640        self.pipeline = utils.OperationNode(
641            "apply_transform", parents=[self], c="#4cc9f0")
642        return self

Apply a transform to the scalars in the volume.

Arguments:
  • T : (LinearTransform, NonLinearTransform) The transformation to be applied
  • fit : (bool) fit/adapt the old bounding box to the modified geometry
  • interpolation : (str) one of the following: "linear", "nearest", "cubic"
def imagedata(self):
644    def imagedata(self):
645        """
646        DEPRECATED:
647        Use `Volume.dataset` instead.
648
649        Return the underlying `vtkImagaData` object.
650        """
651        print("Volume.imagedata() is deprecated, use Volume.dataset instead")
652        return self.dataset

DEPRECATED: Use Volume.dataset instead.

Return the underlying vtkImagaData object.

def modified(self):
654    def modified(self):
655        """
656        Mark the object as modified.
657
658        Example:
659
660        - [numpy2volume0.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/numpy2volume0.py)
661        """
662        self.dataset.GetPointData().GetScalars().Modified()
663        return self

Mark the object as modified.

Example:

def tonumpy(self):
665    def tonumpy(self):
666        """
667        Get read-write access to voxels of a Volume object as a numpy array.
668
669        When you set values in the output image, you don't want numpy to reallocate the array
670        but instead set values in the existing array, so use the [:] operator.
671
672        Example:
673            `arr[:] = arr*2 + 15`
674
675        If the array is modified add a call to:
676        `volume.modified()`
677        when all your modifications are completed.
678        """
679        narray_shape = tuple(reversed(self.dataset.GetDimensions()))
680
681        scals = self.dataset.GetPointData().GetScalars()
682        comps = scals.GetNumberOfComponents()
683        if comps == 1:
684            narray = utils.vtk2numpy(scals).reshape(narray_shape)
685            narray = np.transpose(narray, axes=[2, 1, 0])
686        else:
687            narray = utils.vtk2numpy(scals).reshape(*narray_shape, comps)
688            narray = np.transpose(narray, axes=[2, 1, 0, 3])
689
690        # narray = utils.vtk2numpy(self.dataset.GetPointData().GetScalars()).reshape(narray_shape)
691        # narray = np.transpose(narray, axes=[2, 1, 0])
692
693        return narray

Get read-write access to voxels of a Volume object as a numpy array.

When you set values in the output image, you don't want numpy to reallocate the array but instead set values in the existing array, so use the [:] operator.

Example:

arr[:] = arr*2 + 15

If the array is modified add a call to: volume.modified() when all your modifications are completed.

shape

Return the nr. of voxels in the 3 dimensions.

def dimensions(self):
700    def dimensions(self):
701        """Return the nr. of voxels in the 3 dimensions."""
702        return np.array(self.dataset.GetDimensions())

Return the nr. of voxels in the 3 dimensions.

def scalar_range(self):
704    def scalar_range(self):
705        """Return the range of the scalar values."""
706        return np.array(self.dataset.GetScalarRange())

Return the range of the scalar values.

def spacing(self, s=None):
708    def spacing(self, s=None):
709        """Set/get the voxels size in the 3 dimensions."""
710        if s is not None:
711            self.dataset.SetSpacing(s)
712            return self
713        return np.array(self.dataset.GetSpacing())

Set/get the voxels size in the 3 dimensions.

def origin(self, s=None):
715    def origin(self, s=None):
716        """
717        Set/get the origin of the volumetric dataset.
718
719        The origin is the position in world coordinates of the point index (0,0,0).
720        This point does not have to be part of the dataset, in other words,
721        the dataset extent does not have to start at (0,0,0) and the origin 
722        can be outside of the dataset bounding box. 
723        The origin plus spacing determine the position in space of the points.
724        """
725        if s is not None:
726            self.dataset.SetOrigin(s)
727            return self
728        return np.array(self.dataset.GetOrigin())

Set/get the origin of the volumetric dataset.

The origin is the position in world coordinates of the point index (0,0,0). This point does not have to be part of the dataset, in other words, the dataset extent does not have to start at (0,0,0) and the origin can be outside of the dataset bounding box. The origin plus spacing determine the position in space of the points.

def center(self):
730    def center(self):
731        """Get the center of the volumetric dataset."""
732        # note that this does not have the set method like origin and spacing
733        return np.array(self.dataset.GetCenter())

Get the center of the volumetric dataset.

def get_cell_from_ijk(self, ijk):
735    def get_cell_from_ijk(self, ijk):
736        """
737        Get the voxel id number at the given ijk coordinates.
738
739        Arguments:
740            ijk : (list)
741                the ijk coordinates of the voxel
742        """
743        return self.ComputeCellId(ijk)

Get the voxel id number at the given ijk coordinates.

Arguments:
  • ijk : (list) the ijk coordinates of the voxel
def get_point_from_ijk(self, ijk):
745    def get_point_from_ijk(self, ijk):
746        """
747        Get the point id number at the given ijk coordinates.
748
749        Arguments:
750            ijk : (list)
751                the ijk coordinates of the voxel
752        """
753        return self.ComputePointId(ijk)

Get the point id number at the given ijk coordinates.

Arguments:
  • ijk : (list) the ijk coordinates of the voxel
def permute_axes(self, x, y, z):
755    def permute_axes(self, x, y, z):
756        """
757        Reorder the axes of the Volume by specifying
758        the input axes which are supposed to become the new X, Y, and Z.
759        """
760        imp = vtk.new("ImagePermute")
761        imp.SetFilteredAxes(x, y, z)
762        imp.SetInputData(self.dataset)
763        imp.Update()
764        self._update(imp.GetOutput())
765        self.pipeline = utils.OperationNode(
766            f"permute_axes({(x,y,z)})", parents=[self], c="#4cc9f0"
767        )
768        return self

Reorder the axes of the Volume by specifying the input axes which are supposed to become the new X, Y, and Z.

def resample(self, new_spacing, interpolation=1):
770    def resample(self, new_spacing, interpolation=1):
771        """
772        Resamples a `Volume` to be larger or smaller.
773
774        This method modifies the spacing of the input.
775        Linear interpolation is used to resample the data.
776
777        Arguments:
778            new_spacing : (list)
779                a list of 3 new spacings for the 3 axes
780            interpolation : (int)
781                0=nearest_neighbor, 1=linear, 2=cubic
782        """
783        rsp = vtk.new("ImageResample")
784        oldsp = self.spacing()
785        for i in range(3):
786            if oldsp[i] != new_spacing[i]:
787                rsp.SetAxisOutputSpacing(i, new_spacing[i])
788        rsp.InterpolateOn()
789        rsp.SetInterpolationMode(interpolation)
790        rsp.OptimizationOn()
791        rsp.Update()
792        self._update(rsp.GetOutput())
793        self.pipeline = utils.OperationNode(
794            "resample", comment=f"spacing: {tuple(new_spacing)}", parents=[self], c="#4cc9f0"
795        )
796        return self

Resamples a Volume to be larger or smaller.

This method modifies the spacing of the input. Linear interpolation is used to resample the data.

Arguments:
  • new_spacing : (list) a list of 3 new spacings for the 3 axes
  • interpolation : (int) 0=nearest_neighbor, 1=linear, 2=cubic
def threshold(self, above=None, below=None, replace=None, replace_value=None):
799    def threshold(self, above=None, below=None, replace=None, replace_value=None):
800        """
801        Binary or continuous volume thresholding.
802        Find the voxels that contain a value above/below the input values
803        and replace them with a new value (default is 0).
804        """
805        th = vtk.new("ImageThreshold")
806        th.SetInputData(self.dataset)
807
808        # sanity checks
809        if above is not None and below is not None:
810            if above == below:
811                return self
812            if above > below:
813                vedo.logger.warning("in volume.threshold(), above > below, skip.")
814                return self
815
816        ## cases
817        if below is not None and above is not None:
818            th.ThresholdBetween(above, below)
819
820        elif above is not None:
821            th.ThresholdByUpper(above)
822
823        elif below is not None:
824            th.ThresholdByLower(below)
825
826        ##
827        if replace is not None:
828            th.SetReplaceIn(True)
829            th.SetInValue(replace)
830        else:
831            th.SetReplaceIn(False)
832
833        if replace_value is not None:
834            th.SetReplaceOut(True)
835            th.SetOutValue(replace_value)
836        else:
837            th.SetReplaceOut(False)
838
839        th.Update()
840        self._update(th.GetOutput())
841        self.pipeline = utils.OperationNode("threshold", parents=[self], c="#4cc9f0")
842        return self

Binary or continuous volume thresholding. Find the voxels that contain a value above/below the input values and replace them with a new value (default is 0).

def crop( self, left=None, right=None, back=None, front=None, bottom=None, top=None, VOI=()):
844    def crop(self, left=None, right=None, back=None, front=None, bottom=None, top=None, VOI=()):
845        """
846        Crop a `Volume` object.
847
848        Arguments:
849            left : (float)
850                fraction to crop from the left plane (negative x)
851            right : (float)
852                fraction to crop from the right plane (positive x)
853            back : (float)
854                fraction to crop from the back plane (negative y)
855            front : (float)
856                fraction to crop from the front plane (positive y)
857            bottom : (float)
858                fraction to crop from the bottom plane (negative z)
859            top : (float)
860                fraction to crop from the top plane (positive z)
861            VOI : (list)
862                extract Volume Of Interest expressed in voxel numbers
863
864        Example:
865            `vol.crop(VOI=(xmin, xmax, ymin, ymax, zmin, zmax)) # all integers nrs`
866        """
867        extractVOI = vtk.new("ExtractVOI")
868        extractVOI.SetInputData(self.dataset)
869
870        if VOI:
871            extractVOI.SetVOI(VOI)
872        else:
873            d = self.dataset.GetDimensions()
874            bx0, bx1, by0, by1, bz0, bz1 = 0, d[0]-1, 0, d[1]-1, 0, d[2]-1
875            if left is not None:   bx0 = int((d[0]-1)*left)
876            if right is not None:  bx1 = int((d[0]-1)*(1-right))
877            if back is not None:   by0 = int((d[1]-1)*back)
878            if front is not None:  by1 = int((d[1]-1)*(1-front))
879            if bottom is not None: bz0 = int((d[2]-1)*bottom)
880            if top is not None:    bz1 = int((d[2]-1)*(1-top))
881            extractVOI.SetVOI(bx0, bx1, by0, by1, bz0, bz1)
882        extractVOI.Update()
883        self._update(extractVOI.GetOutput())
884
885        self.pipeline = utils.OperationNode(
886            "crop", parents=[self], c="#4cc9f0", comment=f"dims={tuple(self.dimensions())}"
887        )
888        return self

Crop a Volume object.

Arguments:
  • left : (float) fraction to crop from the left plane (negative x)
  • right : (float) fraction to crop from the right plane (positive x)
  • back : (float) fraction to crop from the back plane (negative y)
  • front : (float) fraction to crop from the front plane (positive y)
  • bottom : (float) fraction to crop from the bottom plane (negative z)
  • top : (float) fraction to crop from the top plane (positive z)
  • VOI : (list) extract Volume Of Interest expressed in voxel numbers
Example:

vol.crop(VOI=(xmin, xmax, ymin, ymax, zmin, zmax)) # all integers nrs

def append(self, volumes, axis='z', preserve_extents=False):
890    def append(self, volumes, axis="z", preserve_extents=False):
891        """
892        Take the components from multiple inputs and merges them into one output.
893        Except for the append axis, all inputs must have the same extent.
894        All inputs must have the same number of scalar components.
895        The output has the same origin and spacing as the first input.
896        The origin and spacing of all other inputs are ignored.
897        All inputs must have the same scalar type.
898
899        Arguments:
900            axis : (int, str)
901                axis expanded to hold the multiple images
902            preserve_extents : (bool)
903                if True, the extent of the inputs is used to place
904                the image in the output. The whole extent of the output is the union of the input
905                whole extents. Any portion of the output not covered by the inputs is set to zero.
906                The origin and spacing is taken from the first input.
907
908        Example:
909            ```python
910            from vedo import Volume, dataurl
911            vol = Volume(dataurl+'embryo.tif')
912            vol.append(vol, axis='x').show().close()
913            ```
914            ![](https://vedo.embl.es/images/feats/volume_append.png)
915        """
916        ima = vtk.new("ImageAppend")
917        ima.SetInputData(self.dataset)
918        if not utils.is_sequence(volumes):
919            volumes = [volumes]
920        for volume in volumes:
921            if isinstance(volume, vtk.vtkImageData):
922                ima.AddInputData(volume)
923            else:
924                ima.AddInputData(volume.dataset)
925        ima.SetPreserveExtents(preserve_extents)
926        if axis == "x":
927            axis = 0
928        elif axis == "y":
929            axis = 1
930        elif axis == "z":
931            axis = 2
932        ima.SetAppendAxis(axis)
933        ima.Update()
934        self._update(ima.GetOutput())
935
936        self.pipeline = utils.OperationNode(
937            "append",
938            parents=[self, *volumes],
939            c="#4cc9f0",
940            comment=f"dims={tuple(self.dimensions())}",
941        )
942        return self

Take the components from multiple inputs and merges them into one output. Except for the append axis, all inputs must have the same extent. All inputs must have the same number of scalar components. The output has the same origin and spacing as the first input. The origin and spacing of all other inputs are ignored. All inputs must have the same scalar type.

Arguments:
  • axis : (int, str) axis expanded to hold the multiple images
  • preserve_extents : (bool) if True, the extent of the inputs is used to place the image in the output. The whole extent of the output is the union of the input whole extents. Any portion of the output not covered by the inputs is set to zero. The origin and spacing is taken from the first input.
Example:
from vedo import Volume, dataurl
vol = Volume(dataurl+'embryo.tif')
vol.append(vol, axis='x').show().close()

def pad(self, voxels=10, value=0):
944    def pad(self, voxels=10, value=0):
945        """
946        Add the specified number of voxels at the `Volume` borders.
947        Voxels can be a list formatted as `[nx0, nx1, ny0, ny1, nz0, nz1]`.
948
949        Arguments:
950            voxels : (int, list)
951                number of voxels to be added (or a list of length 4)
952            value : (int)
953                intensity value (gray-scale color) of the padding
954
955        Example:
956            ```python
957            from vedo import Volume, dataurl, show
958            iso = Volume(dataurl+'embryo.tif').isosurface()
959            vol = iso.binarize(spacing=(100, 100, 100)).pad(10)
960            vol.dilate([15,15,15])
961            show(iso, vol.isosurface(), N=2, axes=1)
962            ```
963            ![](https://vedo.embl.es/images/volumetric/volume_pad.png)
964        """
965        x0, x1, y0, y1, z0, z1 = self.dataset.GetExtent()
966        pf = vtk.new("ImageConstantPad")
967        pf.SetInputData(self.dataset)
968        pf.SetConstant(value)
969        if utils.is_sequence(voxels):
970            pf.SetOutputWholeExtent(
971                x0 - voxels[0], x1 + voxels[1],
972                y0 - voxels[2], y1 + voxels[3],
973                z0 - voxels[4], z1 + voxels[5],
974            )
975        else:
976            pf.SetOutputWholeExtent(
977                x0 - voxels, x1 + voxels,
978                y0 - voxels, y1 + voxels,
979                z0 - voxels, z1 + voxels,
980            )
981        pf.Update()
982        self._update(pf.GetOutput())
983        self.pipeline = utils.OperationNode(
984            "pad", comment=f"{voxels} voxels", parents=[self], c="#f28482"
985        )
986        return self

Add the specified number of voxels at the Volume borders. Voxels can be a list formatted as [nx0, nx1, ny0, ny1, nz0, nz1].

Arguments:
  • voxels : (int, list) number of voxels to be added (or a list of length 4)
  • value : (int) intensity value (gray-scale color) of the padding
Example:
from vedo import Volume, dataurl, show
iso = Volume(dataurl+'embryo.tif').isosurface()
vol = iso.binarize(spacing=(100, 100, 100)).pad(10)
vol.dilate([15,15,15])
show(iso, vol.isosurface(), N=2, axes=1)

def resize(self, newdims):
 988    def resize(self, newdims):
 989        """Increase or reduce the number of voxels of a Volume with interpolation."""
 990        rsz = vtk.new("ImageResize")
 991        rsz.SetResizeMethodToOutputDimensions()
 992        rsz.SetInputData(self.dataset)
 993        rsz.SetOutputDimensions(newdims)
 994        rsz.Update()
 995        self.dataset = rsz.GetOutput()
 996        self._update(self.dataset)
 997        self.pipeline = utils.OperationNode(
 998            "resize", parents=[self], c="#4cc9f0", comment=f"dims={tuple(self.dimensions())}"
 999        )
1000        return self

Increase or reduce the number of voxels of a Volume with interpolation.

def normalize(self):
1002    def normalize(self):
1003        """Normalize that scalar components for each point."""
1004        norm = vtk.new("ImageNormalize")
1005        norm.SetInputData(self.dataset)
1006        norm.Update()
1007        self._update(norm.GetOutput())
1008        self.pipeline = utils.OperationNode("normalize", parents=[self], c="#4cc9f0")
1009        return self

Normalize that scalar components for each point.

def mirror(self, axis='x'):
1011    def mirror(self, axis="x"):
1012        """
1013        Mirror flip along one of the cartesian axes.
1014        """
1015        img = self.dataset
1016
1017        ff = vtk.new("ImageFlip")
1018        ff.SetInputData(img)
1019        if axis.lower() == "x":
1020            ff.SetFilteredAxis(0)
1021        elif axis.lower() == "y":
1022            ff.SetFilteredAxis(1)
1023        elif axis.lower() == "z":
1024            ff.SetFilteredAxis(2)
1025        else:
1026            vedo.logger.error("mirror must be set to either x, y, z or n")
1027            raise RuntimeError()
1028        ff.Update()
1029        self._update(ff.GetOutput())
1030        self.pipeline = utils.OperationNode(f"mirror {axis}", parents=[self], c="#4cc9f0")
1031        return self

Mirror flip along one of the cartesian axes.

def operation(self, operation, volume2=None):
1033    def operation(self, operation, volume2=None):
1034        """
1035        Perform operations with `Volume` objects.
1036        Keyword `volume2` can be a constant float.
1037
1038        Possible operations are:
1039        ```
1040        +, -, /, 1/x, sin, cos, exp, log,
1041        abs, **2, sqrt, min, max, atan, atan2, median,
1042        mag, dot, gradient, divergence, laplacian.
1043        ```
1044
1045        Examples:
1046            - [volume_operations.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/volume_operations.py)
1047        """
1048        op = operation.lower()
1049        image1 = self.dataset
1050
1051        mf = None
1052        if op in ["median"]:
1053            mf = vtk.new("ImageMedian3D")
1054            mf.SetInputData(image1)
1055        elif op in ["mag"]:
1056            mf = vtk.new("ImageMagnitude")
1057            mf.SetInputData(image1)
1058        elif op in ["dot", "dotproduct"]:
1059            mf = vtk.new("ImageDotProduct")
1060            mf.SetInput1Data(image1)
1061            mf.SetInput2Data(volume2.dataset)
1062        elif op in ["grad", "gradient"]:
1063            mf = vtk.new("ImageGradient")
1064            mf.SetDimensionality(3)
1065            mf.SetInputData(image1)
1066        elif op in ["div", "divergence"]:
1067            mf = vtk.new("ImageDivergence")
1068            mf.SetInputData(image1)
1069        elif op in ["laplacian"]:
1070            mf = vtk.new("ImageLaplacian")
1071            mf.SetDimensionality(3)
1072            mf.SetInputData(image1)
1073
1074        if mf is not None:
1075            mf.Update()
1076            vol = Volume(mf.GetOutput())
1077            vol.pipeline = utils.OperationNode(
1078                "operation", comment=f"{op}", parents=[self], c="#4cc9f0", shape="cylinder"
1079            )
1080            return vol  ###########################
1081
1082        mat = vtk.new("ImageMathematics")
1083        mat.SetInput1Data(image1)
1084
1085        K = None
1086
1087        if utils.is_number(volume2):
1088            K = volume2
1089            mat.SetConstantK(K)
1090            mat.SetConstantC(K)
1091
1092        elif volume2 is not None:  # assume image2 is a constant value
1093            mat.SetInput2Data(volume2.dataset)
1094
1095        # ###########################
1096        if op in ["+", "add", "plus"]:
1097            if K:
1098                mat.SetOperationToAddConstant()
1099            else:
1100                mat.SetOperationToAdd()
1101
1102        elif op in ["-", "subtract", "minus"]:
1103            if K:
1104                mat.SetConstantC(-float(K))
1105                mat.SetOperationToAddConstant()
1106            else:
1107                mat.SetOperationToSubtract()
1108
1109        elif op in ["*", "multiply", "times"]:
1110            if K:
1111                mat.SetOperationToMultiplyByK()
1112            else:
1113                mat.SetOperationToMultiply()
1114
1115        elif op in ["/", "divide"]:
1116            if K:
1117                mat.SetConstantK(1.0 / K)
1118                mat.SetOperationToMultiplyByK()
1119            else:
1120                mat.SetOperationToDivide()
1121
1122        elif op in ["1/x", "invert"]:
1123            mat.SetOperationToInvert()
1124        elif op in ["sin"]:
1125            mat.SetOperationToSin()
1126        elif op in ["cos"]:
1127            mat.SetOperationToCos()
1128        elif op in ["exp"]:
1129            mat.SetOperationToExp()
1130        elif op in ["log"]:
1131            mat.SetOperationToLog()
1132        elif op in ["abs"]:
1133            mat.SetOperationToAbsoluteValue()
1134        elif op in ["**2", "square"]:
1135            mat.SetOperationToSquare()
1136        elif op in ["sqrt", "sqr"]:
1137            mat.SetOperationToSquareRoot()
1138        elif op in ["min"]:
1139            mat.SetOperationToMin()
1140        elif op in ["max"]:
1141            mat.SetOperationToMax()
1142        elif op in ["atan"]:
1143            mat.SetOperationToATAN()
1144        elif op in ["atan2"]:
1145            mat.SetOperationToATAN2()
1146        else:
1147            vedo.logger.error(f"unknown operation {operation}")
1148            raise RuntimeError()
1149        mat.Update()
1150
1151        self._update(mat.GetOutput())
1152
1153        self.pipeline = utils.OperationNode(
1154            "operation", comment=f"{op}", parents=[self, volume2], shape="cylinder", c="#4cc9f0"
1155        )
1156        return self

Perform operations with Volume objects. Keyword volume2 can be a constant float.

Possible operations are:

+, -, /, 1/x, sin, cos, exp, log,
abs, **2, sqrt, min, max, atan, atan2, median,
mag, dot, gradient, divergence, laplacian.
Examples:
def frequency_pass_filter(self, low_cutoff=None, high_cutoff=None, order=1):
1158    def frequency_pass_filter(self, low_cutoff=None, high_cutoff=None, order=1):
1159        """
1160        Low-pass and high-pass filtering become trivial in the frequency domain.
1161        A portion of the pixels/voxels are simply masked or attenuated.
1162        This function applies a high pass Butterworth filter that attenuates the
1163        frequency domain image.
1164
1165        The gradual attenuation of the filter is important.
1166        A simple high-pass filter would simply mask a set of pixels in the frequency domain,
1167        but the abrupt transition would cause a ringing effect in the spatial domain.
1168
1169        Arguments:
1170            low_cutoff : (list)
1171                the cutoff frequencies for x, y and z
1172            high_cutoff : (list)
1173                the cutoff frequencies for x, y and z
1174            order : (int)
1175                order determines sharpness of the cutoff curve
1176        """
1177        # https://lorensen.github.io/VTKExamples/site/Cxx/ImageProcessing/IdealHighPass
1178        fft = vtk.new("ImageFFT")
1179        fft.SetInputData(self.dataset)
1180        fft.Update()
1181        out = fft.GetOutput()
1182
1183        if high_cutoff:
1184            blp = vtk.new("ImageButterworthLowPass")
1185            blp.SetInputData(out)
1186            blp.SetCutOff(high_cutoff)
1187            blp.SetOrder(order)
1188            blp.Update()
1189            out = blp.GetOutput()
1190
1191        if low_cutoff:
1192            bhp = vtk.new("ImageButterworthHighPass")
1193            bhp.SetInputData(out)
1194            bhp.SetCutOff(low_cutoff)
1195            bhp.SetOrder(order)
1196            bhp.Update()
1197            out = bhp.GetOutput()
1198
1199        rfft = vtk.new("ImageRFFT")
1200        rfft.SetInputData(out)
1201        rfft.Update()
1202
1203        ecomp = vtk.new("ImageExtractComponents")
1204        ecomp.SetInputData(rfft.GetOutput())
1205        ecomp.SetComponents(0)
1206        ecomp.Update()
1207        self._update(ecomp.GetOutput())
1208        self.pipeline = utils.OperationNode("frequency_pass_filter", parents=[self], c="#4cc9f0")
1209        return self

Low-pass and high-pass filtering become trivial in the frequency domain. A portion of the pixels/voxels are simply masked or attenuated. This function applies a high pass Butterworth filter that attenuates the frequency domain image.

The gradual attenuation of the filter is important. A simple high-pass filter would simply mask a set of pixels in the frequency domain, but the abrupt transition would cause a ringing effect in the spatial domain.

Arguments:
  • low_cutoff : (list) the cutoff frequencies for x, y and z
  • high_cutoff : (list) the cutoff frequencies for x, y and z
  • order : (int) order determines sharpness of the cutoff curve
def smooth_gaussian(self, sigma=(2, 2, 2), radius=None):
1211    def smooth_gaussian(self, sigma=(2, 2, 2), radius=None):
1212        """
1213        Performs a convolution of the input Volume with a gaussian.
1214
1215        Arguments:
1216            sigma : (float, list)
1217                standard deviation(s) in voxel units.
1218                A list can be given to smooth in the three direction differently.
1219            radius : (float, list)
1220                radius factor(s) determine how far out the gaussian
1221                kernel will go before being clamped to zero. A list can be given too.
1222        """
1223        gsf = vtk.new("ImageGaussianSmooth")
1224        gsf.SetDimensionality(3)
1225        gsf.SetInputData(self.dataset)
1226        if utils.is_sequence(sigma):
1227            gsf.SetStandardDeviations(sigma)
1228        else:
1229            gsf.SetStandardDeviation(sigma)
1230        if radius is not None:
1231            if utils.is_sequence(radius):
1232                gsf.SetRadiusFactors(radius)
1233            else:
1234                gsf.SetRadiusFactor(radius)
1235        gsf.Update()
1236        self._update(gsf.GetOutput())
1237        self.pipeline = utils.OperationNode("smooth_gaussian", parents=[self], c="#4cc9f0")
1238        return self

Performs a convolution of the input Volume with a gaussian.

Arguments:
  • sigma : (float, list) standard deviation(s) in voxel units. A list can be given to smooth in the three direction differently.
  • radius : (float, list) radius factor(s) determine how far out the gaussian kernel will go before being clamped to zero. A list can be given too.
def smooth_median(self, neighbours=(2, 2, 2)):
1240    def smooth_median(self, neighbours=(2, 2, 2)):
1241        """
1242        Median filter that replaces each pixel with the median value
1243        from a rectangular neighborhood around that pixel.
1244        """
1245        imgm = vtk.new("ImageMedian3D")
1246        imgm.SetInputData(self.dataset)
1247        if utils.is_sequence(neighbours):
1248            imgm.SetKernelSize(neighbours[0], neighbours[1], neighbours[2])
1249        else:
1250            imgm.SetKernelSize(neighbours, neighbours, neighbours)
1251        imgm.Update()
1252        self._update(imgm.GetOutput())
1253        self.pipeline = utils.OperationNode("smooth_median", parents=[self], c="#4cc9f0")
1254        return self

Median filter that replaces each pixel with the median value from a rectangular neighborhood around that pixel.

def erode(self, neighbours=(2, 2, 2)):
1256    def erode(self, neighbours=(2, 2, 2)):
1257        """
1258        Replace a voxel with the minimum over an ellipsoidal neighborhood of voxels.
1259        If `neighbours` of an axis is 1, no processing is done on that axis.
1260
1261        Examples:
1262            - [erode_dilate.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/erode_dilate.py)
1263
1264                ![](https://vedo.embl.es/images/volumetric/erode_dilate.png)
1265        """
1266        ver = vtk.new("ImageContinuousErode3D")
1267        ver.SetInputData(self.dataset)
1268        ver.SetKernelSize(neighbours[0], neighbours[1], neighbours[2])
1269        ver.Update()
1270        self._update(ver.GetOutput())
1271        self.pipeline = utils.OperationNode("erode", parents=[self], c="#4cc9f0")
1272        return self

Replace a voxel with the minimum over an ellipsoidal neighborhood of voxels. If neighbours of an axis is 1, no processing is done on that axis.

Examples:
def dilate(self, neighbours=(2, 2, 2)):
1274    def dilate(self, neighbours=(2, 2, 2)):
1275        """
1276        Replace a voxel with the maximum over an ellipsoidal neighborhood of voxels.
1277        If `neighbours` of an axis is 1, no processing is done on that axis.
1278
1279        Check also `erode()` and `pad()`.
1280
1281        Examples:
1282            - [erode_dilate.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/erode_dilate.py)
1283        """
1284        ver = vtk.new("ImageContinuousDilate3D")
1285        ver.SetInputData(self.dataset)
1286        ver.SetKernelSize(neighbours[0], neighbours[1], neighbours[2])
1287        ver.Update()
1288        self._update(ver.GetOutput())
1289        self.pipeline = utils.OperationNode("dilate", parents=[self], c="#4cc9f0")
1290        return self

Replace a voxel with the maximum over an ellipsoidal neighborhood of voxels. If neighbours of an axis is 1, no processing is done on that axis.

Check also erode() and pad().

Examples:
def magnitude(self):
1292    def magnitude(self):
1293        """Colapses components with magnitude function."""
1294        imgm = vtk.new("ImageMagnitude")
1295        imgm.SetInputData(self.dataset)
1296        imgm.Update()
1297        self._update(imgm.GetOutput())
1298        self.pipeline = utils.OperationNode("magnitude", parents=[self], c="#4cc9f0")
1299        return self

Colapses components with magnitude function.

def topoints(self):
1301    def topoints(self):
1302        """
1303        Extract all image voxels as points.
1304        This function takes an input `Volume` and creates an `Mesh`
1305        that contains the points and the point attributes.
1306
1307        Examples:
1308            - [vol2points.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/vol2points.py)
1309        """
1310        v2p = vtk.new("ImageToPoints")
1311        v2p.SetInputData(self.dataset)
1312        v2p.Update()
1313        mpts = vedo.Points(v2p.GetOutput())
1314        mpts.pipeline = utils.OperationNode("topoints", parents=[self], c="#4cc9f0:#e9c46a")
1315        return mpts

Extract all image voxels as points. This function takes an input Volume and creates an Mesh that contains the points and the point attributes.

Examples:
def euclidean_distance(self, anisotropy=False, max_distance=None):
1317    def euclidean_distance(self, anisotropy=False, max_distance=None):
1318        """
1319        Implementation of the Euclidean DT (Distance Transform) using Saito's algorithm.
1320        The distance map produced contains the square of the Euclidean distance values.
1321        The algorithm has a O(n^(D+1)) complexity over n x n x...x n images in D dimensions.
1322
1323        Check out also: https://en.wikipedia.org/wiki/Distance_transform
1324
1325        Arguments:
1326            anisotropy : bool
1327                used to define whether Spacing should be used in the
1328                computation of the distances.
1329            max_distance : bool
1330                any distance bigger than max_distance will not be
1331                computed but set to this specified value instead.
1332
1333        Examples:
1334            - [euclidian_dist.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/euclidian_dist.py)
1335        """
1336        euv = vtk.new("ImageEuclideanDistance")
1337        euv.SetInputData(self.dataset)
1338        euv.SetConsiderAnisotropy(anisotropy)
1339        if max_distance is not None:
1340            euv.InitializeOn()
1341            euv.SetMaximumDistance(max_distance)
1342        euv.SetAlgorithmToSaito()
1343        euv.Update()
1344        vol = Volume(euv.GetOutput())
1345        vol.pipeline = utils.OperationNode("euclidean_distance", parents=[self], c="#4cc9f0")
1346        return vol

Implementation of the Euclidean DT (Distance Transform) using Saito's algorithm. The distance map produced contains the square of the Euclidean distance values. The algorithm has a O(n^(D+1)) complexity over n x n x...x n images in D dimensions.

Check out also: https://en.wikipedia.org/wiki/Distance_transform

Arguments:
  • anisotropy : bool used to define whether Spacing should be used in the computation of the distances.
  • max_distance : bool any distance bigger than max_distance will not be computed but set to this specified value instead.
Examples:
def correlation_with(self, vol2, dim=2):
1348    def correlation_with(self, vol2, dim=2):
1349        """
1350        Find the correlation between two volumetric data sets.
1351        Keyword `dim` determines whether the correlation will be 3D, 2D or 1D.
1352        The default is a 2D Correlation.
1353
1354        The output size will match the size of the first input.
1355        The second input is considered the correlation kernel.
1356        """
1357        imc = vtk.new("ImageCorrelation")
1358        imc.SetInput1Data(self.dataset)
1359        imc.SetInput2Data(vol2.dataset)
1360        imc.SetDimensionality(dim)
1361        imc.Update()
1362        vol = Volume(imc.GetOutput())
1363
1364        vol.pipeline = utils.OperationNode("correlation_with", parents=[self, vol2], c="#4cc9f0")
1365        return vol

Find the correlation between two volumetric data sets. Keyword dim determines whether the correlation will be 3D, 2D or 1D. The default is a 2D Correlation.

The output size will match the size of the first input. The second input is considered the correlation kernel.

def scale_voxels(self, scale=1):
1367    def scale_voxels(self, scale=1):
1368        """Scale the voxel content by factor `scale`."""
1369        rsl = vtk.new("ImageReslice")
1370        rsl.SetInputData(self.dataset)
1371        rsl.SetScalarScale(scale)
1372        rsl.Update()
1373        self._update(rsl.GetOutput())
1374        self.pipeline = utils.OperationNode(
1375            "scale_voxels", comment=f"scale={scale}", parents=[self], c="#4cc9f0"
1376        )
1377        return self

Scale the voxel content by factor scale.