私の問題の背景は、3D 構造が .vtk ファイルに保存されており、それを操作 (拡張、浸食など) する必要があることです。次のコード スニペットは、順番に実行するように設計されています。つまり、1 つずつ実行しても問題はありません (私が言及したものを除いて!)。
私は VTK に非常に慣れていないので、非常に基本的な間違いをお詫びします!
問題
私の問題は、UnstructuredGrid または PolyData を読み取ることができない SimpleITK の問題から生じています。
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 画像に変換します
それを 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(PolyData に適用される ShepardMethod をほとんど見たように)。
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)
上記の派手な作業と同様に、これにより ParaView にバウンディング ボックスが提供されました。輪郭を適用すると、2 つの三角形の構造が提供されます。つまり、ほとんど何も書き込まれていないようです。または、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 出力の変更
元のデータは、ParaView を使用して開いた構造化グリッド .vtk ファイルから取得し、クリップしてメッシュに不要な構造を削除します。出力を保存すると unstructuredGrid が保存されますが、それを変更できるかどうかを判断できず、そもそもこの混乱を回避できません!