VTKの3次元データを断面抽出してmatplotlibで画像化する手順【Python、スカラー値】

pip install vtk
pip install numpy
pip install scipy
pip install matplotlib

# === 読み込み部分 ===
reader = vtk.vtkDataSetReader() # インスタンスの作成
reader.SetFileName(filename) # 読み込むファイルのPATHを指定
reader.Update() # 実際にファイルを読み込む
dataset = reader.GetOutput() # 読み込んだ3Dデータ（点、セル、スカラー値など）を含むオブジェクトを取り出す

# === スライス平面を定義 ===
plane = vtk.vtkPlane() # インスタンスの作成
plane.SetOrigin(0, 0, 25.0)  # 平面が通る点(X,Y,Z)
plane.SetNormal(0, 0, 1)  # 平面の法線ベクトル(Nx, Ny, Nz)

# === スライスを作成 ===
cutter = vtk.vtkCutter() # インスタンスの作成
cutter.SetCutFunction(plane) # スライスに使う平面を指定
cutter.SetInputData(dataset) # スライスされるデータを指定
cutter.Update() # スライス処理を実行
cut_output = cutter.GetOutput() # スライスした結果をオブジェクトとして取り出す

points = cut_output.GetPoints().GetData()
scalars = cut_output.GetPointData().GetArray(field_name)

cell_centers = vtk.vtkCellCenters() # インスタンスを作成
cell_centers.SetInputData(cut_output) # セル中心を計算したいデータを指定
cell_centers.Update()  # セルの重心を計算

points = cell_centers.GetOutput().GetPoints().GetData()
scalars = cut_output.GetCellData().GetArray(field_name)

coords = numpy_support.vtk_to_numpy(points)
values = numpy_support.vtk_to_numpy(scalars)

from scipy.interpolate import griddata
grid_z = griddata((x, y), values, (grid_x, grid_y), method="linear")

# スライス面の座標を取得
x, y, z = coords[:, 0], coords[:, 1], coords[:, 2]

# グリッドの解像度を指定
grid_x, grid_y = np.mgrid[min(x) : max(x) : 200j, min(y) : max(y) : 200j]
grid_z = griddata((x, y), values, (grid_x, grid_y), method="linear")

import matplotlib.pyplot as plt

plt.figure(figsize=(6, 5)) # 描画ウィンドウを作成。サイズを幅6インチ、高さ5インチに設定
plt.imshow(
    grid_z.T, # カラーマップの値。imshowは行=Y、列=Xとして扱うので転置
    extent=(min(x), max(x), min(y), max(y)), # 座標のスケーリングを元の空間のx-yに合わせる
    origin="lower", # 座標の原点を左下にする（デフォルトは左上）
    cmap="viridis", # 色のグラデーションを指定
    aspect="equal", # XとYを同じ単位にする（アスペクト比=1）
)
plt.colorbar(label="Scalar Value") # カラーバーのラベルの設定
plt.title("Interpolated Slice as Color Map") # 図のタイトル
plt.xlabel("X") # 横軸のラベル
plt.ylabel("Y") # 縦軸のラベル
plt.tight_layout() # ラベルやタイトルがはみ出ないようにレイアウトを自動調整する
plt.savefig("interpolated_slice.png", dpi=300) # 解像度を300dpiにしてpng画像として保存する
plt.show() # 画面に画像を表示する

import vtk
import numpy as np

# データ定義位置（point or cell）
field_position = "point"
# field_position = "cell"

# 構造格子のサイズ
nx, ny, nz = 50, 50, 50

# スカラー値（例：中心からの距離の2乗）
scalars = np.zeros((nx, ny, nz))
for k in range(nz):
    for j in range(ny):
        for i in range(nx):
            x = i - nx / 2
            y = j - ny / 2
            z = k - nz / 2
            scalars[i, j, k] = x**2 + y**2 + z**2

# vtkImageDataに変換
image_data = vtk.vtkImageData()
if field_position == "point":
    image_data.SetDimensions(nx, ny, nz)
else:
    image_data.SetDimensions(nx + 1, ny + 1, nz + 1)  # セル数より1点多くなることに注意
image_data.SetSpacing(1.0, 1.0, 1.0)

# VTK用スカラー配列に変換
flat_scalars = scalars.ravel(order="F")  # VTKはFortran順
vtk_array = vtk.vtkFloatArray()
vtk_array.SetName("Distance2")
vtk_array.SetNumberOfValues(len(flat_scalars))
for i, val in enumerate(flat_scalars):
    vtk_array.SetValue(i, val)

# スカラーをセット
if field_position == "point":
    image_data.GetPointData().SetScalars(vtk_array)
else:
    image_data.GetCellData().SetScalars(vtk_array)


# ファイルに保存
writer = vtk.vtkStructuredPointsWriter()
writer.SetFileName("sample.vtk")
writer.SetInputData(image_data)
writer.Write()

import vtk
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import griddata
from vtk.util import numpy_support

# VTKファイルのパスを指定
filename = "sample.vtk"

# フィールド名を指定（可視化するデータの名前）
field_name = "Distance2"

# フィールドの定義点を指定（各点にデータがあるのでpointを指定）
field_position = "point"
# field_position = "cell"

# === 読み込み部分 ===
reader = vtk.vtkDataSetReader()  # インスタンスの作成
reader.SetFileName(filename)  # 読み込むファイルのPATHを指定
reader.Update()  # 実際にファイルを読み込む
dataset = (
    reader.GetOutput()
)  # 読み込んだ3Dデータ（点、セル、スカラー値など）を含むオブジェクトを取り出す

# === スライス平面を定義 ===
plane = vtk.vtkPlane()  # インスタンスの作成
plane.SetOrigin(0, 0, 25.0)  # 平面が通る点(X,Y,Z)
plane.SetNormal(0, 0, 1)  # 平面の法線ベクトル(Nx, Ny, Nz)

# === スライスを作成 ===
cutter = vtk.vtkCutter()  # インスタンスの作成
cutter.SetCutFunction(plane)  # スライスに使う平面を指定
cutter.SetInputData(dataset)  # スライスされるデータを指定
cutter.Update()  # スライス処理を実行
cut_output = cutter.GetOutput()  # スライスした結果をオブジェクトとして取り出す

# === スライス結果から座標とスカラー値を取得 ===
if field_position == "cell":
    cell_centers = vtk.vtkCellCenters()  # インスタンスを作成
    cell_centers.SetInputData(cut_output)  # セル中心を計算したいデータを指定
    cell_centers.Update()  # セルの重心を計算
    points = cell_centers.GetOutput().GetPoints().GetData()
    scalars = cut_output.GetCellData().GetArray(field_name)

elif field_position == "point":
    points = cut_output.GetPoints().GetData()
    scalars = cut_output.GetPointData().GetArray(field_name)
else:
    assert 0

if points is None or scalars is None:
    raise ValueError("スライス結果が空です。平面の位置を確認してください。")

# NumPy に変換
coords = numpy_support.vtk_to_numpy(points)
values = numpy_support.vtk_to_numpy(scalars)

# 座標を取得
x, y, z = coords[:, 0], coords[:, 1], coords[:, 2]

# グリッドの解像度を指定して線形補間、構造格子状のデータに整形
grid_x, grid_y = np.mgrid[min(x) : max(x) : 200j, min(y) : max(y) : 200j]
grid_z = griddata((x, y), values, (grid_x, grid_y), method="linear")

# === カラーマップとして描画 ===
plt.figure(figsize=(6, 5))  # 描画ウィンドウを作成。サイズを幅6インチ、高さ5インチに設定
plt.imshow(
    grid_z.T,  # カラーマップの値。imshowは行=Y、列=Xとして扱うので転置
    extent=(
        min(x),
        max(x),
        min(y),
        max(y),
    ),  # 座標のスケーリングを元の空間のx-yに合わせる
    origin="lower",  # 座標の原点を左下にする（デフォルトは左上）
    cmap="viridis",  # 色のグラデーションを指定
    aspect="equal",  # XとYを同じ単位にする（アスペクト比=1）
)
plt.colorbar(label="Scalar Value")  # カラーバーのラベルの設定
plt.title("Interpolated Slice as Color Map")  # 図のタイトル
plt.xlabel("X")  # 横軸のラベル
plt.ylabel("Y")  # 縦軸のラベル
plt.tight_layout()  # ラベルやタイトルがはみ出ないようにレイアウトを自動調整する
plt.savefig("interpolated_slice.png", dpi=300)  # png画像として保存