vedo.volume

   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