将3D numpy导出为VTK文件,以便在Paraview / Mayavi中查看

问题描述 投票:6回答:5

对于那些想要将简单的3D numpy数组(以及轴)导出到.vtk(或.vtr)文件以进行后期处理并在Paraview或Mayavi中显示的人,有一个名为PyEVTK的小模块就是这样做的。该模块支持结构化和非结构化数据等。不幸的是,即使代码在基于unix的系统中工作正常,我也无法在任何简单地使事情变得复杂的Windows安装上工作(不断崩溃)。我联系了开发人员,但他的建议没有用

因此我的问题是:如何使用from vtk.util import numpy_support函数将3D数组(函数本身不支持3D数组)导出到.vtk文件?有没有一种简单的方法可以在不创建vtkDatasets等的情况下完成它?

非常感谢!

numpy vtk paraview
5个回答
7
投票

这是永远的,我完全忘记了这个问题,但我最终搞清楚了。我在我的博客(PyScience)上写了一篇关于它的文章,提供了如何在NumPy和VTK之间进行转换的教程。如果感兴趣,请看看:

pyscience.wordpress.com/2014/09/06/numpy-to-vtk-converting-your-numpy-arrays-to-vtk-arrays-and-files/


4
投票

这不是你的问题的直接答案,但如果你有tvtk(如果你有mayavi,你应该拥有它),你可以用它来写你的数据到vtk格式。 (见:http://code.enthought.com/projects/files/ETS3_API/enthought.tvtk.misc.html

它不使用PyEVTK,它支持广泛的数据源(不仅仅是结构化和非结构化网格),因此它可能适用于其他事情不适用的地方。

作为一个简单的例子(Mayavi的mlab界面可以减少冗长,特别是如果你已经在使用它。):

import numpy as np
from enthought.tvtk.api import tvtk, write_data

data = np.random.random((10,10,10))

grid = tvtk.ImageData(spacing=(10, 5, -10), origin=(100, 350, 200), 
                      dimensions=data.shape)
grid.point_data.scalars = np.ravel(order='F')
grid.point_data.scalars.name = 'Test Data'

# Writes legacy ".vtk" format if filename ends with "vtk", otherwise
# this will write data using the newer xml-based format.
write_data(grid, 'test.vtk')

以及输出文件的一部分:

# vtk DataFile Version 3.0
vtk output
ASCII
DATASET STRUCTURED_POINTS
DIMENSIONS 10 10 10
SPACING 10 5 -10
ORIGIN 100 350 200
POINT_DATA 1000
SCALARS Test%20Data double
LOOKUP_TABLE default
0.598189 0.228948 0.346975 0.948916 0.0109774 0.30281 0.643976 0.17398 0.374673 
0.295613 0.664072 0.307974 0.802966 0.836823 0.827732 0.895217 0.104437 0.292796 
0.604939 0.96141 0.0837524 0.498616 0.608173 0.446545 0.364019 0.222914 0.514992 
...
...

2
投票

我知道这有点晚了,我很喜欢你的教程@ somada141。这应该也有效。

def numpy2VTK(img, spacing=[1.0, 1.0, 1.0]):
 # evolved from code from Stou S.,
 # on http://www.siafoo.net/snippet/314
 # This function, as the name suggests, converts numpy array to VTK
 importer = vtk.vtkImageImport()

 img_data = img.astype('uint8')
 img_string = img_data.tostring()  # type short
 dim = img.shape

 importer.CopyImportVoidPointer(img_string, len(img_string))
 importer.SetDataScalarType(VTK_UNSIGNED_CHAR)
 importer.SetNumberOfScalarComponents(1)

 extent = importer.GetDataExtent()
 importer.SetDataExtent(extent[0], extent[0] + dim[2] - 1,
                       extent[2], extent[2] + dim[1] - 1,
                       extent[4], extent[4] + dim[0] - 1)
 importer.SetWholeExtent(extent[0], extent[0] + dim[2] - 1,
                        extent[2], extent[2] + dim[1] - 1,
                        extent[4], extent[4] + dim[0] - 1)

 importer.SetDataSpacing(spacing[0], spacing[1], spacing[2])
 importer.SetDataOrigin(0, 0, 0)


 return importer

希望能帮助到你!


1
投票

Mayavi的TVTK有一种编写vtk文件的漂亮方式。这是我在@Joe和tvtk文档之后为自己编写的测试示例。它对evtk的优势在于对ascii和html的支持。希望它能帮助其他人。

from tvtk.api import tvtk, write_data
import numpy as np

#data = np.random.random((3, 3, 3))
#
#i = tvtk.ImageData(spacing=(1, 1, 1), origin=(0, 0, 0))
#i.point_data.scalars = data.ravel()
#i.point_data.scalars.name = 'scalars'
#i.dimensions = data.shape
#
#w = tvtk.XMLImageDataWriter(input=i, file_name='spoints3d.vti')
#w.write()

points = np.array([[0,0,0], [1,0,0], [1,1,0], [0,1,0]], 'f')
(n1, n2)  = points.shape
poly_edge = np.array([[0,1,2,3]])

print n1, n2
## Scalar Data
#temperature = np.array([10., 20., 30., 40.])
#pressure = np.random.rand(n1)
#
## Vector Data
#velocity = np.random.rand(n1,n2)
#force     = np.random.rand(n1,n2)
#
##Tensor Data with 
comp = 5
stress = np.random.rand(n1,comp)
#
#print stress.shape
## The TVTK dataset.
mesh = tvtk.PolyData(points=points, polys=poly_edge)
#
## Data 0 # scalar data
#mesh.point_data.scalars = temperature
#mesh.point_data.scalars.name = 'Temperature'
#
## Data 1 # additional scalar data
#mesh.point_data.add_array(pressure)
#mesh.point_data.get_array(1).name = 'Pressure'
#mesh.update()
#
## Data 2 # Vector data
#mesh.point_data.vectors = velocity
#mesh.point_data.vectors.name = 'Velocity'
#mesh.update()
#
## Data 3 additional vector data
#mesh.point_data.add_array( force)
#mesh.point_data.get_array(3).name = 'Force'
#mesh.update()

mesh.point_data.tensors = stress
mesh.point_data.tensors.name = 'Stress'

# Data 4 additional tensor Data
#mesh.point_data.add_array(stress)
#mesh.point_data.get_array(4).name = 'Stress'
#mesh.update()

write_data(mesh, 'polydata.vtk')

# XML format 
# Method 1
#write_data(mesh, 'polydata')

# Method 2
#w = tvtk.XMLPolyDataWriter(input=mesh, file_name='polydata.vtk')
#w.write()

1
投票

这是一个SimpleITK版本,功能load_itk取自here

import SimpleITK as sitk
import numpy as np


if len(sys.argv)<3:
    print('Wrong number of arguments.', file=sys.stderr)
    print('Usage: ' + __file__ + ' input_sitk_file' + ' output_sitk_file', file=sys.stderr)
    sys.exit(1)


def quick_read(filename):
    # Read image information without reading the bulk data.
    file_reader = sitk.ImageFileReader()
    file_reader.SetFileName(filename)
    file_reader.ReadImageInformation()
    print('image size: {0}\nimage spacing: {1}'.format(file_reader.GetSize(), file_reader.GetSpacing()))
    # Some files have a rich meta-data dictionary (e.g. DICOM)
    for key in file_reader.GetMetaDataKeys():
        print(key + ': ' + file_reader.GetMetaData(key))



def load_itk(filename):
    # Reads the image using SimpleITK
    itkimage = sitk.ReadImage(filename)
    # Convert the image to a  numpy array first and then shuffle the dimensions to get axis in the order z,y,x
    data = sitk.GetArrayFromImage(itkimage)
    # Read the origin of the ct_scan, will be used to convert the coordinates from world to voxel and vice versa.
    origin = np.array(list(reversed(itkimage.GetOrigin())))
    # Read the spacing along each dimension
    spacing = np.array(list(reversed(itkimage.GetSpacing())))
    return data, origin, spacing


def convert(data, output_filename):
    image = sitk.GetImageFromArray(data)
    writer = sitk.ImageFileWriter()
    writer.SetFileName(output_filename)
    writer.Execute(image)


def wait():
    print('Press Enter to load & convert or exit using Ctrl+C')
    input()


quick_read(sys.argv[1])
print('-'*20)
wait()
data, origin, spacing = load_itk(sys.argv[1])
convert(sys.argv[2])
© www.soinside.com 2019 - 2024. All rights reserved.