vedo.volume

   1import glob
   2import os
   3from deprecated import deprecated
   4
   5import numpy as np
   6
   7try:
   8    import vedo.vtkclasses as vtk
   9except ImportError:
  10    import vtkmodules.all as vtk
  11
  12import vedo
  13from vedo import utils
  14from vedo.base import Base3DProp
  15from vedo.base import BaseGrid
  16from vedo.mesh import Mesh
  17
  18
  19__docformat__ = "google"
  20
  21__doc__ = """
  22Work with volumetric datasets (voxel data).
  23
  24![](https://vedo.embl.es/images/volumetric/slicePlane2.png)
  25"""
  26
  27__all__ = [
  28    "BaseVolume",  # included to generate documentation in pydoc
  29    "Volume",
  30    "VolumeSlice",
  31]
  32
  33
  34##########################################################################
  35class BaseVolume:
  36    """
  37    Base class. Do not instantiate.
  38    """
  39    def __init__(self):
  40        """Base class. Do not instantiate."""
  41
  42        self._data = None
  43        self._mapper = None
  44        self.transform = None
  45
  46    def _repr_html_(self):
  47        """
  48        HTML representation of the Volume object for Jupyter Notebooks.
  49        
  50        Returns:
  51            HTML text with the image and some properties.
  52        """
  53        import io
  54        import base64
  55        from PIL import Image
  56
  57        library_name = "vedo.volume.Volume"
  58        help_url = "https://vedo.embl.es/docs/vedo/volume.html"
  59
  60        arr = self.thumbnail(azimuth=0, elevation=-60, zoom=1.4, axes=True)
  61
  62        im = Image.fromarray(arr)
  63        buffered = io.BytesIO()
  64        im.save(buffered, format="PNG", quality=100)
  65        encoded = base64.b64encode(buffered.getvalue()).decode("utf-8")
  66        url = "data:image/png;base64," + encoded
  67        image = f"<img src='{url}'></img>"
  68
  69        # statisitics
  70        bounds = "<br/>".join(
  71            [
  72                utils.precision(min_x,4) + " ... " + utils.precision(max_x,4)
  73                for min_x, max_x in zip(self.bounds()[::2], self.bounds()[1::2])
  74            ]
  75        )
  76
  77        help_text = ""
  78        if self.name:
  79            help_text += f"<b> {self.name}: &nbsp&nbsp</b>"
  80        help_text += '<b><a href="' + help_url + '" target="_blank">' + library_name + "</a></b>" 
  81        if self.filename:
  82            dots = ""
  83            if len(self.filename) > 30:
  84                dots = "..."
  85            help_text += f"<br/><code><i>({dots}{self.filename[-30:]})</i></code>"
  86        
  87        pdata = ""
  88        if self._data.GetPointData().GetScalars():
  89            if self._data.GetPointData().GetScalars().GetName():
  90                name = self._data.GetPointData().GetScalars().GetName()
  91                pdata = "<tr><td><b> point data array </b></td><td>" + name + "</td></tr>"
  92
  93        cdata = ""
  94        if self._data.GetCellData().GetScalars():
  95            if self._data.GetCellData().GetScalars().GetName():
  96                name = self._data.GetCellData().GetScalars().GetName()
  97                cdata = "<tr><td><b> voxel data array </b></td><td>" + name + "</td></tr>"
  98
  99        img = self.GetMapper().GetInput()
 100
 101        all = [
 102            "<table>",
 103            "<tr>", 
 104            "<td>", image, "</td>",
 105            "<td style='text-align: center; vertical-align: center;'><br/>", help_text,
 106            "<table>",
 107            "<tr><td><b> bounds </b> <br/> (x/y/z) </td><td>" + str(bounds) + "</td></tr>",
 108            "<tr><td><b> dimensions </b></td><td>" + str(img.GetDimensions()) + "</td></tr>",
 109            "<tr><td><b> voxel spacing </b></td><td>" + utils.precision(img.GetSpacing(),3) + "</td></tr>",
 110            "<tr><td><b> in memory size </b></td><td>" + str(int(img.GetActualMemorySize() / 1024)) + "MB</td></tr>",
 111            pdata,
 112            cdata,
 113            "<tr><td><b> scalar range </b></td><td>" + utils.precision(img.GetScalarRange(), 4) + "</td></tr>",
 114            "</table>",
 115            "</table>",
 116        ]
 117        return "\n".join(all)
 118
 119    def _update(self, img):
 120        self._data = img
 121        self._data.GetPointData().Modified()
 122        self._mapper.SetInputData(img)
 123        self._mapper.Modified()
 124        self._mapper.Update()
 125        return self
 126
 127    def clone(self):
 128        """Return a clone copy of the Volume."""
 129        newimg = vtk.vtkImageData()
 130        newimg.CopyStructure(self._data)
 131        newimg.CopyAttributes(self._data)
 132
 133        newvol = Volume(newimg)
 134        prop = vtk.vtkVolumeProperty()
 135        prop.DeepCopy(self.GetProperty())
 136        newvol.SetProperty(prop)
 137        newvol.SetOrigin(self.GetOrigin())
 138        newvol.SetScale(self.GetScale())
 139        newvol.SetOrientation(self.GetOrientation())
 140        newvol.SetPosition(self.GetPosition())
 141        newvol.pipeline = utils.OperationNode(
 142            "clone", parents=[self], c="#bbd0ff", shape="diamond")
 143        return newvol
 144
 145    def imagedata(self):
 146        """Return the underlying `vtkImagaData` object."""
 147        return self._data
 148
 149    def tonumpy(self):
 150        """
 151        Get read-write access to voxels of a Volume object as a numpy array.
 152
 153        When you set values in the output image, you don't want numpy to reallocate the array
 154        but instead set values in the existing array, so use the [:] operator.
 155
 156        Example:
 157            `arr[:] = arr*2 + 15`
 158
 159        If the array is modified add a call to:
 160        `volume.imagedata().GetPointData().GetScalars().Modified()`
 161        when all your modifications are completed.
 162        """
 163        narray_shape = tuple(reversed(self._data.GetDimensions()))
 164
 165        scals = self._data.GetPointData().GetScalars()
 166        comps = scals.GetNumberOfComponents()
 167        if comps == 1:
 168            narray = utils.vtk2numpy(scals).reshape(narray_shape)
 169            narray = np.transpose(narray, axes=[2, 1, 0])
 170        else:
 171            narray = utils.vtk2numpy(scals).reshape(*narray_shape, comps)
 172            narray = np.transpose(narray, axes=[2, 1, 0, 3])
 173
 174        # narray = utils.vtk2numpy(self._data.GetPointData().GetScalars()).reshape(narray_shape)
 175        # narray = np.transpose(narray, axes=[2, 1, 0])
 176
 177        return narray
 178
 179    def dimensions(self):
 180        """Return the nr. of voxels in the 3 dimensions."""
 181        return np.array(self._data.GetDimensions())
 182
 183    @deprecated(reason=vedo.colors.red + "Please use scalar_range()" + vedo.colors.reset)
 184    def scalarRange(self):
 185        "Deprecated. Please use `scalar_range()`."
 186        return self.scalar_range()
 187
 188    def scalar_range(self):
 189        """Return the range of the scalar values."""
 190        return np.array(self._data.GetScalarRange())
 191
 192    def spacing(self, s=None):
 193        """Set/get the voxels size in the 3 dimensions."""
 194        if s is not None:
 195            self._data.SetSpacing(s)
 196            return self
 197        return np.array(self._data.GetSpacing())
 198
 199    # def pos(self, p=None): # conflicts with API of base.x()
 200    #     """Same effect as calling `origin()`."""
 201    #     return self.origin(p)
 202
 203    def origin(self, s=None):
 204        """Set/get the origin of the volumetric dataset."""
 205        ### supersedes base.origin()
 206        ### DIFFERENT from base.origin()!
 207        if s is not None:
 208            self._data.SetOrigin(s)
 209            return self
 210        return np.array(self._data.GetOrigin())
 211
 212    def center(self, p=None):
 213        """Set/get the center of the volumetric dataset."""
 214        if p is not None:
 215            self._data.SetCenter(p)
 216            return self
 217        return np.array(self._data.GetCenter())
 218
 219    def permute_axes(self, x, y, z):
 220        """
 221        Reorder the axes of the Volume by specifying
 222        the input axes which are supposed to become the new X, Y, and Z.
 223        """
 224        imp = vtk.vtkImagePermute()
 225        imp.SetFilteredAxes(x, y, z)
 226        imp.SetInputData(self.imagedata())
 227        imp.Update()
 228        self._update(imp.GetOutput())
 229        self.pipeline = utils.OperationNode(
 230            f"permute_axes\n{tuple(x,y,z)}", parents=[self], c="#4cc9f0")
 231        return self
 232
 233    def resample(self, new_spacing, interpolation=1):
 234        """
 235        Resamples a `Volume` to be larger or smaller.
 236
 237        This method modifies the spacing of the input.
 238        Linear interpolation is used to resample the data.
 239
 240        Arguments:
 241            new_spacing : (list)
 242                a list of 3 new spacings for the 3 axes
 243            interpolation : (int)
 244                0=nearest_neighbor, 1=linear, 2=cubic
 245        """
 246        rsp = vtk.vtkImageResample()
 247        oldsp = self.spacing()
 248        for i in range(3):
 249            if oldsp[i] != new_spacing[i]:
 250                rsp.SetAxisOutputSpacing(i, new_spacing[i])
 251        rsp.InterpolateOn()
 252        rsp.SetInterpolationMode(interpolation)
 253        rsp.OptimizationOn()
 254        rsp.Update()
 255        self._update(rsp.GetOutput())
 256        self.pipeline = utils.OperationNode(
 257            f"resample\n{tuple(new_spacing)}", parents=[self], c="#4cc9f0")
 258        return self
 259
 260    def interpolation(self, itype):
 261        """
 262        Set interpolation type.
 263
 264        0=nearest neighbour, 1=linear
 265        """
 266        self.property.SetInterpolationType(itype)
 267        return self
 268
 269    def threshold(self, above=None, below=None, replace=None, replace_value=None):
 270        """
 271        Binary or continuous volume thresholding.
 272        Find the voxels that contain a value above/below the input values
 273        and replace them with a new value (default is 0).
 274        """
 275        th = vtk.vtkImageThreshold()
 276        th.SetInputData(self.imagedata())
 277
 278        # sanity checks
 279        if above is not None and below is not None:
 280            if above == below:
 281                return self
 282            if above > below:
 283                vedo.logger.warning("in volume.threshold(), above > below, skip.")
 284                return self
 285
 286        ## cases
 287        if below is not None and above is not None:
 288            th.ThresholdBetween(above, below)
 289
 290        elif above is not None:
 291            th.ThresholdByUpper(above)
 292
 293        elif below is not None:
 294            th.ThresholdByLower(below)
 295
 296        ##
 297        if replace is not None:
 298            th.SetReplaceIn(True)
 299            th.SetInValue(replace)
 300        else:
 301            th.SetReplaceIn(False)
 302
 303        if replace_value is not None:
 304            th.SetReplaceOut(True)
 305            th.SetOutValue(replace_value)
 306        else:
 307            th.SetReplaceOut(False)
 308
 309        th.Update()
 310        self._update(th.GetOutput())
 311        self.pipeline = utils.OperationNode("threshold", parents=[self], c="#4cc9f0")
 312        return self
 313
 314    def crop(
 315            self,
 316            left=None, right=None,
 317            back=None, front=None,
 318            bottom=None, top=None,
 319            VOI=()
 320        ):
 321        """
 322        Crop a `Volume` object.
 323
 324        Arguments:
 325            left : (float)
 326                fraction to crop from the left plane (negative x)
 327            right : (float)
 328                fraction to crop from the right plane (positive x)
 329            back : (float)
 330                fraction to crop from the back plane (negative y)
 331            front : (float)
 332                fraction to crop from the front plane (positive y)
 333            bottom : (float)
 334                fraction to crop from the bottom plane (negative z)
 335            top : (float)
 336                fraction to crop from the top plane (positive z)
 337            VOI : (list)
 338                extract Volume Of Interest expressed in voxel numbers
 339
 340        Example:
 341            `vol.crop(VOI=(xmin, xmax, ymin, ymax, zmin, zmax)) # all integers nrs`
 342        """
 343        extractVOI = vtk.vtkExtractVOI()
 344        extractVOI.SetInputData(self.imagedata())
 345
 346        if VOI:
 347            extractVOI.SetVOI(VOI)
 348        else:
 349            d = self.imagedata().GetDimensions()
 350            bx0, bx1, by0, by1, bz0, bz1 = 0, d[0]-1, 0, d[1]-1, 0, d[2]-1
 351            if left is not None:   bx0 = int((d[0]-1)*left)
 352            if right is not None:  bx1 = int((d[0]-1)*(1-right))
 353            if back is not None:   by0 = int((d[1]-1)*back)
 354            if front is not None:  by1 = int((d[1]-1)*(1-front))
 355            if bottom is not None: bz0 = int((d[2]-1)*bottom)
 356            if top is not None:    bz1 = int((d[2]-1)*(1-top))
 357            extractVOI.SetVOI(bx0, bx1, by0, by1, bz0, bz1)
 358        extractVOI.Update()
 359        self._update(extractVOI.GetOutput())
 360
 361        self.pipeline = utils.OperationNode(
 362            "crop", parents=[self], c="#4cc9f0",
 363            comment=f"dims={tuple(self.dimensions())}",
 364        )
 365        return self
 366
 367    def append(self, volumes, axis="z", preserve_extents=False):
 368        """
 369        Take the components from multiple inputs and merges them into one output.
 370        Except for the append axis, all inputs must have the same extent.
 371        All inputs must have the same number of scalar components.
 372        The output has the same origin and spacing as the first input.
 373        The origin and spacing of all other inputs are ignored.
 374        All inputs must have the same scalar type.
 375
 376        Arguments:
 377            axis : (int, str)
 378                axis expanded to hold the multiple images
 379            preserve_extents : (bool)
 380                if True, the extent of the inputs is used to place
 381                the image in the output. The whole extent of the output is the union of the input
 382                whole extents. Any portion of the output not covered by the inputs is set to zero.
 383                The origin and spacing is taken from the first input.
 384
 385        Example:
 386            ```python
 387            from vedo import Volume, dataurl
 388            vol = Volume(dataurl+'embryo.tif')
 389            vol.append(vol, axis='x').show().close()
 390            ```
 391            ![](https://vedo.embl.es/images/feats/volume_append.png)
 392        """
 393        ima = vtk.vtkImageAppend()
 394        ima.SetInputData(self.imagedata())
 395        if not utils.is_sequence(volumes):
 396            volumes = [volumes]
 397        for volume in volumes:
 398            if isinstance(volume, vtk.vtkImageData):
 399                ima.AddInputData(volume)
 400            else:
 401                ima.AddInputData(volume.imagedata())
 402        ima.SetPreserveExtents(preserve_extents)
 403        if axis == "x":
 404            axis = 0
 405        elif axis == "y":
 406            axis = 1
 407        elif axis == "z":
 408            axis = 2
 409        ima.SetAppendAxis(axis)
 410        ima.Update()
 411        self._update(ima.GetOutput())
 412        
 413        self.pipeline = utils.OperationNode(
 414            "append", parents=[self, *volumes], c="#4cc9f0",
 415            comment=f"dims={tuple(self.dimensions())}",
 416        )
 417        return self
 418
 419    def resize(self, *newdims):
 420        """Increase or reduce the number of voxels of a Volume with interpolation."""
 421        old_dims = np.array(self.imagedata().GetDimensions())
 422        old_spac = np.array(self.imagedata().GetSpacing())
 423        rsz = vtk.vtkImageResize()
 424        rsz.SetResizeMethodToOutputDimensions()
 425        rsz.SetInputData(self.imagedata())
 426        rsz.SetOutputDimensions(newdims)
 427        rsz.Update()
 428        self._data = rsz.GetOutput()
 429        new_spac = old_spac * old_dims / newdims  # keep aspect ratio
 430        self._data.SetSpacing(new_spac)
 431        self._update(self._data)
 432        self.pipeline = utils.OperationNode(
 433            "resize", parents=[self], c="#4cc9f0",
 434            comment=f"dims={tuple(self.dimensions())}",
 435        )
 436        return self
 437
 438    def normalize(self):
 439        """Normalize that scalar components for each point."""
 440        norm = vtk.vtkImageNormalize()
 441        norm.SetInputData(self.imagedata())
 442        norm.Update()
 443        self._update(norm.GetOutput())
 444        self.pipeline = utils.OperationNode("normalize", parents=[self], c="#4cc9f0")
 445        return self
 446
 447    def mirror(self, axis="x"):
 448        """
 449        Mirror flip along one of the cartesian axes.
 450        """
 451        img = self.imagedata()
 452
 453        ff = vtk.vtkImageFlip()
 454        ff.SetInputData(img)
 455        if axis.lower() == "x":
 456            ff.SetFilteredAxis(0)
 457        elif axis.lower() == "y":
 458            ff.SetFilteredAxis(1)
 459        elif axis.lower() == "z":
 460            ff.SetFilteredAxis(2)
 461        else:
 462            vedo.logger.error("mirror must be set to either x, y, z or n")
 463            raise RuntimeError()
 464        ff.Update()
 465        self._update(ff.GetOutput())
 466        self.pipeline = utils.OperationNode(f"mirror {axis}", parents=[self], c="#4cc9f0")
 467        return self
 468
 469    def operation(self, operation, volume2=None):
 470        """
 471        Perform operations with `Volume` objects.
 472        Keyword `volume2` can be a constant float.
 473
 474        Possible operations are:
 475        ```
 476        +, -, /, 1/x, sin, cos, exp, log,
 477        abs, **2, sqrt, min, max, atan, atan2, median,
 478        mag, dot, gradient, divergence, laplacian.
 479        ```
 480
 481        Examples:
 482            - [volume_operations.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/volume_operations.py)
 483        """
 484        op = operation.lower()
 485        image1 = self._data
 486
 487        mf = None
 488        if op in ["median"]:
 489            mf = vtk.vtkImageMedian3D()
 490            mf.SetInputData(image1)
 491        elif op in ["mag"]:
 492            mf = vtk.vtkImageMagnitude()
 493            mf.SetInputData(image1)
 494        elif op in ["dot", "dotproduct"]:
 495            mf = vtk.vtkImageDotProduct()
 496            mf.SetInput1Data(image1)
 497            mf.SetInput2Data(volume2.imagedata())
 498        elif op in ["grad", "gradient"]:
 499            mf = vtk.vtkImageGradient()
 500            mf.SetDimensionality(3)
 501            mf.SetInputData(image1)
 502        elif op in ["div", "divergence"]:
 503            mf = vtk.vtkImageDivergence()
 504            mf.SetInputData(image1)
 505        elif op in ["laplacian"]:
 506            mf = vtk.vtkImageLaplacian()
 507            mf.SetDimensionality(3)
 508            mf.SetInputData(image1)
 509        
 510        if mf is not None:
 511            mf.Update()
 512            vol = Volume(mf.GetOutput())
 513            vol.pipeline = utils.OperationNode(
 514                f"operation\n{op}", parents=[self], c="#4cc9f0",
 515                shape="cylinder",
 516            )
 517            return vol ###########################
 518
 519        mat = vtk.vtkImageMathematics()
 520        mat.SetInput1Data(image1)
 521
 522        K = None
 523
 524        if utils.is_number(volume2):
 525            K = volume2
 526            mat.SetConstantK(K)
 527            mat.SetConstantC(K)
 528
 529        elif volume2 is not None:  # assume image2 is a constant value
 530            mat.SetInput2Data(volume2.imagedata())
 531
 532        # ###########################
 533        if op in ["+", "add", "plus"]:
 534            if K:
 535                mat.SetOperationToAddConstant()
 536            else:
 537                mat.SetOperationToAdd()
 538
 539        elif op in ["-", "subtract", "minus"]:
 540            if K:
 541                mat.SetConstantC(-float(K))
 542                mat.SetOperationToAddConstant()
 543            else:
 544                mat.SetOperationToSubtract()
 545
 546        elif op in ["*", "multiply", "times"]:
 547            if K:
 548                mat.SetOperationToMultiplyByK()
 549            else:
 550                mat.SetOperationToMultiply()
 551
 552        elif op in ["/", "divide"]:
 553            if K:
 554                mat.SetConstantK(1.0 / K)
 555                mat.SetOperationToMultiplyByK()
 556            else:
 557                mat.SetOperationToDivide()
 558
 559        elif op in ["1/x", "invert"]:
 560            mat.SetOperationToInvert()
 561        elif op in ["sin"]:
 562            mat.SetOperationToSin()
 563        elif op in ["cos"]:
 564            mat.SetOperationToCos()
 565        elif op in ["exp"]:
 566            mat.SetOperationToExp()
 567        elif op in ["log"]:
 568            mat.SetOperationToLog()
 569        elif op in ["abs"]:
 570            mat.SetOperationToAbsoluteValue()
 571        elif op in ["**2", "square"]:
 572            mat.SetOperationToSquare()
 573        elif op in ["sqrt", "sqr"]:
 574            mat.SetOperationToSquareRoot()
 575        elif op in ["min"]:
 576            mat.SetOperationToMin()
 577        elif op in ["max"]:
 578            mat.SetOperationToMax()
 579        elif op in ["atan"]:
 580            mat.SetOperationToATAN()
 581        elif op in ["atan2"]:
 582            mat.SetOperationToATAN2()
 583        else:
 584            vedo.logger.error(f"unknown operation {operation}")
 585            raise RuntimeError()
 586        mat.Update()
 587
 588        self._update(mat.GetOutput())
 589
 590        self.pipeline = utils.OperationNode(
 591            f"operation\n{op}", parents=[self, volume2],
 592            shape="cylinder", c="#4cc9f0",
 593        )
 594        return self
 595
 596    def frequency_pass_filter(self, low_cutoff=None, high_cutoff=None, order=1):
 597        """
 598        Low-pass and high-pass filtering become trivial in the frequency domain.
 599        A portion of the pixels/voxels are simply masked or attenuated.
 600        This function applies a high pass Butterworth filter that attenuates the
 601        frequency domain image.
 602
 603        The gradual attenuation of the filter is important.
 604        A simple high-pass filter would simply mask a set of pixels in the frequency domain,
 605        but the abrupt transition would cause a ringing effect in the spatial domain.
 606
 607        Arguments:
 608            low_cutoff : (list)
 609                the cutoff frequencies for x, y and z
 610            high_cutoff : (list)
 611                the cutoff frequencies for x, y and z
 612            order : (int)
 613                order determines sharpness of the cutoff curve
 614        """
 615        # https://lorensen.github.io/VTKExamples/site/Cxx/ImageProcessing/IdealHighPass
 616        fft = vtk.vtkImageFFT()
 617        fft.SetInputData(self._data)
 618        fft.Update()
 619        out = fft.GetOutput()
 620
 621        if high_cutoff:
 622            blp = vtk.vtkImageButterworthLowPass()
 623            blp.SetInputData(out)
 624            blp.SetCutOff(high_cutoff)
 625            blp.SetOrder(order)
 626            blp.Update()
 627            out = blp.GetOutput()
 628
 629        if low_cutoff:
 630            bhp = vtk.vtkImageButterworthHighPass()
 631            bhp.SetInputData(out)
 632            bhp.SetCutOff(low_cutoff)
 633            bhp.SetOrder(order)
 634            bhp.Update()
 635            out = bhp.GetOutput()
 636
 637        rfft = vtk.vtkImageRFFT()
 638        rfft.SetInputData(out)
 639        rfft.Update()
 640
 641        ecomp = vtk.vtkImageExtractComponents()
 642        ecomp.SetInputData(rfft.GetOutput())
 643        ecomp.SetComponents(0)
 644        ecomp.Update()
 645        self._update(ecomp.GetOutput())
 646        self.pipeline = utils.OperationNode("frequency_pass_filter", parents=[self], c="#4cc9f0")
 647        return self
 648
 649    def smooth_gaussian(self, sigma=(2, 2, 2), radius=None):
 650        """
 651        Performs a convolution of the input Volume with a gaussian.
 652
 653        Arguments:
 654            sigma : (float, list)
 655                standard deviation(s) in voxel units.
 656                A list can be given to smooth in the three direction differently.
 657            radius : (float, list)
 658                radius factor(s) determine how far out the gaussian
 659                kernel will go before being clamped to zero. A list can be given too.
 660        """
 661        gsf = vtk.vtkImageGaussianSmooth()
 662        gsf.SetDimensionality(3)
 663        gsf.SetInputData(self.imagedata())
 664        if utils.is_sequence(sigma):
 665            gsf.SetStandardDeviations(sigma)
 666        else:
 667            gsf.SetStandardDeviation(sigma)
 668        if radius is not None:
 669            if utils.is_sequence(radius):
 670                gsf.SetRadiusFactors(radius)
 671            else:
 672                gsf.SetRadiusFactor(radius)
 673        gsf.Update()
 674        self._update(gsf.GetOutput())
 675        self.pipeline = utils.OperationNode("smooth_gaussian", parents=[self], c="#4cc9f0")
 676        return self
 677
 678    def smooth_median(self, neighbours=(2, 2, 2)):
 679        """
 680        Median filter that replaces each pixel with the median value
 681        from a rectangular neighborhood around that pixel.
 682        """
 683        imgm = vtk.vtkImageMedian3D()
 684        imgm.SetInputData(self.imagedata())
 685        if utils.is_sequence(neighbours):
 686            imgm.SetKernelSize(neighbours[0], neighbours[1], neighbours[2])
 687        else:
 688            imgm.SetKernelSize(neighbours, neighbours, neighbours)
 689        imgm.Update()
 690        self._update(imgm.GetOutput())
 691        self.pipeline = utils.OperationNode("smooth_median", parents=[self], c="#4cc9f0")
 692        return self
 693
 694    def erode(self, neighbours=(2, 2, 2)):
 695        """
 696        Replace a voxel with the minimum over an ellipsoidal neighborhood of voxels.
 697        If `neighbours` of an axis is 1, no processing is done on that axis.
 698
 699        Examples:
 700            - [erode_dilate.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/erode_dilate.py)
 701        
 702                ![](https://vedo.embl.es/images/volumetric/erode_dilate.png)
 703        """
 704        ver = vtk.vtkImageContinuousErode3D()
 705        ver.SetInputData(self._data)
 706        ver.SetKernelSize(neighbours[0], neighbours[1], neighbours[2])
 707        ver.Update()
 708        self._update(ver.GetOutput())
 709        self.pipeline = utils.OperationNode("erode", parents=[self], c="#4cc9f0")
 710        return self
 711
 712    def dilate(self, neighbours=(2, 2, 2)):
 713        """
 714        Replace a voxel with the maximum over an ellipsoidal neighborhood of voxels.
 715        If `neighbours` of an axis is 1, no processing is done on that axis.
 716
 717        Examples:
 718            - [erode_dilate.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/erode_dilate.py)
 719        """
 720        ver = vtk.vtkImageContinuousDilate3D()
 721        ver.SetInputData(self._data)
 722        ver.SetKernelSize(neighbours[0], neighbours[1], neighbours[2])
 723        ver.Update()
 724        self._update(ver.GetOutput())
 725        self.pipeline = utils.OperationNode("dilate", parents=[self], c="#4cc9f0")
 726        return self
 727
 728    def magnitude(self):
 729        """Colapses components with magnitude function."""
 730        imgm = vtk.vtkImageMagnitude()
 731        imgm.SetInputData(self.imagedata())
 732        imgm.Update()
 733        self._update(imgm.GetOutput())
 734        self.pipeline = utils.OperationNode("magnitude", parents=[self], c="#4cc9f0")
 735        return self
 736
 737    def topoints(self):
 738        """
 739        Extract all image voxels as points.
 740        This function takes an input `Volume` and creates an `Mesh`
 741        that contains the points and the point attributes.
 742
 743        Examples:
 744            - [vol2points.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/vol2points.py)
 745        """
 746        v2p = vtk.vtkImageToPoints()
 747        v2p.SetInputData(self.imagedata())
 748        v2p.Update()
 749        mpts = vedo.Points(v2p.GetOutput())
 750        mpts.pipeline = utils.OperationNode("topoints", parents=[self], c="#4cc9f0:#e9c46a")
 751        return mpts
 752
 753    def euclidean_distance(self, anisotropy=False, max_distance=None):
 754        """
 755        Implementation of the Euclidean DT (Distance Transform) using Saito's algorithm.
 756        The distance map produced contains the square of the Euclidean distance values.
 757        The algorithm has a O(n^(D+1)) complexity over n x n x...x n images in D dimensions.
 758
 759        Check out also: https://en.wikipedia.org/wiki/Distance_transform
 760
 761        Arguments:
 762            anisotropy : bool
 763                used to define whether Spacing should be used in the
 764                computation of the distances.
 765            max_distance : bool
 766                any distance bigger than max_distance will not be
 767                computed but set to this specified value instead.
 768
 769        Examples:
 770            - [euclidian_dist.py](examples/volumetric/euclidian_dist.py)
 771        """
 772        euv = vtk.vtkImageEuclideanDistance()
 773        euv.SetInputData(self._data)
 774        euv.SetConsiderAnisotropy(anisotropy)
 775        if max_distance is not None:
 776            euv.InitializeOn()
 777            euv.SetMaximumDistance(max_distance)
 778        euv.SetAlgorithmToSaito()
 779        euv.Update()
 780        vol = Volume(euv.GetOutput())
 781        vol.pipeline = utils.OperationNode("euclidean_distance", parents=[self], c="#4cc9f0")
 782        return vol
 783
 784    def correlation_with(self, vol2, dim=2):
 785        """
 786        Find the correlation between two volumetric data sets.
 787        Keyword `dim` determines whether the correlation will be 3D, 2D or 1D.
 788        The default is a 2D Correlation.
 789
 790        The output size will match the size of the first input.
 791        The second input is considered the correlation kernel.
 792        """
 793        imc = vtk.vtkImageCorrelation()
 794        imc.SetInput1Data(self._data)
 795        imc.SetInput2Data(vol2.imagedata())
 796        imc.SetDimensionality(dim)
 797        imc.Update()
 798        vol = Volume(imc.GetOutput())
 799
 800        vol.pipeline = utils.OperationNode(
 801            "correlation_with", parents=[self, vol2], c="#4cc9f0",
 802        )
 803        return vol
 804
 805    def scale_voxels(self, scale=1):
 806        """Scale the voxel content by factor `scale`."""
 807        rsl = vtk.vtkImageReslice()
 808        rsl.SetInputData(self.imagedata())
 809        rsl.SetScalarScale(scale)
 810        rsl.Update()
 811        self._update(rsl.GetOutput())
 812        self.pipeline = utils.OperationNode(
 813            f"scale_voxels\nscale={scale}", parents=[self], c="#4cc9f0")
 814        return self
 815
 816
 817##########################################################################
 818class Volume(BaseVolume, BaseGrid, vtk.vtkVolume):
 819    """
 820    Class to describe dataset that are defined on "voxels":
 821    the 3D equivalent of 2D pixels.
 822    """
 823    def __init__(
 824        self,
 825        inputobj=None,
 826        c="RdBu_r",
 827        alpha=(0.0, 0.0, 0.2, 0.4, 0.8, 1.0),
 828        alpha_gradient=None,
 829        alpha_unit=1,
 830        mode=0,
 831        spacing=None,
 832        dims=None,
 833        origin=None,
 834        mapper="smart",
 835    ):
 836        """
 837        This class can be initialized with a numpy object, a `vtkImageData`
 838        or a list of 2D bmp files.
 839
 840        Arguments:
 841            c : (list, str)
 842                sets colors along the scalar range, or a matplotlib color map name
 843            alphas : (float, list)
 844                sets transparencies along the scalar range
 845            alpha_unit : (float)
 846                low values make composite rendering look brighter and denser
 847            origin : (list)
 848                set volume origin coordinates
 849            spacing : (list)
 850                voxel dimensions in x, y and z.
 851            dims : (list)
 852                specify the dimensions of the volume.
 853            mapper : (str)
 854                either 'gpu', 'opengl_gpu', 'fixed' or 'smart'
 855            mode : (int)
 856                define the volumetric rendering style:
 857                    - 0, composite rendering
 858                    - 1, maximum projection
 859                    - 2, minimum projection
 860                    - 3, average projection
 861                    - 4, additive mode
 862
 863                <br>The default mode is "composite" where the scalar values are sampled through
 864                the volume and composited in a front-to-back scheme through alpha blending.
 865                The final color and opacity is determined using the color and opacity transfer
 866                functions specified in alpha keyword.
 867
 868                Maximum and minimum intensity blend modes use the maximum and minimum
 869                scalar values, respectively, along the sampling ray.
 870                The final color and opacity is determined by passing the resultant value
 871                through the color and opacity transfer functions.
 872
 873                Additive blend mode accumulates scalar values by passing each value
 874                through the opacity transfer function and then adding up the product
 875                of the value and its opacity. In other words, the scalar values are scaled
 876                using the opacity transfer function and summed to derive the final color.
 877                Note that the resulting image is always grayscale i.e. aggregated values
 878                are not passed through the color transfer function.
 879                This is because the final value is a derived value and not a real data value
 880                along the sampling ray.
 881
 882                Average intensity blend mode works similar to the additive blend mode where
 883                the scalar values are multiplied by opacity calculated from the opacity
 884                transfer function and then added.
 885                The additional step here is to divide the sum by the number of samples
 886                taken through the volume.
 887                As is the case with the additive intensity projection, the final image will
 888                always be grayscale i.e. the aggregated values are not passed through the
 889                color transfer function.
 890
 891        Example:
 892            ```python
 893            from vedo import Volume
 894            vol = Volume("path/to/mydata/rec*.bmp")
 895            vol.show()
 896            ```
 897
 898        Examples:
 899            - [numpy2volume1.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/numpy2volume1.py)
 900
 901                ![](https://vedo.embl.es/images/volumetric/numpy2volume1.png)
 902
 903            - [read_volume2.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/read_volume2.py)
 904
 905                ![](https://vedo.embl.es/images/volumetric/read_volume2.png)
 906
 907        .. note:: 
 908            if a `list` of values is used for `alphas` this is interpreted
 909            as a transfer function along the range of the scalar.
 910        """
 911        vtk.vtkVolume.__init__(self)
 912        BaseGrid.__init__(self)
 913        BaseVolume.__init__(self)
 914        # super().__init__()
 915        
 916        ###################
 917        if isinstance(inputobj, str):
 918
 919            if "https://" in inputobj:
 920                inputobj = vedo.io.download(inputobj, verbose=False)  # fpath
 921            elif os.path.isfile(inputobj):
 922                pass
 923            else:
 924                inputobj = sorted(glob.glob(inputobj))
 925
 926        ###################
 927        if "gpu" in mapper:
 928            self._mapper = vtk.vtkGPUVolumeRayCastMapper()
 929        elif "opengl_gpu" in mapper:
 930            self._mapper = vtk.vtkOpenGLGPUVolumeRayCastMapper()
 931        elif "smart" in mapper:
 932            self._mapper = vtk.vtkSmartVolumeMapper()
 933        elif "fixed" in mapper:
 934            self._mapper = vtk.vtkFixedPointVolumeRayCastMapper()
 935        elif isinstance(mapper, vtk.vtkMapper):
 936            self._mapper = mapper
 937        else:
 938            print("Error unknown mapper type", [mapper])
 939            raise RuntimeError()
 940        self.SetMapper(self._mapper)
 941
 942        ###################
 943        inputtype = str(type(inputobj))
 944
 945        # print('Volume inputtype', inputtype, c='b')
 946
 947        if inputobj is None:
 948            img = vtk.vtkImageData()
 949
 950        elif utils.is_sequence(inputobj):
 951
 952            if isinstance(inputobj[0], str) and ".bmp" in inputobj[0].lower():
 953                # scan sequence of BMP files
 954                ima = vtk.vtkImageAppend()
 955                ima.SetAppendAxis(2)
 956                pb = utils.ProgressBar(0, len(inputobj))
 957                for i in pb.range():
 958                    f = inputobj[i]
 959                    if "_rec_spr.bmp" in f:
 960                        continue
 961                    picr = vtk.vtkBMPReader()
 962                    picr.SetFileName(f)
 963                    picr.Update()
 964                    mgf = vtk.vtkImageMagnitude()
 965                    mgf.SetInputData(picr.GetOutput())
 966                    mgf.Update()
 967                    ima.AddInputData(mgf.GetOutput())
 968                    pb.print("loading...")
 969                ima.Update()
 970                img = ima.GetOutput()
 971
 972            else:
 973
 974                if len(inputobj.shape) == 1:
 975                    varr = utils.numpy2vtk(inputobj)
 976                else:
 977                    varr = utils.numpy2vtk(inputobj.ravel(order="F"))
 978                varr.SetName("input_scalars")
 979
 980                img = vtk.vtkImageData()
 981                if dims is not None:
 982                    img.SetDimensions(dims[2],dims[1],dims[0])
 983                else:
 984                    if len(inputobj.shape)==1:
 985                        vedo.logger.error("must set dimensions (dims keyword) in Volume")
 986                        raise RuntimeError()
 987                    img.SetDimensions(inputobj.shape)
 988                img.GetPointData().AddArray(varr)
 989                img.GetPointData().SetActiveScalars(varr.GetName())
 990
 991        elif "ImageData" in inputtype:
 992            img = inputobj
 993
 994        elif isinstance(inputobj, Volume):
 995            img = inputobj.inputdata()
 996
 997        elif "UniformGrid" in inputtype:
 998            img = inputobj
 999
1000        elif hasattr(inputobj, "GetOutput"):  # passing vtk object, try extract imagdedata
1001            if hasattr(inputobj, "Update"):
1002                inputobj.Update()
1003            img = inputobj.GetOutput()
1004
1005        elif isinstance(inputobj, str):
1006            if "https://" in inputobj:
1007                inputobj = vedo.io.download(inputobj, verbose=False)
1008            img = vedo.io.loadImageData(inputobj)
1009
1010        else:
1011            vedo.logger.error(f"cannot understand input type {inputtype}")
1012            return
1013
1014        if dims is not None:
1015            img.SetDimensions(dims)
1016
1017        if origin is not None:
1018            img.SetOrigin(origin)  ### DIFFERENT from volume.origin()!
1019
1020        if spacing is not None:
1021            img.SetSpacing(spacing)
1022
1023        self._data = img
1024        self._mapper.SetInputData(img)
1025
1026        if img.GetPointData().GetScalars():
1027            if img.GetPointData().GetScalars().GetNumberOfComponents() == 1:
1028                self.mode(mode).color(c).alpha(alpha).alpha_gradient(alpha_gradient)
1029                self.GetProperty().SetShade(True)
1030                self.GetProperty().SetInterpolationType(1)
1031                self.GetProperty().SetScalarOpacityUnitDistance(alpha_unit)
1032
1033        # remember stuff:
1034        self._mode = mode
1035        self._color = c
1036        self._alpha = alpha
1037        self._alpha_grad = alpha_gradient
1038        self._alpha_unit = alpha_unit
1039        
1040        self.pipeline = utils.OperationNode("Volume",
1041            comment=f"dims={tuple(self.dimensions())}", c="#4cc9f0",
1042        )
1043        #######################################################################
1044
1045
1046    def _update(self, data):
1047        self._data = data
1048        self._data.GetPointData().Modified()
1049        self._mapper.SetInputData(data)
1050        self._mapper.Modified()
1051        self._mapper.Update()
1052        return self
1053
1054    def mode(self, mode=None):
1055        """
1056        Define the volumetric rendering mode following this:
1057            - 0, composite rendering
1058            - 1, maximum projection rendering
1059            - 2, minimum projection rendering
1060            - 3, average projection rendering
1061            - 4, additive mode
1062
1063        The default mode is "composite" where the scalar values are sampled through
1064        the volume and composited in a front-to-back scheme through alpha blending.
1065        The final color and opacity is determined using the color and opacity transfer
1066        functions specified in alpha keyword.
1067
1068        Maximum and minimum intensity blend modes use the maximum and minimum
1069        scalar values, respectively, along the sampling ray.
1070        The final color and opacity is determined by passing the resultant value
1071        through the color and opacity transfer functions.
1072
1073        Additive blend mode accumulates scalar values by passing each value
1074        through the opacity transfer function and then adding up the product
1075        of the value and its opacity. In other words, the scalar values are scaled
1076        using the opacity transfer function and summed to derive the final color.
1077        Note that the resulting image is always grayscale i.e. aggregated values
1078        are not passed through the color transfer function.
1079        This is because the final value is a derived value and not a real data value
1080        along the sampling ray.
1081
1082        Average intensity blend mode works similar to the additive blend mode where
1083        the scalar values are multiplied by opacity calculated from the opacity
1084        transfer function and then added.
1085        The additional step here is to divide the sum by the number of samples
1086        taken through the volume.
1087        As is the case with the additive intensity projection, the final image will
1088        always be grayscale i.e. the aggregated values are not passed through the
1089        color transfer function.
1090        """
1091        if mode is None:
1092            return self._mapper.GetBlendMode()
1093
1094        if isinstance(mode, str):
1095            if "comp" in mode:
1096                mode = 0
1097            elif "proj" in mode:
1098                if "max" in mode:
1099                    mode = 1
1100                elif "min" in mode:
1101                    mode = 2
1102                elif "ave" in mode:
1103                    mode = 3
1104                else:
1105                    vedo.logger.warning(f"unknown mode {mode}")
1106                    mode = 0
1107            elif "add" in mode:
1108                mode = 4
1109            else:
1110                vedo.logger.warning(f"unknown mode {mode}")
1111                mode = 0
1112
1113        self._mapper.SetBlendMode(mode)
1114        self._mode = mode
1115        return self
1116
1117    def shade(self, status=None):
1118        """
1119        Set/Get the shading of a Volume.
1120        Shading can be further controlled with `volume.lighting()` method.
1121
1122        If shading is turned on, the mapper may perform shading calculations.
1123        In some cases shading does not apply
1124        (for example, in maximum intensity projection mode).
1125        """
1126        if status is None:
1127            return self.GetProperty().GetShade()
1128        self.GetProperty().SetShade(status)
1129        return self
1130
1131    def cmap(self, c, alpha=None, vmin=None, vmax=None):
1132        """Same as `color()`.
1133
1134        Arguments:
1135            alpha : (list)
1136                use a list to specify transparencies along the scalar range
1137            vmin : (float)
1138                force the min of the scalar range to be this value
1139            vmax : (float)
1140                force the max of the scalar range to be this value
1141        """
1142        return self.color(c, alpha, vmin, vmax)
1143
1144    def jittering(self, status=None):
1145        """
1146        If `True`, each ray traversal direction will be perturbed slightly
1147        using a noise-texture to get rid of wood-grain effects.
1148        """
1149        if hasattr(self._mapper, "SetUseJittering"):  # tetmesh doesnt have it
1150            if status is None:
1151                return self._mapper.GetUseJittering()
1152            self._mapper.SetUseJittering(status)
1153        return self
1154
1155    def mask(self, data):
1156        """
1157        Mask a volume visualization with a binary value.
1158        Needs to specify keyword mapper='gpu'.
1159
1160        Example:
1161        ```python
1162            from vedo import np, Volume, show
1163            data_matrix = np.zeros([75, 75, 75], dtype=np.uint8)
1164            # all voxels have value zero except:
1165            data_matrix[0:35,   0:35,  0:35] = 1
1166            data_matrix[35:55, 35:55, 35:55] = 2
1167            data_matrix[55:74, 55:74, 55:74] = 3
1168            vol = Volume(data_matrix, c=['white','b','g','r'], mapper='gpu')
1169            data_mask = np.zeros_like(data_matrix)
1170            data_mask[10:65, 10:45, 20:75] = 1
1171            vol.mask(data_mask)
1172            show(vol, axes=1).close()
1173        ```
1174        """
1175        mask = Volume(data.astype(np.uint8))
1176        try:
1177            self.mapper().SetMaskTypeToBinary()
1178            self.mapper().SetMaskInput(mask.inputdata())
1179        except AttributeError:
1180            vedo.logger.error("volume.mask() must create the volume with Volume(..., mapper='gpu')")
1181        return self
1182
1183    def alpha_gradient(self, alpha_grad, vmin=None, vmax=None):
1184        """
1185        Assign a set of tranparencies to a volume's gradient
1186        along the range of the scalar value.
1187        A single constant value can also be assigned.
1188        The gradient function is used to decrease the opacity
1189        in the "flat" regions of the volume while maintaining the opacity
1190        at the boundaries between material types.  The gradient is measured
1191        as the amount by which the intensity changes over unit distance.
1192
1193        The format for alpha_grad is the same as for method `volume.alpha()`.
1194        """
1195        if vmin is None:
1196            vmin, _ = self._data.GetScalarRange()
1197        if vmax is None:
1198            _, vmax = self._data.GetScalarRange()
1199        self._alpha_grad = alpha_grad
1200        volumeProperty = self.GetProperty()
1201
1202        if alpha_grad is None:
1203            volumeProperty.DisableGradientOpacityOn()
1204            return self
1205
1206        volumeProperty.DisableGradientOpacityOff()
1207
1208        gotf = volumeProperty.GetGradientOpacity()
1209        if utils.is_sequence(alpha_grad):
1210            alpha_grad = np.array(alpha_grad)
1211            if len(alpha_grad.shape)==1: # user passing a flat list e.g. (0.0, 0.3, 0.9, 1)
1212                for i, al in enumerate(alpha_grad):
1213                    xalpha = vmin + (vmax - vmin) * i / (len(alpha_grad) - 1)
1214                    # Create transfer mapping scalar value to gradient opacity
1215                    gotf.AddPoint(xalpha, al)
1216            elif len(alpha_grad.shape) == 2:  # user passing [(x0,alpha0), ...]
1217                gotf.AddPoint(vmin, alpha_grad[0][1])
1218                for xalpha, al in alpha_grad:
1219                    # Create transfer mapping scalar value to opacity
1220                    gotf.AddPoint(xalpha, al)
1221                gotf.AddPoint(vmax, alpha_grad[-1][1])
1222            # print("alpha_grad at", round(xalpha, 1), "\tset to", al)
1223        else:
1224            gotf.AddPoint(vmin, alpha_grad)  # constant alpha_grad
1225            gotf.AddPoint(vmax, alpha_grad)
1226        return self
1227
1228    def component_weight(self, i, weight):
1229        """Set the scalar component weight in range [0,1]."""
1230        self.GetProperty().SetComponentWeight(i, weight)
1231        return self
1232
1233    def xslice(self, i):
1234        """Extract the slice at index `i` of volume along x-axis."""
1235        vslice = vtk.vtkImageDataGeometryFilter()
1236        vslice.SetInputData(self.imagedata())
1237        nx, ny, nz = self.imagedata().GetDimensions()
1238        if i > nx - 1:
1239            i = nx - 1
1240        vslice.SetExtent(i, i, 0, ny, 0, nz)
1241        vslice.Update()
1242        m = Mesh(vslice.GetOutput())
1243        m.pipeline = utils.OperationNode(f"xslice {i}", parents=[self], c="#4cc9f0:#e9c46a")
1244        return m
1245
1246    def yslice(self, j):
1247        """Extract the slice at index `j` of volume along y-axis."""
1248        vslice = vtk.vtkImageDataGeometryFilter()
1249        vslice.SetInputData(self.imagedata())
1250        nx, ny, nz = self.imagedata().GetDimensions()
1251        if j > ny - 1:
1252            j = ny - 1
1253        vslice.SetExtent(0, nx, j, j, 0, nz)
1254        vslice.Update()
1255        m = Mesh(vslice.GetOutput())
1256        m.pipeline = utils.OperationNode(f"yslice {j}", parents=[self], c="#4cc9f0:#e9c46a")
1257        return m
1258
1259    def zslice(self, k):
1260        """Extract the slice at index `i` of volume along z-axis."""
1261        vslice = vtk.vtkImageDataGeometryFilter()
1262        vslice.SetInputData(self.imagedata())
1263        nx, ny, nz = self.imagedata().GetDimensions()
1264        if k > nz - 1:
1265            k = nz - 1
1266        vslice.SetExtent(0, nx, 0, ny, k, k)
1267        vslice.Update()
1268        m = Mesh(vslice.GetOutput())
1269        m.pipeline = utils.OperationNode(f"zslice {k}", parents=[self], c="#4cc9f0:#e9c46a")
1270        return m
1271
1272    @deprecated(reason=vedo.colors.red + "Please use slice_plane()" + vedo.colors.reset)
1273    def slicePlane(self, *a, **b):
1274        "Deprecated. Please use `scalar_range()`"
1275        return self.slice_plane(*a, **b)
1276
1277    def slice_plane(self, origin=(0, 0, 0), normal=(1, 1, 1), autocrop=False):
1278        """
1279        Extract the slice along a given plane position and normal.
1280
1281        Example:
1282            - [slicePlane1.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/slicePlane1.py)
1283
1284                ![](https://vedo.embl.es/images/volumetric/slicePlane1.gif)
1285        """
1286        reslice = vtk.vtkImageReslice()
1287        reslice.SetInputData(self._data)
1288        reslice.SetOutputDimensionality(2)
1289        newaxis = utils.versor(normal)
1290        pos = np.array(origin)
1291        initaxis = (0, 0, 1)
1292        crossvec = np.cross(initaxis, newaxis)
1293        angle = np.arccos(np.dot(initaxis, newaxis))
1294        T = vtk.vtkTransform()
1295        T.PostMultiply()
1296        T.RotateWXYZ(np.rad2deg(angle), crossvec)
1297        T.Translate(pos)
1298        M = T.GetMatrix()
1299        reslice.SetResliceAxes(M)
1300        reslice.SetInterpolationModeToLinear()
1301        reslice.SetAutoCropOutput(not autocrop)
1302        reslice.Update()
1303        vslice = vtk.vtkImageDataGeometryFilter()
1304        vslice.SetInputData(reslice.GetOutput())
1305        vslice.Update()
1306        msh = Mesh(vslice.GetOutput())
1307        msh.SetOrientation(T.GetOrientation())
1308        msh.SetPosition(pos)
1309        msh.pipeline = utils.OperationNode(
1310            "slice_plane", parents=[self], c="#4cc9f0:#e9c46a")
1311        return msh
1312
1313    def warp(self, source, target, sigma=1, mode="3d", fit=False):
1314        """
1315        Warp volume scalars within a Volume by specifying
1316        source and target sets of points.
1317
1318        Arguments:
1319            source : (Points, list)
1320                the list of source points
1321            target : (Points, list)
1322                the list of target points
1323            fit : (bool)
1324                fit/adapt the old bounding box to the warped geometry
1325        """
1326        if isinstance(source, vedo.Points):
1327            source = source.points()
1328        if isinstance(target, vedo.Points):
1329            target = target.points()
1330
1331        ns = len(source)
1332        ptsou = vtk.vtkPoints()
1333        ptsou.SetNumberOfPoints(ns)
1334        for i in range(ns):
1335            ptsou.SetPoint(i, source[i])
1336
1337        nt = len(target)
1338        if ns != nt:
1339            vedo.logger.error(f"#source {ns} != {nt} #target points")
1340            raise RuntimeError()
1341
1342        pttar = vtk.vtkPoints()
1343        pttar.SetNumberOfPoints(nt)
1344        for i in range(ns):
1345            pttar.SetPoint(i, target[i])
1346
1347        T = vtk.vtkThinPlateSplineTransform()
1348        if mode.lower() == "3d":
1349            T.SetBasisToR()
1350        elif mode.lower() == "2d":
1351            T.SetBasisToR2LogR()
1352        else:
1353            vedo.logger.error(f"unknown mode {mode}")
1354            raise RuntimeError()
1355
1356        T.SetSigma(sigma)
1357        T.SetSourceLandmarks(ptsou)
1358        T.SetTargetLandmarks(pttar)
1359        T.Inverse()
1360        self.transform = T
1361        self.apply_transform(T, fit=fit)
1362        self.pipeline = utils.OperationNode("warp", parents=[self], c="#4cc9f0")
1363        return self
1364
1365    def apply_transform(self, T, fit=False):
1366        """
1367        Apply a transform to the scalars in the volume.
1368
1369        Arguments:
1370            T : (vtkTransform, matrix)
1371                The transformation to be applied
1372            fit : (bool)
1373                fit/adapt the old bounding box to the warped geometry
1374        """
1375        if isinstance(T, vtk.vtkMatrix4x4):
1376            tr = vtk.vtkTransform()
1377            tr.SetMatrix(T)
1378            T = tr
1379
1380        elif utils.is_sequence(T):
1381            M = vtk.vtkMatrix4x4()
1382            n = len(T[0])
1383            for i in range(n):
1384                for j in range(n):
1385                    M.SetElement(i, j, T[i][j])
1386            tr = vtk.vtkTransform()
1387            tr.SetMatrix(M)
1388            T = tr
1389
1390        reslice = vtk.vtkImageReslice()
1391        reslice.SetInputData(self._data)
1392        reslice.SetResliceTransform(T)
1393        reslice.SetOutputDimensionality(3)
1394        reslice.SetInterpolationModeToLinear()
1395
1396        spacing = self._data.GetSpacing()
1397        origin = self._data.GetOrigin()
1398
1399        if fit:
1400            bb = self.box()
1401            if isinstance(T, vtk.vtkThinPlateSplineTransform):
1402                TI = vtk.vtkThinPlateSplineTransform()
1403                TI.DeepCopy(T)
1404                TI.Inverse()
1405            else:
1406                TI = vtk.vtkTransform()
1407                TI.DeepCopy(T)
1408            bb.apply_transform(TI)
1409            bounds = bb.GetBounds()
1410            bounds = (
1411                bounds[0]/spacing[0], bounds[1]/spacing[0],
1412                bounds[2]/spacing[1], bounds[3]/spacing[1],
1413                bounds[4]/spacing[2], bounds[5]/spacing[2],
1414            )
1415            bounds = np.round(bounds).astype(int)
1416            reslice.SetOutputExtent(bounds)
1417            reslice.SetOutputSpacing(spacing[0], spacing[1], spacing[2])
1418            reslice.SetOutputOrigin(origin[0], origin[1], origin[2])
1419
1420        reslice.Update()
1421        self._update(reslice.GetOutput())
1422        self.pipeline = utils.OperationNode("apply_transform", parents=[self], c="#4cc9f0")
1423        return self
1424
1425
1426##########################################################################
1427class VolumeSlice(BaseVolume, Base3DProp, vtk.vtkImageSlice):
1428    """
1429    Derived class of `vtkImageSlice`.
1430    """
1431    def __init__(self, inputobj=None):
1432        """
1433        This class is equivalent to `Volume` except for its representation.
1434        The main purpose of this class is to be used in conjunction with `Volume`
1435        for visualization using `mode="image"`.
1436        """
1437        vtk.vtkImageSlice.__init__(self)
1438        Base3DProp.__init__(self)
1439        BaseVolume.__init__(self)
1440        # super().__init__()
1441
1442        self._mapper = vtk.vtkImageResliceMapper()
1443        self._mapper.SliceFacesCameraOn()
1444        self._mapper.SliceAtFocalPointOn()
1445        self._mapper.SetAutoAdjustImageQuality(False)
1446        self._mapper.BorderOff()
1447
1448        self.lut = None
1449
1450        self.property = vtk.vtkImageProperty()
1451        self.property.SetInterpolationTypeToLinear()
1452        self.SetProperty(self.property)
1453
1454        ###################
1455        if isinstance(inputobj, str):
1456            if "https://" in inputobj:
1457                inputobj = vedo.io.download(inputobj, verbose=False)  # fpath
1458            elif os.path.isfile(inputobj):
1459                pass
1460            else:
1461                inputobj = sorted(glob.glob(inputobj))
1462
1463        ###################
1464        inputtype = str(type(inputobj))
1465
1466        if inputobj is None:
1467            img = vtk.vtkImageData()
1468
1469        if isinstance(inputobj, Volume):
1470            img = inputobj.imagedata()
1471            self.lut = utils.ctf2lut(inputobj)
1472
1473        elif utils.is_sequence(inputobj):
1474
1475            if isinstance(inputobj[0], str):  # scan sequence of BMP files
1476                ima = vtk.vtkImageAppend()
1477                ima.SetAppendAxis(2)
1478                pb = utils.ProgressBar(0, len(inputobj))
1479                for i in pb.range():
1480                    f = inputobj[i]
1481                    picr = vtk.vtkBMPReader()
1482                    picr.SetFileName(f)
1483                    picr.Update()
1484                    mgf = vtk.vtkImageMagnitude()
1485                    mgf.SetInputData(picr.GetOutput())
1486                    mgf.Update()
1487                    ima.AddInputData(mgf.GetOutput())
1488                    pb.print("loading...")
1489                ima.Update()
1490                img = ima.GetOutput()
1491
1492            else:
1493                if "ndarray" not in inputtype:
1494                    inputobj = np.array(inputobj)
1495
1496                if len(inputobj.shape) == 1:
1497                    varr = utils.numpy2vtk(inputobj, dtype=float)
1498                else:
1499                    if len(inputobj.shape) > 2:
1500                        inputobj = np.transpose(inputobj, axes=[2, 1, 0])
1501                    varr = utils.numpy2vtk(inputobj.ravel(order="F"), dtype=float)
1502                varr.SetName("input_scalars")
1503
1504                img = vtk.vtkImageData()
1505                img.SetDimensions(inputobj.shape)
1506                img.GetPointData().AddArray(varr)
1507                img.GetPointData().SetActiveScalars(varr.GetName())
1508
1509        elif "ImageData" in inputtype:
1510            img = inputobj
1511
1512        elif isinstance(inputobj, Volume):
1513            img = inputobj.inputdata()
1514
1515        elif "UniformGrid" in inputtype:
1516            img = inputobj
1517
1518        elif hasattr(inputobj, "GetOutput"): # passing vtk object, try extract imagdedata
1519            if hasattr(inputobj, "Update"):
1520                inputobj.Update()
1521            img = inputobj.GetOutput()
1522
1523        elif isinstance(inputobj, str):
1524            if "https://" in inputobj:
1525                inputobj = vedo.io.download(inputobj, verbose=False)
1526            img = vedo.io.loadImageData(inputobj)
1527
1528        else:
1529            vedo.logger.error(f"cannot understand input type {inputtype}")
1530            return
1531
1532        self._data = img
1533        self._mapper.SetInputData(img)
1534        self.SetMapper(self._mapper)
1535
1536    def bounds(self):
1537        """Return the bounding box as [x0,x1, y0,y1, z0,z1]"""
1538        bns = [0, 0, 0, 0, 0, 0]
1539        self.GetBounds(bns)
1540        return bns
1541
1542    def colorize(self, lut=None, fix_scalar_range=False):
1543        """
1544        Assign a LUT (Look Up Table) to colorize the slice, leave it `None`
1545        to reuse an existing Volume color map.
1546        Use "bw" for automatic black and white.
1547        """
1548        if lut is None and self.lut:
1549            self.property.SetLookupTable(self.lut)
1550        elif isinstance(lut, vtk.vtkLookupTable):
1551            self.property.SetLookupTable(lut)
1552        elif lut == "bw":
1553            self.property.SetLookupTable(None)
1554        self.property.SetUseLookupTableScalarRange(fix_scalar_range)
1555        return self
1556
1557    def alpha(self, value):
1558        """Set opacity to the slice"""
1559        self.property.SetOpacity(value)
1560        return self
1561
1562    def auto_adjust_quality(self, value=True):
1563        """Automatically reduce the rendering quality for greater speed when interacting"""
1564        self._mapper.SetAutoAdjustImageQuality(value)
1565        return self
1566
1567    def slab(self, thickness=0, mode=0, sample_factor=2):
1568        """
1569        Make a thick slice (slab).
1570
1571        Arguments:
1572            thickness : (float)
1573                set the slab thickness, for thick slicing
1574            mode : (int)
1575                The slab type:
1576                    0 = min
1577                    1 = max
1578                    2 = mean
1579                    3 = sum
1580            sample_factor : (float)
1581                Set the number of slab samples to use as a factor of the number of input slices
1582                within the slab thickness. The default value is 2, but 1 will increase speed
1583                with very little loss of quality.
1584        """
1585        self._mapper.SetSlabThickness(thickness)
1586        self._mapper.SetSlabType(mode)
1587        self._mapper.SetSlabSampleFactor(sample_factor)
1588        return self
1589
1590    def face_camera(self, value=True):
1591        """Make the slice always face the camera or not."""
1592        self._mapper.SetSliceFacesCameraOn(value)
1593        return self
1594
1595    def jump_to_nearest_slice(self, value=True):
1596        """
1597        This causes the slicing to occur at the closest slice to the focal point,
1598        instead of the default behavior where a new slice is interpolated between
1599        the original slices.
1600        Nothing happens if the plane is oblique to the original slices."""
1601        self.SetJumpToNearestSlice(value)
1602        return self
1603
1604    def fill_background(self, value=True):
1605        """
1606        Instead of rendering only to the image border, 
1607        render out to the viewport boundary with the background color. 
1608        The background color will be the lowest color on the lookup
1609        table that is being used for the image."""
1610        self._mapper.SetBackground(value)
1611        return self
1612
1613    def lighting(self, window, level, ambient=1.0, diffuse=0.0):
1614        """Assign the values for window and color level."""
1615        self.property.SetColorWindow(window)
1616        self.property.SetColorLevel(level)
1617        self.property.SetAmbient(ambient)
1618        self.property.SetDiffuse(diffuse)
1619        return self