我的问题的背景是我有一个保存在 .vtk 文件中的 3D 结构，我需要对其进行操作（扩张、侵蚀等）。以下代码片段被设计为按顺序运行，即如果您一个接一个地运行它们，应该没有问题（除了我提到的那些！）。
我对VTK很陌生，所以对任何非常基本的错误表示歉意！
问题
我的问题源于 SimpleITK 的问题，它无法读取 UnstructuredGrid 或 PolyData：
In [1]: import SimpleITK as sitk
In [2]: img_vtk = sitk.ReadImage(file_vtk)
Traceback (most recent call last):
File "<ipython-input-52-435ce999db50>", line 1, in <module>
img_vtk = sitk.ReadImage(file_vtk)
File "/usr/local/lib/python3.5/dist-packages/SimpleITK/SimpleITK.py", line 8614, in ReadImage
return _SimpleITK.ReadImage(*args)
RuntimeError: Exception thrown in SimpleITK ReadImage: /tmp/SimpleITK/Code/IO/src/sitkImageReaderBase.cxx:97:
sitk::ERROR: Unable to determine ImageIO reader for "/data/ROMPA_MRIandSeg/09S/Analysis/1_model/clip_dilate.vtk"
然而，SimpleITK 可以读取 StructuredGrid，所以我尝试通过使用 VTK 读取和转换来解决这个问题。
import vtk
reader = vtk.vtkGenericDataObjectReader() # Using generic to allow it to match either Unstructured or PolyData
reader.SetFileName(file_vtk)
reader.Update()
output = reader.GetOutput()
然而，从那时起，我尝试过的每一种方法似乎都失败了。
建议的解决方案
转换为numpy，然后转换为sitk image
我试图将它转换为一个 numpy 数组 ()，然后插入一个规则网格，虚拟变量为 1 以指定结构上的值。
from vtk.utils import numpy_support
import scipy.interpolate
import numpy as np
nparray = numpy_support.vtk_to_numpy(output.GetPointData().GetArray(0))
output_bounds = output.GetBounds()
x_grid = range(math.floor(output_bounds[0]),math.ceil(output_bounds[1]),1)
y_grid = range(math.floor(output_bounds[2]),math.ceil(output_bounds[3]),1)
z_grid = range(math.floor(output_bounds[4]),math.ceil(output_bounds[5]),1)
grid = list()
for x in x_grid:
for y in y_grid:
for z in z_grid:
grid.append((x,y,z))
dummy = np.array([1 for i in range(nparray.shape[0])])
npgrid = scipy.interpolate.griddata(nparray,dummy,grid,fill_value=0)
npgrid.reshape(len(x_grid),len(y_grid),len(z_grid))
img = sitk.GetImageFromArray(npgrid)
sitk.WriteImage(img,file_out)
但是，当我在 ParaView 中加载它时，会为输出显示一个边界框，但输出的轮廓是空的。
使用ShepardMethod
在将 UnstructuredGrid 转换为 PolyData 之后，我尝试使用内置的 进行插值ShepardMethod（因为我经常看到 ShepardMethod 被应用于 PolyData）：
bounds = output.GetBounds()
spacings = [1.0,1.0,1.0] # arbitrary spacing
dimensions = [0,0,0]
for i,spacing in enumerate(spacings):
dimensions[i] = int(math.ceil((bounds[i*2 + 1]-bounds[i*2])/spacing))
vtkPoints = vtk.vtkPoints()
for i in range(0,nparray.shape[0]):
x=nparray[i,0]
y=nparray[i,1]
z=nparray[i,2]
p=[x,y,z]
vtkPoints.InsertNextPoint(p)
poly = vtk.vtkPolyData()
poly.SetPoints(vtkPoints)
shepard = vtk.vtkShepardMethod()
shepard.SetInputData(poly)
shepard.SetSampleDimensions(dimensions)
shepard.SetModelBounds(output.GetBounds())
shepard.Update()
shepard_data = shepard.GetOutput().GetPointData().GetArray(0)
shepard_numpy = numpy_support.vtk_to_numpy(shepard_data)
shepard_numpy = shepard_numpy.reshape(dimensions[0],dimensions[1],dimensions[2])
shepard_img = sitk.GetImageFromArray(shepard_numpy)
sitk.WriteImage(shepard_img,file_out)
与上面的 numpy 工作一样，这在 ParaView 中提供了一个边界框。应用轮廓提供了两个三角形的结构，即几乎没有任何东西似乎已成功写入。或者，我尝试使用 VTK 直接编写输出。
shepard_data = shepard.GetOutput()
shepard_grid = vtk.vtkImageToStructuredGrid()
shepard_grid.SetInputData(shepard_data)
shepard_grid.Update()
writer = vtk.vtkStructuredGridWriter()
writer.SetFileName(file_out)
writer.SetInputData(shepard_grid.GetOutput())
writer.Write()
这产生了与以前相同的输出。
使用ProbeFilter
我尝试使用上述方法ProbeFilter（同时转换为numpy直接写入）。不幸的是，输出与上面相同。
mesh = vtk.vtkStructuredGrid()
mesh.SetDimensions(dimensions)
probe = vtk.vtkProbeFilter()
probe.SetInputData(mesh)
probe.SetSourceData(output)
probe.Update()
probe_out = probe.GetOutput()
writer = vtk.vtkStructuredGridWriter()
writer.SetFileName(file_out)
writer.SetInputData(probe.GetOutput())
writer.Write()
probe_data = probe.GetOutput().GetPointData().GetArray(0)
probe_numpy = numpy_support.vtk_to_numpy(probe_data)
probe_numpy = probe_numpy.reshape(dimensions[0],dimensions[1],dimensions[2])
probe_img = sitk.GetImageFromArray(probe_numpy)
sitk.WriteImage(probe_img,file_out)
但是，这似乎没有产生可行的输出（vtkStructuredGridWriter产生一个空文件，并且probe_numpy是空的）。
更改 ParaView 输出
我的原始数据来自一个结构化网格 .vtk 文件，我使用 ParaView 打开该文件，然后剪辑以删除网格中不需要的结构。保存输出保存了一个非结构化网格，我一直无法弄清楚我是否可以改变它，并首先避免这种混乱！