使用python将多个旧ASCII.vtk文件合并到一个文件中

2024-10-01 09:40:42 发布

您现在位置:Python中文网/ 问答频道 /正文

在计算过程中,我保存了一堆与单个时间点相关的.vtk文件。每一个文件都使用多数据ASCII文件格式描述一个多面体(一个C++类已经定义)。多面体由polyhedron类的一个简单成员函数编写。在

<>为了避免为我需要可视化和扰乱我的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:

^{pr2}$

文件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个坐标和单元格(面),根本没有添加。在

是我还是vtkArray和类似的vtkObjects的python包装器没有被迭代的选项?在


Tags: 文件theinforcellreaderwritervtk
2条回答

您应该考虑使用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()

另外,请注意,第二次调用读取器将用第二个输入文件覆盖第一个输入文件中的输出数据集。如果你想使用同一个阅读器,你必须做一个浅拷贝(如上面的脚本)。在

为了这个目的,我放弃了python-vtk-api,这里有一个小的快速编写的类,有人可能会发现它很有用:

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

相关问题 更多 >