Python vtk.vtkImageData() Examples

def getDepthMapData(self, viewId):

        mapId = self.reader.GetCurrentMapId(viewId)
        if mapId < 0:
            return None, None

        depthImage = vtk.vtkImageData()
        transform = vtk.vtkTransform()
        self.reader.GetDataForMapId(viewId, mapId, depthImage, transform)

        dims = depthImage.GetDimensions()
        d = vnp.getNumpyFromVtk(depthImage, 'ImageScalars')
        d = d.reshape(dims[1], dims[0])
        t = np.array([[transform.GetMatrix().GetElement(r, c) for c in range(4)] for r in range(4)])

        return d, t 

def add_image(self, image):
        img_mapper = vtk.vtkImageMapper()
        img_actor = vtk.vtkActor2D()
        img_data = vtk.vtkImageData()
        img_data.SetDimensions(image.shape[0], image.shape[1], 1)
        img_data.AllocateScalars(vtk.VTK_UNSIGNED_CHAR, 3)
        for x in range(0, image.shape[0]):
            for y in range(0, image.shape[1]):
                pixel = img_data.GetScalarPointer(x, y, 0)
                pixel = np.array(image[x, y, :])
        img_mapper.SetInputData(img_data)
        img_mapper.SetColorWindow(255)
        img_mapper.SetColorLevel(127.5)
        img_actor.SetMapper(img_mapper)
        self.renderer.AddActor(img_actor) 

def _build_image_data(self, img):
        """Internal helper to consturct the output"""
        if self.__needToUpdateOutput:
            # Clean out the output data object
            img.DeepCopy(vtk.vtkImageData())
            self.__needToUpdateOutput = False
        ext = self.__extent
        dims = self.__dims
        nx, ny, nz = ext[dims[0]], ext[dims[1]], ext[dims[2]]
        if not self.__usePointData:
            nx += 1
            ny += 1
            nz += 1
        sx, sy, sz = self.__spacing[0],self.__spacing[1],self.__spacing[2]
        ox, oy, oz = self.__origin[0],self.__origin[1],self.__origin[2]
        img.SetDimensions(nx, ny, nz)
        img.SetSpacing(sx, sy, sz)
        img.SetOrigin(ox, oy, oz)
        return img 

def __init__(self, extent=(10, 10, 10, 1), order='C',
                 spacing=(1.0, 1.0, 1.0), origin=(0.0, 0.0, 0.0),
                 dims=(0, 1, 2, 3), dt=1.0, points=False, **kwargs):
        FilterBase.__init__(self, nInputPorts=1, nOutputPorts=1,
                            inputType='vtkTable', outputType='vtkImageData', **kwargs)
        if len(extent) != 4:
            raise _helpers.PVGeoError('`extent` must be of length 4.')
        self.__extent = list(extent)
        self.__dims = list(dims) # these are indexes for the filter to use on the reshape.
        # NOTE: self.__dims[0] is the x axis index, etc., self.__dims[3] is the time axis
        self.__spacing = list(spacing) # image data spacing
        self.__origin = list(origin) # image data origin
        self.__order = list(order) # unpacking order: 'C' or 'F'
        self.__data = None # this is where we hold the data so entire filter does
        # not execute on every time step. Data will be a disctionary of 4D arrays
        # each 4D array will be in (nx, ny, nz, nt) shape
        self.__needToRun = True
        self.__timesteps = None
        self.__dt = dt
        # Optional parameter to switch between cell and point data
        self.__usePointData = points
        self.__needToUpdateOutput = True 

def RequestData(self, request, inInfo, outInfo):
        """Used by pipeline to generate output"""
        # Get output:
        output = vtk.vtkImageData.GetData(outInfo, 0)
        # Perform Read if needed
        if self.__raster is None:
            self._read_up_front()
        self._get_raw_data()
        self._build_image_data(output)
        # Now add the data based on what the user has selected
        for name, band in self.__raster.bands.items():
            data = band.data
            output.GetPointData().AddArray(interface.convert_array(data.flatten(), name=name))
        if self.__scheme is not None:
            colors = self.__raster.GetRGB(scheme=self.__scheme).reshape((-1,3))
            output.GetPointData().SetScalars(interface.convert_array(colors, name=self.__scheme))
            output.GetPointData().SetActiveScalars(self.__scheme)
        return 1 

def _get_volume_data(self):
        if self.object is None:
            return None
        elif isinstance(self.object, np.ndarray):
            return self._volume_from_array(self._subsample_array(self.object))
        else:
            available_serializer = [v for k, v in VTKVolume._serializers.items() if isinstance(self.object, k)]
            if not available_serializer:
                import vtk
                from vtk.util import numpy_support

                def volume_serializer(inst):
                    imageData = inst.object
                    array = numpy_support.vtk_to_numpy(imageData.GetPointData().GetScalars())
                    dims = imageData.GetDimensions()[::-1]
                    inst.spacing = imageData.GetSpacing()[::-1]
                    inst.origin = imageData.GetOrigin()
                    return inst._volume_from_array(inst._subsample_array(array.reshape(dims, order='C')))

                VTKVolume.register_serializer(vtk.vtkImageData, volume_serializer)
                serializer = volume_serializer
            else:
                serializer = available_serializer[0]
            return serializer(self) 

def numpy_to_image(numpy_array):
    """
    @brief Convert a numpy 2D or 3D array to a vtkImageData object
    @param numpy_array 2D or 3D numpy array containing image data
    @return vtkImageData with the numpy_array content
    """

    shape = numpy_array.shape
    if len(shape) < 2:
        raise Exception('numpy array must have dimensionality of at least 2')

    h, w = shape[0], shape[1]
    c = 1
    if len(shape) == 3:
        c = shape[2]

    # Reshape 2D image to 1D array suitable for conversion to a
    # vtkArray with numpy_support.numpy_to_vtk()
    linear_array = np.reshape(numpy_array, (w*h, c))
    vtk_array = numpy_to_vtk(linear_array)

    image = vtkImageData()
    image.SetDimensions(w, h, 1)
    image.AllocateScalars()
    image.GetPointData().GetScalars().DeepCopy(vtk_array)

    return image 

def test_simple_run(self):
        """`ExtractTopography`: Test extraction on simple data"""
        # Produce some input data
        data = vtk.vtkImageData()
        data.SetOrigin(0.0, 0.0, 0.0)
        data.SetSpacing(5.0, 5.0, 5.0)
        data.SetDimensions(20, 20, 20)
        x = y = np.linspace(0, 100, num=50, dtype=float)
        g = np.meshgrid(x, y)
        # Convert to XYZ points
        points = np.vstack(map(np.ravel, g)).T
        z = np.reshape(np.full(len(points), 55.0), (len(points), -1))
        #z = np.reshape(np.random.uniform(low=55.0, high=65.0, size=(len(points),)), (len(points), -1))
        points = np.concatenate((points, z), axis=1)
        topo = interface.points_to_poly_data(points)


        # # apply filter
        f = ExtractTopography()
        grd = f.apply(data, topo)

        # Test the output
        self.assertIsNotNone(grd)
        self.assertEqual(grd.GetDimensions(), data.GetDimensions())
        self.assertEqual(grd.GetSpacing(), data.GetSpacing())
        self.assertEqual(grd.GetOrigin(), data.GetOrigin())
        # Now check the active topo cell data?
        # TODO: implement
        active = dsa.WrapDataObject(grd).CellData['Extracted']
        for i in range(grd.GetNumberOfCells()):
            cell = grd.GetCell(i)
            bounds = cell.GetBounds()
            z = bounds[5]#(bounds[4]+bounds[5])/2.0
            if z <= 55.0:
                self.assertTrue(active[i])
            elif z > 55.0:
                self.assertFalse(active[i])
            else:
                self.fail(msg='Non-testable cell encountered')

        return 

def setUp(self):
        TestBase.setUp(self)
        # Create some input tables
        self.idi = vtk.vtkImageData()
        self.idi.SetDimensions(20, 2, 10)
        self.idi.SetSpacing(1, 1, 1)
        self.idi.SetOrigin(100, 100, 100)
        # Populate the tables
        self.n = 400
        self.title = 'Array 0'
        self.arr = np.random.random(self.n)
        self.idi.GetPointData().AddArray(interface.convert_array(self.arr, self.title))
        return 

def __init__(self, corner=1):
        FilterBase.__init__(self,
                            nInputPorts=1, inputType='vtkImageData',
                            nOutputPorts=1, outputType='vtkImageData')
        self.__corner = corner 

def set_origin(self, x0, y0, z0):
        """Set the origin of the output `vtkImageData`"""
        if self.__origin != [x0, y0, z0]:
            self.__origin = [x0, y0, z0]
            self.Modified() 

def _build_image_data(self, output):
        """Properly builds the output ``vtkImageData`` object"""
        if self.__raster is None:
            raise _helpers.PVGeoError('Raster invalid.')
        raster = self.__raster
        output.SetDimensions(raster.nsamps, raster.nlines, 1)
        output.SetSpacing(raster.pixel_size.x, raster.pixel_size.y, 1)
        corner = raster.global_metadata.projection_information.corner_point[0]
        output.SetOrigin(corner.x, corner.y, 0)
        return output


    #### Pipeline Methods #### 

def __init__(self, outputType='vtkImageData', **kwargs):
        DelimitedTextReader.__init__(self, outputType=outputType, **kwargs)
        # These are attributes the derived from file contents:
        self.set_delimiter(' ')
        self.__nx = None
        self.__ny = None
        self.__xo = None
        self.__yo = None
        self.__cellsize = None
        self.__data_name = kwargs.get('data_name', 'Data')
        self.NODATA_VALUE = -9999 

def perform_write_out(self, input_data_object, filename, object_name):
        """Writes an input ``vtkImageData`` object to a file"""
        img = input_data_object

        # Check dims: make sure 2D
        # TODO: handle any orientation
        nx, ny, nz = img.GetDimensions()
        if nx == ny == 1 and nz != 1:
            raise RuntimeError('Only 2D data on the XY plane is supported at this time.')

        ox, oy, oz = img.GetOrigin()
        dx, dy, dz = img.GetSpacing()

        # Get data ranges
        xmin, xmax, ymin, ymax, zmin, zmax = img.GetBounds()

        # Note user has to select a single array to save out
        field, name = self.__input_array[0], self.__input_array[1]
        vtkarr = _helpers.get_vtk_array(img, field, name)
        arr = interface.convert_array(vtkarr)
        dmin, dmax = arr.min(), arr.max()

        # arr = arr.reshape((nx, ny), order='F')

        meta = 'DSAA\n%d %d\n%f %f\n%f %f\n%f %f' % (nx, ny, xmin, xmax,
                                                     ymin, ymax, dmin, dmax)
        # Now write out the data!
        np.savetxt(filename, arr, header=meta, comments='', fmt=self.get_format())


        return 1 

def __init__(self):
        WriterBase.__init__(self, inputType='vtkImageData', ext='grd')
        self.__input_array = [None, None] 

def __init__(self, outputType='vtkImageData', **kwargs):
        ReaderBase.__init__(self, outputType=outputType, **kwargs)
        self.__grids = None
        self.__data_name = kwargs.get('data_name', 'Data') 

def to_vtk(self, output=None, z=0.0, dz=1.0, data_name='Data'):
        """Convert to a ``vtkImageData`` object"""
        self.mask()
        self.validate()
        if output is None:
            output = vtk.vtkImageData()
        # Build the data object
        output.SetOrigin(self.xll, self.yll, z)
        output.SetSpacing(self.dx, self.dy, dz)
        output.SetDimensions(self.nx, self.ny, 1)
        vtkarr = interface.convert_array(self.data, name=data_name)
        output.GetPointData().AddArray(vtkarr)
        return output 

def make_image_data():
    image_data = vtk.vtkImageData()
    image_data.SetDimensions(3, 4, 5)
    image_data.AllocateScalars(vtk.VTK_DOUBLE, 1)

    dims = image_data.GetDimensions()

    # Fill every entry of the image data with random double
    for z in range(dims[2]):
        for y in range(dims[1]):
            for x in range(dims[0]):
                image_data.SetScalarComponentFromDouble(x, y, z, 0, np.random.rand())
    return image_data 

def register_serializer(cls, class_type, serializer):
        """
        Register a seriliazer for a given type of class.
        A serializer is a function which take an instance of `class_type`
        (like a vtk.vtkImageData) as input and return a numpy array of the data
        """
        cls._serializers.update({class_type:serializer}) 

def applies(cls, obj):
        if ((isinstance(obj, np.ndarray) and obj.ndim == 3) or
            any([isinstance(obj, k) for k in cls._serializers.keys()])):
            return True
        elif 'vtk' not in sys.modules:
            return False
        else:
            import vtk
            return isinstance(obj, vtk.vtkImageData) 

def points(self, points):
        """Points cannot be set.

        This setter overrides the base class's setter to ensure a user does not
        attempt to set them. See https://github.com/pyvista/pyvista/issues/713.

        """
        raise AttributeError("The points cannot be set. The points of "
            "`UniformGrid`/`vtkImageData` are implicitly defined by the "
            "`origin`, `spacing`, and `dimensions` of the grid."
            ) 

def _clip_with_function(dataset, function, invert=True, value=0.0):
        """Clip using an implicit function (internal helper)."""
        if isinstance(dataset, vtk.vtkPolyData):
            alg = vtk.vtkClipPolyData()
        # elif isinstance(dataset, vtk.vtkImageData):
        #     alg = vtk.vtkClipVolume()
        #     alg.SetMixed3DCellGeneration(True)
        else:
            alg = vtk.vtkTableBasedClipDataSet()
        alg.SetInputDataObject(dataset) # Use the grid as the data we desire to cut
        alg.SetValue(value)
        alg.SetClipFunction(function) # the implicit function
        alg.SetInsideOut(invert) # invert the clip if needed
        alg.Update() # Perform the Cut
        return _get_output(alg) 

def __init__(self, *args, **kwargs):
        """Initialize the texture."""
        assert_empty_kwargs(**kwargs)

        if len(args) == 1:
            if isinstance(args[0], vtk.vtkTexture):
                self._from_texture(args[0])
            elif isinstance(args[0], np.ndarray):
                self._from_array(args[0])
            elif isinstance(args[0], vtk.vtkImageData):
                self._from_image_data(args[0])
            elif isinstance(args[0], str):
                self._from_texture(pyvista.read_texture(args[0]))
            else:
                raise TypeError('Table unable to be made from ({})'.format(type(args[0]))) 

def dump_image_data(dataset_dir, data_dir, dataset, color_array_info, root=None, compress=True):
    """Dump image data object to vtkjs."""
    if root is None:
        root = {}
    root['vtkClass'] = 'vtkImageData'
    container = root

    container['spacing'] = dataset.GetSpacing()
    container['origin'] = dataset.GetOrigin()
    container['extent'] = dataset.GetExtent()

    dump_all_arrays(dataset_dir, data_dir, dataset, container, compress)

    return root


# ----------------------------------------------------------------------------- 

def test_multi_block_init_vtk():
    multi = vtk.vtkMultiBlockDataSet()
    multi.SetBlock(0, vtk.vtkRectilinearGrid())
    multi.SetBlock(1, vtk.vtkStructuredGrid())
    multi = MultiBlock(multi)
    assert isinstance(multi, MultiBlock)
    assert multi.n_blocks == 2
    assert isinstance(multi.GetBlock(0), RectilinearGrid)
    assert isinstance(multi.GetBlock(1), StructuredGrid)
    multi = vtk.vtkMultiBlockDataSet()
    multi.SetBlock(0, vtk.vtkRectilinearGrid())
    multi.SetBlock(1, vtk.vtkStructuredGrid())
    multi = MultiBlock(multi, deep=True)
    assert isinstance(multi, MultiBlock)
    assert multi.n_blocks == 2
    assert isinstance(multi.GetBlock(0), RectilinearGrid)
    assert isinstance(multi.GetBlock(1), StructuredGrid)
    # Test nested structure
    multi = vtk.vtkMultiBlockDataSet()
    multi.SetBlock(0, vtk.vtkRectilinearGrid())
    multi.SetBlock(1, vtk.vtkImageData())
    nested = vtk.vtkMultiBlockDataSet()
    nested.SetBlock(0, vtk.vtkUnstructuredGrid())
    nested.SetBlock(1, vtk.vtkStructuredGrid())
    multi.SetBlock(2, nested)
    # Wrap the nested structure
    multi = MultiBlock(multi)
    assert isinstance(multi, MultiBlock)
    assert multi.n_blocks == 3
    assert isinstance(multi.GetBlock(0), RectilinearGrid)
    assert isinstance(multi.GetBlock(1), UniformGrid)
    assert isinstance(multi.GetBlock(2), MultiBlock) 

def clone(self):
        """
        Returns a Clone of an ImageDataBase instance
        """
        imagedata = vtkImageData()
        imagedata.DeepCopy(self._imagedata)
        clone = copy.copy(self)
        clone.set_imagedata(imagedata)
        clone.set_actor(actor_from_imagedata(self._imagedata))
        return clone 

def warp_by_scalar(dataset, scalars=None, factor=1.0, normal=None,
                       inplace=False, **kwargs):
        """Warp the dataset's points by a point data scalars array's values.

        This modifies point coordinates by moving points along point normals by
        the scalar amount times the scale factor.

        Parameters
        ----------
        scalars : str, optional
            Name of scalars to warp by. Defaults to currently active scalars.

        factor : float, optional
            A scaling factor to increase the scaling effect. Alias
            ``scale_factor`` also accepted - if present, overrides ``factor``.

        normal : np.array, list, tuple of length 3
            User specified normal. If given, data normals will be ignored and
            the given normal will be used to project the warp.

        inplace : bool
            If True, the points of the give dataset will be updated.

        """
        factor = kwargs.pop('scale_factor', factor)
        assert_empty_kwargs(**kwargs)
        if scalars is None:
            field, scalars = dataset.active_scalars_info
        arr, field = get_array(dataset, scalars, preference='point', info=True)
        if field != FieldAssociation.POINT:
            raise TypeError('Dataset can only by warped by a point data array.')
        # Run the algorithm
        alg = vtk.vtkWarpScalar()
        alg.SetInputDataObject(dataset)
        alg.SetInputArrayToProcess(0, 0, 0, field.value, scalars) # args: (idx, port, connection, field, name)
        alg.SetScaleFactor(factor)
        if normal is not None:
            alg.SetNormal(normal)
            alg.SetUseNormal(True)
        alg.Update()
        output = _get_output(alg)
        if inplace:
            if isinstance(dataset, (vtk.vtkImageData, vtk.vtkRectilinearGrid)):
                raise TypeError("This filter cannot be applied inplace for this mesh type.")
            dataset.overwrite(output)
            return
        return output 

def _translate(self, pdi, pdo):
        """Internal helper to translate the inputs origin"""
        if pdo is None:
            pdo = vtk.vtkImageData()

        [nx, ny, nz] = pdi.GetDimensions()
        [sx, sy, sz] = pdi.GetSpacing()
        [ox, oy, oz] = pdi.GetOrigin()

        pdo.DeepCopy(pdi)

        xx,yy,zz = 0.0,0.0,0.0

        if self.__corner == 1:
            # South East Bottom
            xx = ox - (nx-1)*sx
            yy = oy
            zz = oz
        elif self.__corner == 2:
            # North West Bottom
            xx = ox
            yy = oy - (ny-1)*sy
            zz = oz
        elif self.__corner == 3:
            # North East Bottom
            xx = ox - (nx-1)*sx
            yy = oy - (ny-1)*sy
            zz = oz
        elif self.__corner == 4:
            # South West Top
            xx = ox
            yy = oy
            zz = oz - (nz-1)*sz
        elif self.__corner == 5:
            # South East Top
            xx = ox - (nx-1)*sx
            yy = oy
            zz = oz - (nz-1)*sz
        elif self.__corner == 6:
            # North West Top
            xx = ox
            yy = oy - (ny-1)*sy
            zz = oz - (nz-1)*sz
        elif self.__corner == 7:
            # North East Top
            xx = ox - (nx-1)*sx
            yy = oy - (ny-1)*sy
            zz = oz - (nz-1)*sz

        pdo.SetOrigin(xx, yy, zz)

        return pdo 

def test(self):
        """`ReverseImageDataAxii`: check that it works on point and cell data"""
        # Create input vtkImageData:
        nx, ny, nz = 10, 11, 12
        arr = np.random.random((nz, ny, nx))
        arrCells = np.random.random((nz-1, ny-1, nx-1))
        image = vtk.vtkImageData()
        image.SetDimensions(nx, ny, nz)
        image.SetSpacing(2, 2, 2)
        image.SetOrigin(0, 0, 0)
        data = interface.convert_array(arr.flatten(), name='Data', deep=True)
        cellData = interface.convert_array(arrCells.flatten(), name='CellData', deep=True)
        image.GetPointData().AddArray(data)
        image.GetCellData().AddArray(cellData)
        # Now perfrom the reverse for only X:
        f = ReverseImageDataAxii()
        f.SetInputDataObject(image)
        f.set_flip_x(True)
        f.set_flip_y(False)
        f.set_flip_z(False)
        f.Update()
        ido = f.GetOutput()
        self.assertIsNotNone(ido)
        self.assertEqual(ido.GetNumberOfPoints(), len(arr.flatten()))
        # Check that x-axis was reversed
        wido = dsa.WrapDataObject(ido)
        test = wido.PointData['Data']
        testCells = wido.CellData['CellData']
        self.assertTrue(np.allclose(test, np.flip(arr, axis=2).flatten(), rtol=RTOL))
        self.assertTrue(np.allclose(testCells, np.flip(arrCells, axis=2).flatten(), rtol=RTOL))
        # Now perfrom the reverse for all axii:
        f.set_flip_x(True)
        f.set_flip_y(True)
        f.set_flip_z(True)
        f.Update()
        ido = f.GetOutput()
        self.assertIsNotNone(ido)
        self.assertEqual(ido.GetNumberOfPoints(), len(arr.flatten()))
        # Check that x-axis was reversed
        wido = dsa.WrapDataObject(ido)
        test = wido.PointData['Data']
        testCells = wido.CellData['CellData']
        tarr = np.flip(arr, axis=0)
        tarr = np.flip(tarr, axis=1)
        tarr = np.flip(tarr, axis=2)
        tarrCells = np.flip(arrCells, axis=0)
        tarrCells = np.flip(tarrCells, axis=1)
        tarrCells = np.flip(tarrCells, axis=2)
        self.assertTrue(np.allclose(test, tarr.flatten(), rtol=RTOL))
        self.assertTrue(np.allclose(testCells, tarrCells.flatten(), rtol=RTOL))


############################################################################### 

def to_itk_image(image_like):

    if isinstance(image_like, itk.Image):
        return image_like

    if is_arraylike(image_like):
        array = np.asarray(image_like)
        case_use_view = array.flags['OWNDATA']
        if have_dask and isinstance(image_like, dask.array.core.Array):
            case_use_view = False
        array = np.ascontiguousarray(array)
        if case_use_view:
            image_from_array = itk.image_view_from_array(array)
        else:
            image_from_array = itk.image_from_array(array)
        return image_from_array
    elif have_vtk and isinstance(image_like, vtk.vtkImageData):
        from vtk.util import numpy_support as vtk_numpy_support
        array = vtk_numpy_support.vtk_to_numpy(
            image_like.GetPointData().GetScalars())
        array.shape = tuple(image_like.GetDimensions())[::-1]
        image_from_array = itk.image_view_from_array(array)
        image_from_array.SetSpacing(image_like.GetSpacing())
        image_from_array.SetOrigin(image_like.GetOrigin())
        return image_from_array
    elif have_simpleitk and isinstance(image_like, sitk.Image):
        array = sitk.GetArrayViewFromImage(image_like)
        image_from_array = itk.image_view_from_array(array)
        image_from_array.SetSpacing(image_like.GetSpacing())
        image_from_array.SetOrigin(image_like.GetOrigin())
        direction = image_like.GetDirection()
        npdirection = np.asarray(direction)
        npdirection = np.reshape(npdirection, (-1, 3))
        itkdirection = itk.matrix_from_array(npdirection)
        image_from_array.SetDirection(itkdirection)
        return image_from_array
    elif have_imagej:
        import imglyb
        if isinstance(image_like,
                      imglyb.util.ReferenceGuardingRandomAccessibleInterval):
            array = imglyb.to_numpy(image_like)
            image_from_array = itk.image_view_from_array(array)
            return image_from_array
    elif isinstance(image_like, itk.ProcessObject):
        return itk.output(image_like)

    return None