TECHNOLOGY

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

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

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

目次

  1. VTKファイルの入力方法と出力方法
    1. サンプルコード
    2. 説明
      1. 読み込み
      2. 書き込み
  2. スカラーデータの編集方法
    1. サンプルコード
    2. 説明
      1. Point上のデータを編集
      2. Cell上のデータを編集
    3. 結果
  3. ベクトルデータの編集方法
    1. サンプルコード１
      1. 説明
      2. 結果
    2. サンプルコード２
      1. 説明
      2. 結果
  4. まとめ

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

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ファイルを編集する上での基礎になるものだと思います。

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

-TECHNOLOGY
-,