Python vtk.util.numpy_support.vtk_to_numpy() Examples

def _read_stl(file_name):
    '''Internal function to read stl mesh files
    '''

    reader = vtk.vtkSTLReader()
    reader.SetFileName(file_name)
    reader.Update()

    mesh_data = PanelMesh(file_name)
    mesh_data.orig_type = 'Stereolithography (.stl)'
    mesh_data.num_faces = int(reader.GetOutput().GetNumberOfCells())
    mesh_data.num_points = mesh_data.num_faces * 3
    for i in range(mesh_data.num_faces):
        n = i*3
        mesh_data.faces.append(np.array([n,n+1,n+2,n+2]))
        mesh_data.points.append(np.array(vtk_to_numpy(reader.GetOutput().GetCell(i).GetPoints().GetData())))
    mesh_data.points = np.array(mesh_data.points).reshape([mesh_data.num_faces*3,3])

    return mesh_data 

def test_shallow_copy_back_propagation():
    """Test that the original data object's points get modified after a
    shallow copy.

    Reference: https://github.com/pyvista/pyvista/issues/375#issuecomment-531691483
    """
    # Case 1
    points = vtk.vtkPoints()
    points.InsertNextPoint(0.0, 0.0, 0.0)
    points.InsertNextPoint(1.0, 0.0, 0.0)
    points.InsertNextPoint(2.0, 0.0, 0.0)
    original = vtk.vtkPolyData()
    original.SetPoints(points)
    wrapped = pyvista.PolyData(original, deep=False)
    wrapped.points[:] = 2.8
    orig_points = vtk_to_numpy(original.GetPoints().GetData())
    assert np.allclose(orig_points, wrapped.points)
    # Case 2
    original = vtk.vtkPolyData()
    wrapped = pyvista.PolyData(original, deep=False)
    wrapped.points = np.random.rand(5, 3)
    orig_points = vtk_to_numpy(original.GetPoints().GetData())
    assert np.allclose(orig_points, wrapped.points) 

def _read_stl(file_name):
    '''Internal function to read stl mesh files
    '''

    reader = vtk.vtkSTLReader()
    reader.SetFileName(file_name)
    reader.Update()

    mesh_data = PanelMesh(file_name)
    mesh_data.orig_type = 'Stereolithography (.stl)'
    mesh_data.num_faces = int(reader.GetOutput().GetNumberOfCells())
    mesh_data.num_points = mesh_data.num_faces * 3
    for i in range(mesh_data.num_faces):
        n = i*3
        mesh_data.faces.append(np.array([n,n+1,n+2,n+2]))
        mesh_data.points.append(np.array(vtk_to_numpy(reader.GetOutput().GetCell(i).GetPoints().GetData())))
    mesh_data.points = np.array(mesh_data.points).reshape([mesh_data.num_faces*3,3])

    return mesh_data 

def _read_stl(file_name):
    '''Internal function to read stl mesh files
    '''

    reader = vtk.vtkSTLReader()
    reader.SetFileName(file_name)
    reader.Update()

    mesh_data = PanelMesh(file_name)
    mesh_data.orig_type = 'Stereolithography (.stl)'
    mesh_data.num_faces = int(reader.GetOutput().GetNumberOfCells())
    mesh_data.num_points = mesh_data.num_faces * 3
    for i in range(mesh_data.num_faces):
        n = i*3
        mesh_data.faces.append(np.array([n,n+1,n+2,n+2]))
        mesh_data.points.append(np.array(vtk_to_numpy(reader.GetOutput().GetCell(i).GetPoints().GetData())))
    mesh_data.points = np.array(mesh_data.points).reshape([mesh_data.num_faces*3,3])

    return mesh_data 

def _read_stl(file_name):
    '''Internal function to read stl mesh files
    '''

    reader = vtk.vtkSTLReader()
    reader.SetFileName(file_name)
    reader.Update()

    mesh_data = PanelMesh(file_name)
    mesh_data.orig_type = 'Stereolithography (.stl)'
    mesh_data.num_faces = int(reader.GetOutput().GetNumberOfCells())
    mesh_data.num_points = mesh_data.num_faces * 3
    for i in range(mesh_data.num_faces):
        n = i*3
        mesh_data.faces.append(np.array([n,n+1,n+2,n+2]))
        mesh_data.points.append(np.array(vtk_to_numpy(reader.GetOutput().GetCell(i).GetPoints().GetData())))
    mesh_data.points = np.array(mesh_data.points).reshape([mesh_data.num_faces*3,3])

    return mesh_data 

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 test_VTK_object_conversion(self):
            mesh = self.mesh
            vec = np.arange(mesh.nC)
            models = {'arange': vec}

            vtkObj = mesh.to_vtk(models)

            self.assertEqual(mesh.nC, vtkObj.GetNumberOfCells())
            # TODO: this is actually different?: self.assertEqual(mesh.nN, vtkObj.GetNumberOfPoints())
            # Remember that the tree vtk conversion adds two arrays
            self.assertEqual(len(models.keys())+2, vtkObj.GetCellData().GetNumberOfArrays())
            bnds = vtkObj.GetBounds()
            self.assertEqual(mesh.x0[0], bnds[0])
            self.assertEqual(mesh.x0[1], bnds[2])
            self.assertEqual(mesh.x0[2], bnds[4])

            names = list(models.keys())
            for name in names:
                arr = nps.vtk_to_numpy(vtkObj.GetCellData().GetArray(name))
                arr = arr.flatten(order='F')
                self.assertTrue(np.allclose(models[name], arr)) 

def test_VTK_object_conversion(self):
            mesh = self.mesh
            vec = np.arange(mesh.nC)
            models = {'arange': vec}

            vtkObj = mesh.to_vtk(models)

            self.assertEqual(mesh.nC, vtkObj.GetNumberOfCells())
            self.assertEqual(mesh.nN, vtkObj.GetNumberOfPoints())
            self.assertEqual(len(models.keys()), vtkObj.GetCellData().GetNumberOfArrays())
            bnds = vtkObj.GetBounds()
            self.assertEqual(mesh.x0[0], bnds[0])
            self.assertEqual(mesh.x0[1], bnds[2])
            self.assertEqual(mesh.x0[2], bnds[4])

            names = list(models.keys())
            for i in range(vtkObj.GetCellData().GetNumberOfArrays()):
                name = names[i]
                self.assertEqual(name, vtkObj.GetCellData().GetArrayName(i))
                arr = nps.vtk_to_numpy(vtkObj.GetCellData().GetArray(i))
                arr = arr.flatten(order='F')
                self.assertTrue(np.allclose(models[name], arr)) 

def test_VTK_object_conversion(self):
            mesh = self.mesh
            vec = np.arange(mesh.nC)
            models = {'arange': vec}

            vtkObj = mesh.to_vtk(models)

            self.assertEqual(mesh.nC, vtkObj.GetNumberOfCells())
            self.assertEqual(mesh.nN, vtkObj.GetNumberOfPoints())
            self.assertEqual(len(models.keys()), vtkObj.GetCellData().GetNumberOfArrays())
            bnds = vtkObj.GetBounds()
            self.assertEqual(mesh.x0[0], bnds[0])
            self.assertEqual(mesh.x0[1], bnds[2])
            self.assertEqual(mesh.x0[2], bnds[4])

            for i in range(vtkObj.GetCellData().GetNumberOfArrays()):
                name = list(models.keys())[i]
                self.assertEqual(name, vtkObj.GetCellData().GetArrayName(i))
                arr = nps.vtk_to_numpy(vtkObj.GetCellData().GetArray(i))
                arr = arr.flatten(order='F')
                self.assertTrue(np.allclose(models[name], arr)) 

def _check_data(self, grid, data):
        arr = nps.vtk_to_numpy(grid.GetCellData().GetArray(0))
        self.assertTrue(np.allclose(data, arr, rtol=RTOL))
        return 

def readData(filename, x = all, y = all, z = all):
    """Read data from raw/mhd image
    
    Arguments:
        filename (str): file name as regular expression
        x,y,z (tuple): data range specifications
    
    Returns:
        array: image data
    """   
    
    imr = vtk.vtkMetaImageReader()
    imr.SetFileName(filename);
    imr.Update()
    
    im = imr.GetOutput()
    dims = im.GetDimensions()
    
    print dims    
    
    sc = im.GetPointData().GetScalars()
    img = vtk_to_numpy(sc)
    #print img.shape
    
    dims = list(dims);
    dims[0:3] = [dims[2], dims[1], dims[0]];
    
    imgs = list(img.shape);
    if len(imgs) > 1:
        imgs.pop(0);
        dims = dims + imgs;
    
    img = img.reshape(dims)
    #img = img.transpose([1,2,0]);
    tp = [2,1,0];
    tp = tp + [i for i in range(3, len(dims))];
    img = img.transpose(tp);
    
    return io.dataToRange(img, x = x, y = y, z = z); 

def test_numpy_to_idarr_bool():
    mask = np.ones(10, np.bool_)
    idarr = pyvista.utilities.cells.numpy_to_idarr(mask)
    assert np.allclose(mask.nonzero()[0], vtk_to_numpy(idarr)) 

def cells(self):
        """Return a numpy array of the cells."""
        return vtk_to_numpy(self.GetData()).ravel() 

def _vtk_to_vtkjs(data_array):
    from vtk.util.numpy_support import vtk_to_numpy
    # From vtkType.h
    _vtk_data_type_to_vtkjs_type = {
        2: 'Int8Array',
        15: 'Int8Array',
        3: 'Uint8Array',
        4: 'Int16Array',
        5: 'Uint16Array',
        6: 'Int32Array',
        7: 'Uint32Array',
        8: 'BigInt64Array',
        9: 'BigUint64Array',
        10: 'Float32Array',
        11: 'Float64Array',
        16: 'BigInt64Array',
        17: 'BigUint64Array',
    }
    vtk_data_type = data_array.GetDataType()
    data_type = _vtk_data_type_to_vtkjs_type[vtk_data_type]
    numpy_array = vtk_to_numpy(data_array)
    if vtk_data_type == 8 or vtk_data_type == 16:
        ii32 = np.iinfo(np.int32)
        value_range = data_array.GetValueRange()
        if value_range[0] < ii32.min or value_range[1] > ii32.max:
            raise ValueError(
                '64 integers are not supported yet by WebGL / vtk.js')
        numpy_array = numpy_array.astype(np.int32)
        data_type = 'Int32Array'
    elif vtk_data_type == 9 or vtk_data_type == 17:
        ui32 = np.iinfo(np.uint32)
        value_range = data_array.GetValueRange()
        if value_range[0] < ui32.min or value_range[1] > ui32.max:
            raise ValueError(
                '64 integers are not supported by WebGL / vtk.js yet')
        numpy_array = numpy_array.astype(np.uint32)
        data_type = 'Uint32Array'

    return data_type, numpy_array 

def getNumpyFromVtk(dataObj, arrayName='Points', arrayType='points'):
    assert arrayType in ('points', 'cells')

    if arrayName == 'Points':
        vtkArray = dataObj.GetPoints().GetData()
    elif arrayType == 'points':
        vtkArray = dataObj.GetPointData().GetArray(arrayName)
    else:
        vtkArray = dataObj.GetCellData().GetArray(arrayName)

    if not vtkArray:
        raise KeyError('Array not found')

    return numpy_support.vtk_to_numpy(vtkArray) 

def getNumpyImageFromVtk(vtkimg, flip=True):
    img = numpy_support.vtk_to_numpy(vtkimg.GetPointData().GetScalars())
    w, h, _ = vtkimg.GetDimensions()
    img.shape = (h, w, -1)
    if img.shape[2] == 1:
        img.shape = img.shape[:2]
    if flip:
        img = np.flipud(img)
    return img 

def reshapeArrs(self, mesh):
            for i in range(self.nt):
                ind_reorder = nps.vtk_to_numpy(
                    mesh.GetCellData().GetArray('index_cell_corner'))
                self.arrs[i] = self.arrs[i][ind_reorder] 

def vtk_to_tensor_mesh(TensorMesh, vtrGrid):
        """Converts a ``vtkRectilinearGrid`` or :class:`pyvista.RectilinearGrid`
        to a :class:`discretize.TensorMesh` object.

        """
        # Sort information
        hx = np.abs(np.diff(_nps.vtk_to_numpy(vtrGrid.GetXCoordinates())))
        xR = _nps.vtk_to_numpy(vtrGrid.GetXCoordinates())[0]
        hy = np.abs(np.diff(_nps.vtk_to_numpy(vtrGrid.GetYCoordinates())))
        yR = _nps.vtk_to_numpy(vtrGrid.GetYCoordinates())[0]
        zD = np.diff(_nps.vtk_to_numpy(vtrGrid.GetZCoordinates()))
        # Check the direction of hz
        if np.all(zD < 0):
            hz = np.abs(zD[::-1])
            zR = _nps.vtk_to_numpy(vtrGrid.GetZCoordinates())[-1]
        else:
            hz = np.abs(zD)
            zR = _nps.vtk_to_numpy(vtrGrid.GetZCoordinates())[0]
        x0 = np.array([xR, yR, zR])

        # Make the object
        tensMsh = TensorMesh([hx, hy, hz], x0=x0)

        # Grap the models
        models = {}
        for i in np.arange(vtrGrid.GetCellData().GetNumberOfArrays()):
            modelName = vtrGrid.GetCellData().GetArrayName(i)
            if np.all(zD < 0):
                modFlip = _nps.vtk_to_numpy(vtrGrid.GetCellData().GetArray(i))
                tM = tensMsh.r(modFlip, 'CC', 'CC', 'M')
                modArr = tensMsh.r(tM[:, :, ::-1], 'CC', 'CC', 'V')
            else:
                modArr = _nps.vtk_to_numpy(vtrGrid.GetCellData().GetArray(i))
            models[modelName] = modArr

        # Return the data
        return tensMsh, models 

def _check_data(self, table, data):
        arr = nps.vtk_to_numpy(table.GetColumn(0))
        self.assertTrue(np.allclose(data, arr, rtol=RTOL))
        return arr

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

def test_data_fidelity(self):
        """`XYZTextReader`: check data fidelity"""
        titles = self.header.split(', ')
        for i in range(self.ncols):
            arr = nps.vtk_to_numpy(self.TABLE.GetColumnByName(titles[i]))
            self.assertTrue(np.allclose(self.data[:,i], arr, rtol=RTOL))
        return 

def test_data(self,):
        """`SGeMSGridReader`: data fidelity"""
        for i, title in enumerate(self.titles):
            arr = nps.vtk_to_numpy(self.GRID.GetCellData().GetArray(i))
            self.assertTrue(np.allclose(self.data[:,i], arr, rtol=RTOL))
        return 

def test_data_fidelity(self):
        """`GSLibReader`: data fidelity"""
        for i, title in enumerate(self.titles):
            arr = nps.vtk_to_numpy(self.TABLE.GetColumn(i))
            self.assertTrue(np.allclose(self.data[:,i], arr, rtol=RTOL))
        return 

def place_model_on_octree_mesh(mesh, model, data_name='Data'):
        """Places model data onto a mesh. This is for the UBC Grid data reaers
        to associate model data with the mesh grid.

        Args:
            mesh (vtkUnstructuredGrid): The ``vtkUnstructuredGrid`` that is the
                mesh to place the model data upon. Needs to have been read in by ubcOcTree
            model (np.ndarray): A NumPy float array that holds all of the data
                to place inside of the mesh's cells.
            data_name (str): The name of the model data array once placed on the
                ``vtkUnstructuredGrid``.

        Return:
            vtkUnstructuredGrid:
                The input ``vtkUnstructuredGrid`` with model data appended.
        """
        if isinstance(model, dict):
            for key in model.keys():
                mesh = OcTreeReader.place_model_on_octree_mesh(mesh, model[key], data_name=key)
            return mesh
        # Make sure this model file fits the dimensions of the mesh
        numCells = mesh.GetNumberOfCells()
        if (numCells < len(model)):
            raise _helpers.PVGeoError('This model file has more data than the given mesh has cells to hold.')
        elif (numCells > len(model)):
            raise _helpers.PVGeoError('This model file does not have enough data to fill the given mesh\'s cells.')

        # This is absolutely crucial!
        # Do not play with unless you know what you are doing!
        # Also note that this assumes ``discretize`` handles addin this array
        ind_reorder = nps.vtk_to_numpy(
            mesh.GetCellData().GetArray('index_cell_corner'))

        model = model[ind_reorder]

        # Convert data to VTK data structure and append to output
        c = interface.convert_array(model, name=data_name, deep=True)
        # THIS IS CELL DATA! Add the model data to CELL data:
        mesh.GetCellData().AddArray(c)
        return mesh 

def y(self):
        """Get the coordinates along the Y-direction."""
        return vtk_to_numpy(self.GetYCoordinates()) 

def z(self):
        """Get the coordinates along the Z-direction."""
        return vtk_to_numpy(self.GetZCoordinates()) 

def vtkPolyData_to_tracts(polydata):

	result = {}
	result['lines'] = ns.vtk_to_numpy(polydata.GetLines().GetData())
	result['points'] = ns.vtk_to_numpy(polydata.GetPoints().GetData())
	result['numberOfLines'] = polydata.GetNumberOfLines()

	data = {}
	if polydata.GetPointData().GetScalars():
		data['ActiveScalars'] = polydata.GetPointData().GetScalars().GetName()
		result['Scalars'] = polydata.GetPointData().GetScalars()
	if polydata.GetPointData().GetVectors():
		data['ActiveVectors'] = polydata.GetPointData().GetVectors().GetName()
	if polydata.GetPointData().GetTensors():
		data['ActiveTensors'] = polydata.GetPointData().GetTensors().GetName()

	for i in xrange(polydata.GetPointData().GetNumberOfArrays()):
		array = polydata.GetPointData().GetArray(i)
		np_array = ns.vtk_to_numpy(array)
		if np_array.ndim == 1:
			np_array = np_array.reshape(len(np_array), 1)
		data[polydata.GetPointData().GetArrayName(i)] = np_array

	result['pointData'] = data

	tracts, data = vtkPolyData_dictionary_to_tracts_and_data(result)
	return tracts, data 

def curvature(poly_data, curv_type='mean'):
        """Return the pointwise curvature of a mesh.

        Parameters
        ----------
        mesh : vtk.polydata
            vtk polydata mesh

        curvature string, optional
            One of the following strings
            Mean
            Gaussian
            Maximum
            Minimum

        Return
        ------
        curvature : np.ndarray
            Curvature values

        """
        curv_type = curv_type.lower()

        # Create curve filter and compute curvature
        curvefilter = vtk.vtkCurvatures()
        curvefilter.SetInputData(poly_data)
        if curv_type == 'mean':
            curvefilter.SetCurvatureTypeToMean()
        elif curv_type == 'gaussian':
            curvefilter.SetCurvatureTypeToGaussian()
        elif curv_type == 'maximum':
            curvefilter.SetCurvatureTypeToMaximum()
        elif curv_type == 'minimum':
            curvefilter.SetCurvatureTypeToMinimum()
        else:
            raise ValueError('Curv_Type must be either "Mean", '
                             '"Gaussian", "Maximum", or "Minimum"')
        curvefilter.Update()

        # Compute and return curvature
        curv = _get_output(curvefilter)
        return vtk_to_numpy(curv.GetPointData().GetScalars()) 

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 

def convert_array(arr, name=None, deep=0, array_type=None):
    """Convert a NumPy array to a vtkDataArray or vice versa.

    Parameters
    -----------
    arr : ndarray or vtkDataArry
        A numpy array or vtkDataArry to convert
    name : str
        The name of the data array for VTK
    deep : bool
        if input is numpy array then deep copy values

    Return
    ------
    vtkDataArray, ndarray, or DataFrame:
        the converted array (if input is a NumPy ndaray then returns
        ``vtkDataArray`` or is input is ``vtkDataArray`` then returns NumPy
        ``ndarray``). If pdf==True and the input is ``vtkDataArry``,
        return a pandas DataFrame.

    """
    if arr is None:
        return
    if isinstance(arr, np.ndarray):
        if arr.dtype is np.dtype('O'):
            arr = arr.astype('|S')
        arr = np.ascontiguousarray(arr)
        if arr.dtype.type in (np.str_, np.bytes_):
            # This handles strings
            vtk_data = convert_string_array(arr)
        else:
            # This will handle numerical data
            arr = np.ascontiguousarray(arr)
            vtk_data = nps.numpy_to_vtk(num_array=arr, deep=deep, array_type=array_type)

        if isinstance(name, str):
            vtk_data.SetName(name)
        return vtk_data
    # Otherwise input must be a vtkDataArray
    if not isinstance(arr, (vtk.vtkDataArray, vtk.vtkBitArray, vtk.vtkStringArray)):
        raise TypeError('Invalid input array type ({}).'.format(type(arr)))
    # Handle booleans
    if isinstance(arr, vtk.vtkBitArray):
        arr = vtk_bit_array_to_char(arr)
    # Handle string arrays
    if isinstance(arr, vtk.vtkStringArray):
        return convert_string_array(arr)
    # Convert from vtkDataArry to NumPy
    return nps.vtk_to_numpy(arr)