使用python-vtk将多个遗留ASCII .vtk文件合并为一个文件
在我的计算过程中，我保存了一堆与某个时间点相关的 .vtk 文件。每个文件描述了一个多面体（这是我定义的一个 C++ 类），使用的是 POLYDATA ASCII 文件格式。多面体是通过多面体类的一个简单成员函数来写入的。
为了避免为我需要可视化的多面体集合定义一个全新的类，从而搞乱我的 C++ 代码，我想把多个 .vtk 文件合并成一个 .vtk 文件。
使用 python-vtk 时我遇到了一些问题：
from vtk import *
reader = vtkPolyDataReader()
reader.SetFileName("file1.vtk")
reader.Update()
polyData1 = reader.GetOutput()
reader.SetFileName('file2.vtk')
reader.Update()
polyData2 = reader.GetOutput()
# Expand the output points
points1 = polyData1.GetPoints()
points2 = polyData2.GetPoints()
insertPosition = points1.GetNumberOfPoints()
for i in xrange(points2.GetNumberOfPoints()):
insertPoint = points2.GetPoint(i)
points1.InsertPoint(insertPosition,
insertPoint[0], insertPoint[1], insertPoint[2])
insertPosition += 1
print points1.GetNumberOfPoints()
# Change the cell ids of every cell in the polydata2 to correspond with
# the new points (appended point array)
increment = points1.GetNumberOfPoints();
for i in xrange(polyData2.GetNumberOfCells()):
cell = polyData2.GetCell(i)
cellIds = cell.GetPointIds()
for j in xrange(cellIds.GetNumberOfIds()):
oldId = cellIds.GetId(j)
cellIds.SetId(j, oldId + increment)
polyData1.Allocate(polyData1.GetNumberOfCells(), 1)
for i in xrange(polyData2.GetNumberOfCells()):
cell = polyData2.GetCell(i)
polyData1.InsertNextCell(cell.GetCellType(), cell.GetPointIds())
writer = vtkPolyDataWriter()
writer.SetFileName("output.vtk")
writer.SetInput(polyData1)
writer.Write()
这样做会导致重复的点，这没关系。问题是这个脚本在以下 .vtk 文件上执行：
文件1：
# vtk DataFile Version 2.0
surface written 2011-12-19T15:30:18
ASCII
DATASET POLYDATA
POINTS 8 float
0.48999999999999999112 0.4000000000000000222 0.5999999999999999778
0.48999999999999999112 0.5 0.5999999999999999778
0.48999999999999999112 0.5 0.69999999999999995559
0.48999999999999999112 0.4000000000000000222 0.69999999999999995559
0.5 0.5 0.5999999999999999778
0.5 0.5 0.69999999999999995559
0.5 0.4000000000000000222 0.69999999999999995559
0.5 0.4000000000000000222 0.5999999999999999778
POLYGONS 6 30
4 0 1 2 3
4 4 5 6 7
4 4 1 2 5
4 6 5 2 3
4 7 0 1 4
4 0 7 6 3
CELL_DATA 6
FIELD attributes 1
zone 1 6 float
1 1 1 1 1 1
文件2：
# vtk DataFile Version 2.0
surface written 2011-12-19T15:30:18
ASCII
DATASET POLYDATA
POINTS 8 float
0.58999999999999996891 0.5999999999999999778 0.5
0.58999999999999996891 0.69999999999999995559 0.5
0.58999999999999996891 0.69999999999999995559 0.5999999999999999778
0.58999999999999996891 0.5999999999999999778 0.5999999999999999778
0.5999999999999999778 0.69999999999999995559 0.5
0.5999999999999999778 0.69999999999999995559 0.5999999999999999778
0.5999999999999999778 0.5999999999999999778 0.5999999999999999778
0.5999999999999999778 0.5999999999999999778 0.5
POLYGONS 6 30
4 0 1 2 3
4 4 5 6 7
4 4 1 2 5
4 6 5 2 3
4 7 0 1 4
4 0 7 6 3
CELL_DATA 6
FIELD attributes 1
zone 1 6 float
1 1 1 1 1 1
结果是点的坐标变成了 8 个，而单元（面）根本没有被添加。
是我一个人遇到这个问题，还是说 python 对 vtkArray 和类似 vtkObjects 的封装没有迭代的选项呢？
2 个回答
2
为了这个目的，我放弃了python-vtk这个接口。这里有一个简单又快速写成的类，可以完成这个任务，希望对某些人有帮助：
class polyDataVtk(object):
"""Class representing the polydata vtk information stored in legacy ASCII .vtk POLYDATA files."""
def __init__(self, fileName = None):
self.__points = []
self.__polygons = []
if fileName is not None:
self.__fileName = fileName
def parse(self, fileName):
"""Parse the POLYDATA information from a .vtk file and append the data to the object.
Does not check for the file consistency."""
file = open(fileName, 'r')
# Use local data first.
points = []
polygons = []
line = ""
while(True):
line = file.readline()
if 'POINTS' in line:
break
nPoints = 0
if (line == ""):
print "No POINTS defined in the .vtk file"
return
# Set the number of points
nPoints = int(line.split()[1])
# Append the numbers.
for i in xrange(nPoints):
points.append(map(lambda x : float(x), file.readline().split()))
# Append polygons.
line = ""
while(True):
line = file.readline()
if 'POLYGONS' in line:
break
if (line == ""):
print "No POLYGONS defined in the .vtk file"
return
# Set the number of polygons.
nPolygons = int(line.split()[1])
# Read the polygons.
for i in xrange(nPolygons):
line = file.readline()
polygons.append(map(lambda x : int(x) + len(self.__points), line.split())[1:])
# File parsed without a problem.
self.__points.extend(points)
self.__polygons.extend(polygons)
file.close()
def write(self,fileName=None, append=False):
# Overwrite the filename provided to the constructor.
if fileName is not None:
self.__fileName = fileName
# No fileName is provided by the constructor or the write method.
if self.__fileName is None:
self.__fileName = "result.vtk"
# Append or overwrite?
if append:
file = open(self.__fileName, 'a')
else:
file = open(self.__fileName, 'w')
file.writelines("# vtk DataFile Version 2.0\n")
file.writelines("appended vtk files\n")
file.writelines("ASCII\n")
file.writelines("DATASET POLYDATA\n")
file.writelines("POINTS %d float \n" % len(self.__points))
for point in self.__points:
file.writelines("%.10f %.10f %.10f\n" % (point[0], point[1], point[2]))
size = 0
for polygon in self.__polygons:
size += len(polygon)
file.writelines("POLYGONS %d %d \n" % (len(self.__polygons),
size
+ len(self.__polygons)))
for polygon in self.__polygons:
file.writelines("%d " % len(polygon))
for label in polygon:
# I don't know how many labels are there in general.
file.writelines("%d " % label)
file.writelines("\n")
file.close()
def fileName(self):
return self.__fileName
3
你可以考虑使用 vtkAppendPolyData 来合并这些文件。
from vtk import *
reader = vtkPolyDataReader()
append = vtkAppendPolyData()
filenames = ['file1.vtk', 'file2.vtk']
for file in filenames:
reader.SetFileName(file)
reader.Update()
polydata = vtkPolyData()
polydata.ShallowCopy(reader.GetOutput())
append.AddInputData(polydata)
append.Update()
writer = vtkPolyDataWriter()
writer.SetFileName('output.vtk')
writer.SetInput(append.GetOutput())
writer.Write()
另外，要注意的是，如果你第二次调用读取器，它会把第一个输入文件的输出数据替换成第二个输入文件的数据。如果你想使用同一个读取器，就必须像上面的脚本那样做一个浅拷贝。