PythonによるVTKファイルのスカラー量とベクトル量の読み込みと編集方法

PythonによりVTKファイルを編集するには、VTKファイルを読み込んで、データを加工した後に、VTKファイルで出力することが必要になります。この記事では、それらの方法について記述します。

本記事のコードを実行するには、pythonにvtkパッケージをインストールする必要があります。（pip install vtkなど）

VTKファイルの入力方法と出力方法

サンプルコード

まず、VTKファイルを読み込んで、そのまま出力するためのサンプルコードを示します。

import vtk

# インスタンスを作成
reader = vtk.vtkDataSetReader()
# 読み込むファイルのパスを設定
reader.SetFileName("./sample.vtk")
# データを読み込む
reader.Update()
# 読み込んだデータの種類を確認
data = reader.GetOutput()

# ... ここにVTKファイルからデータを取得したり、加工したりするプログラムを記述します ...

# インスタンスを作成
writer = vtk.vtkDataSetWriter()
# 出力するファイルのパスを設定
writer.SetFileName("./result.vtk")
# 出力するデータを設定
writer.SetInputData(data)
# データをファイルに書き出す
writer.Write()

説明

読み込み

まず、VTKファイルを読み込むために、vtk.vtkDataSetReaderクラスのインスタンス（リーダー）を作成します。次に、SetFileNameメソッドで、読み込みたいファイルのパスを指定します（この例では、カレントディレクトリにあるsample.vtkファイルを指定しています）。その後、Updateメソッドを実行して、リーダーが実際にファイルを読み込みます。最後に、GetOutputメソッドで読み込んだデータを取得し、それをdataという名前の変数に保存します。

書き込み

まず、VTKファイルを出力するために、vtk.vtkDataSetWriterクラスのインスタンス（ライター）を生成します。次に、SetFileNameメソッドで、書き出すファイルのパスを指定します（この例では、カレントディレクトリにあるresult.vtkファイルを指定しています）。その後、SetInputDataで書き出す対象のデータを指定しています。ここで、指定されているdataは、前段階で読み込んだデータセットです。その後、Writeメソッドを実行して、最終的な書き出し操作を実行します。これにより、指定されたファイルにデータを保存します。

スカラーデータの編集方法

VTKファイルを読み込んだ後に、スカラーデータを編集する方法を説明します。スカラー量というのは値を１つしかもたない量のことです。例えば、温度や圧力、体積などはスカラー量です。

ここでは、サンプルファイルとしてOpenFOAMのCavityを計算して得られたVTKファイルを使用します（Cavity流れの詳細についてはこちらを引用します）。

ざっくりと書くと、箱に入った流体があって、箱の上面を動かした際にできる流れのことです。

サンプルコードでは、圧力を意味するスカラーデータ（変数名：p）を読み込んで、その値に101325 Pa（大気圧）を加算する方法を記載します。

使用するVTKファイル（ZIPに圧縮しております）：

cavity.vtk

1 ファイル 15 KB

ダウンロード

サンプルコード

サンプルコードを示します。「VTKファイルの入力方法と出力方法」の入出力部分の間に、スカラーデータ（変数名：p）を編集するコードを追加しています。

変数 p を編集した後に、新しい変数 p_abs としてVTKファイルにデータを追加しています。

VTKファイルでは、Point上で定義した変数とCell上で定義した変数があります。そのそれぞれについて編集を行っています。

import vtk

# インスタンスを作成
reader = vtk.vtkDataSetReader()
# 読み込むファイルのパスを設定
reader.SetFileName("./cavity.vtk")
# データを読み込む
reader.Update()
# 読み込んだデータの種類を確認
data = reader.GetOutput()

##### Point上のデータを編集 (start)
# 編集したデータを格納する変数を定義する
pnew = vtk.vtkDoubleArray()
pnew.SetName('p_abs')
pnew.SetNumberOfComponents(1)

# データＰの配列を取得
array = data.GetPointData().GetArray('p') 

for i in range(array.GetNumberOfTuples()):
    # 各Pointでの値を取得
    value = array.GetValue(i)
    new_value = value + 101325.0
    # 新しく定義した変数に代入
    pnew.InsertNextValue(new_value)

# データに新しく定義した変数を追加
data.GetPointData().AddArray(pnew)
##### Point上のデータを編集 (end)

##### Cell上のデータを編集 (start)
pnew = vtk.vtkDoubleArray()
pnew.SetName('p_abs')
pnew.SetNumberOfComponents(1)

# データＰの配列を取得
array = data.GetCellData().GetArray('p') 

for i in range(array.GetNumberOfTuples()):
    # 各Cellでの値を取得
    value = array.GetValue(i)
    new_value = value + 101325.0
    # 新しく定義した変数に代入
    pnew.InsertNextValue(new_value)

# データに新しく定義した変数を追加
data.GetCellData().AddArray(pnew)
##### Cell上のデータを編集 (end)

# インスタンスを作成
writer = vtk.vtkDataSetWriter()
# 出力するファイルのパスを設定
writer.SetFileName("./result.vtk")
# 出力するデータを設定
writer.SetInputData(data)
# データをファイルに書き出す
writer.Write()

説明

このコードは、VTKのデータセット内にある Point と Cell のデータを編集し、新しいデータ配列として追加する処理を示しています。具体的には、PointData と CellData に保存されている圧力データ（p）に対して操作を行い、絶対圧（p_abs）として新しいデータ配列を作成しています。

Point上のデータを編集

まず、新しいデータ配列を定義します。

pnew = vtk.vtkDoubleArray()
pnew.SetName('p_abs')
pnew.SetNumberOfComponents(1)

vtkDoubleArray クラスを使って、新しいデータ配列 pnew を作成しています。この配列には、計算された絶対圧 p_abs の値が格納されます。次に、SetNameメソッドを用いて、この新しい配列の名前を 'p_abs' としています。そして、SetNumberOfComponentsでは、各ポイントにつき1つの値（コンポーネント）を持つ配列を定義しています。つまり、スカラー量です。

新しいデータ配列を定義したら、元のデータ配列を取得します。

array = data.GetPointData().GetArray('p')

data オブジェクト内の PointData に保存されている p という名前のデータ配列を取得します。このデータは各ポイントの圧力値を表します。

元のデータを取得したら、各ポイントの圧力値を取得し、それに101325.0 Paを加算した新しい値を計算します。それを新しいデータ配列pnew に追加します。

for i in range(array.GetNumberOfTuples()):
    value = array.GetValue(i)
    new_value = value + 101325.0
    pnew.InsertNextValue(new_value)

最後に、計算された新しい配列pnewを、データセットの'PointDataに追加します。

data.GetPointData().AddArray(pnew)

Cell上のデータを編集

Cell上のデータを編集する場合もPoint上のデータを編集する場合とほとんど同じです。異なるのは、元のデータ配列を取得する際にはGetCellDataを使用します。

array = data.GetCellData().GetArray('p')

また、計算された新しい配列pnewを、データセットのCellDataに追加する場合にも、

data.GetCellData().AddArray(pnew)

を用います。

結果

左がもとのデータ、右が新しく追加したデータになります。値に10325 Paが加算されていることが分かります。

ベクトルデータの編集方法

VTKファイルを読み込んだ後に、ベクトルデータを編集する方法を説明します。ベクトル量というのは値と方向をもつ量のことです。例えば、速度はどの方向に何メートルという形で定義されるのでベクトル量です。

ここでは、サンプルファイルとして先ほどと同じくOpenFOAMのCavityを計算して得られたVTKファイルを使用します。

サンプルコード１

サンプルコード１では、速度を意味するベクトルデータ（変数名：U）を読み込んで、そのベクトルの大きさを計算して新しい変数としてVTKファイルに追加する方法を記載します。

「VTKファイルの入力方法と出力方法」の入出力部分の間に、ベクトルデータ（変数名：U）を編集するコードを追加しています。変数 U を編集した後に、新しい変数としてVTKファイルにデータを追加しています。

VTKファイルでは、Point上で定義した変数とCell上で定義した変数があります。ここでは、PointDataについてサンプルプログラムを示します。

スカラー量で説明したように、CellDataについては、GetPointDataメソッドをGetCellDataメソッドに書き換えることで同じことができます。

ベクトルデータを取得する方法として、GetComponentメソッドを使用する方法とGetTupleメソッドを使用する方法について記述します。

import vtk
import math

# インスタンスを作成
reader = vtk.vtkDataSetReader()
# 読み込むファイルのパスを設定
reader.SetFileName("./cavity.vtk")
# データを読み込む
reader.Update()
# 読み込んだデータの種類を確認
data = reader.GetOutput()

# ベクトルの大きさを計算して追加する(例1)
magU = vtk.vtkDoubleArray()
magU.SetName("magU")
magU.SetNumberOfComponents(1) # スカラーなので1

U = data.GetPointData().GetArray("U")
#U = data.GetCellData().GetArray("U")

for i in range(U.GetNumberOfTuples()):
    Ux = U.GetComponent(i, 0)
    Uy = U.GetComponent(i, 1)
    Uz = U.GetComponent(i, 2)
    magU.InsertNextValue(math.sqrt(Ux*Ux + Uy*Uy + Uz*Uz))

data.GetPointData().AddArray(magU)
#data.GetCellData().AddArray(magU)

# ベクトルの大きさを計算して追加する(例2)
magV = vtk.vtkDoubleArray()
magV.SetName("magV")
magV.SetNumberOfComponents(1) # スカラーなので1

V = data.GetPointData().GetArray("U")
#V = data.GetCellData().GetArray("U")
for i in range(V.GetNumberOfTuples()):
    v = list(V.GetTuple(i))
    magV.InsertNextValue(math.sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]))

data.GetPointData().AddArray(magV)
#data.GetCellData().AddArray(magV)

# インスタンスを作成
writer = vtk.vtkDataSetWriter()
# 出力するファイルのパスを設定
writer.SetFileName("./result.vtk")
# 出力するデータを設定
writer.SetInputData(data)
# データをファイルに書き出す
writer.Write()

説明

スカラー量と異なる点は、ベクトル量では成分が複数ある点です。このため、ベクトルデータの成分を取り出すところについて説明します。

方法の１つめは、GetComponentメソッドを使用する方法です。

for i in range(U.GetNumberOfTuples()):
    Ux = U.GetComponent(i, 0)
    Uy = U.GetComponent(i, 1)
    Uz = U.GetComponent(i, 2)
    magU.InsertNextValue(math.sqrt(Ux*Ux + Uy*Uy + Uz*Uz))

GetComponentの第１引数は、i番目のデータ点上のデータを意味しており、第２引数はベクトルの成分（0=x成分、1=y成分、2=z成分）を意味しています。

これにより、各成分をUx, Uy, Uzに格納した後に、math.sqrtメソッドでベクトルの大きさを求めています。

方法の２つめは、GetTupleメソッドを使用する方法です。

for i in range(V.GetNumberOfTuples()):
    v = list(V.GetTuple(i))
    magV.InsertNextValue(math.sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]))

これにより、i番目のデータ点上のデータをリストとして取得しています。

結果

新しい変数が追加され、次の結果が得られます。

サンプルコード２

サンプルコード２では、速度を意味するベクトルデータ（変数名：U）を読み込んで、方向を反転させてベクトルデータとして保存する例を示します。

物理的な意味はありませんが、ベクトルデータを編集してベクトルデータとして出力する例になります。

VTKファイルでは、Point上で定義した変数とCell上で定義した変数があります。ここでは、PointDataについてサンプルプログラムを示します。

スカラー量で説明したように、CellDataについては、GetPointDataメソッドをGetCellDataメソッドに書き換えることで同じことができます。

import vtk
import math

# インスタンスを作成
reader = vtk.vtkDataSetReader()
# 読み込むファイルのパスを設定
reader.SetFileName("./cavity.vtk")
# データを読み込む
reader.Update()
# 読み込んだデータの種類を確認
data = reader.GetOutput()

invU = vtk.vtkDoubleArray()
invU.SetName("inverse_U")
invU.SetNumberOfComponents(3) # ベクトルなので3
U = data.GetPointData().GetArray("U")
#U = data.GetCellData().GetArray("U")
for i in range(U.GetNumberOfTuples()):
    u = list(U.GetTuple(i))
    for j in range(0, len(u)):
        u[j] = -u[j]
    invU.InsertNextTuple(u)

data.GetPointData().AddArray(invU)
#data.GetCellData().AddArray(invU)

# インスタンスを作成
writer = vtk.vtkDataSetWriter()
# 出力するファイルのパスを設定
writer.SetFileName("./result.vtk")
# 出力するデータを設定
writer.SetInputData(data)
# データをファイルに書き出す
writer.Write()

説明

これまでと同様に編集したベクトルデータを格納するための新しい変数を定義しています。

invU.SetName("inverse_U")
invU.SetNumberOfComponents(3) # ベクトルなので3
U = data.GetPointData().GetArray("U")

異なるのは、ベクトルデータなので成分が3つあります。このため、SetNumberOfComponentsメソッドの引数は3になっています。

ベクトルを反転して、新しい配列に追加する部分は、次のコードで実行しています。

for i in range(U.GetNumberOfTuples()):
    u = list(U.GetTuple(i))
    for j in range(0, len(u)):
        u[j] = -u[j]
    invU.InsertNextTuple(u)

ベクトルデータをもつ変数に値を設定するために、InsertNextTupleメソッドを使用している点がスカラー量とは異なります。

上記以外については、これまでのサンプルコードと同様です。

結果

新しい変数が追加されて、次のようにベクトルが反転します。左側がオリジナルのデータで、右側が編集したデータになります。

まとめ

PythonによりVTKファイルを読み込み、それを編集してVTKファイルとして出力する方法を説明しました。

サンプルファイルでは、スカラー量とベクトル量に分けて、データの取得方法を示しました。

これらの方法はVTKファイルを編集する上での基礎になるものだと思います。

参考になりましたら幸いです。