PyVistaとVTK

PyVistaとVTK#

PyVistaがVTKをどのように使っているか，そして両者の長所をどのように組み合わせられるかを紹介します! (トーク10分，エクササイズ10分)

Tip

このチュートリアルのセクションは，PyVista ドキュメントの Transitioning from VTK to PyVista の章から採用されています．

VTKからPyVistaへの移行#

VTKは主にC++で開発されており，データへのアクセスには連鎖したセッターとゲッターのコマンドを使用します．その代わりに，PyVista は VTK のデータタイプを numpy 配列にラップして，ブラケット構文や派手なインデックスの恩恵を受けられるようにしています．このセクションでは，一連の例で2つのアプローチの違いを説明します．

例えば，VTK Python のバインディングを使用して vtk.vtkImageData データ構造の値をハードコードするには，次のように記述します．

import vtk
from math import cos, sin

300x300の画像データセットの値を作成する

values = vtk.vtkDoubleArray()
values.SetName("values")
values.SetNumberOfComponents(1)
values.SetNumberOfTuples(300*300)

for x in range(300):
   for y in range(300):
      values.SetValue(x*300 + y, 127.5 + (1.0 + sin(x/25.0)*cos(y/25.0)))

画像構造の作成

image_data = vtk.vtkImageData()
image_data.SetOrigin(0, 0, 0)
image_data.SetSpacing(1, 1, 1)
image_data.SetDimensions(300, 300, 1)

画像に値を割り当てる

image_data.GetPointData().SetScalars(values)

ご覧のように，単純な vtk.vtkImageData データセットを作成するためには，かなり多くのボイラープレートが必要です．PyVistaでは，より簡潔で，より "Pythonic "な構文が提供されています．PyVistaで同等のコードは以下の通りです．

import pyvista as pv
import numpy as np

meshgrid 関数を使用して，x および y 値の 2D "グリッド" を作成します．ここでは，実質的にvtkDoubleArrayを置き換えています．

xi = np.arange(300)
x, y = np.meshgrid(xi, xi)
values = 127.5 + (1.0 + np.sin(x/25.0)*np.cos(y/25.0))

グリッドを作成します．値がどのようにFortran順序を使用しなければならないかに注意してください．

grid = pv.ImageData(dimensions=(300, 300, 1))
grid.point_data["values"] = values.flatten(order="F")

ここでは，PyVistaがいくつかのことをしてくれています．

PyVistaでは，データの次元（ numpy.ndarray の形をしている）とデータの値を1行にまとめています．VTKでは，データの形状(空間上の位置)を表すために "タプル" を使用し，データの種類(1=スカラー/スカラーフィールド，2=ベクトル/ベクトルフィールド，n=テンソル/テンソルフィールド)を表すために "コンポーネント" を使用しています．ここでは，1つの変数に形状と値が具体的に格納されています．
pyvista.ImageData は vtk.vtkImageData をラップしたものです．名前が違うだけで，どちらも等間隔の点のコンテナです．どちらも等間隔の点のコンテナです． vtk.vtkImageData で使用するデータは "画像" である必要はありません．

さらに，コンテナが等間隔のデータ用であることがわかっているので， pyvista はデフォルトで原点と間隔を (0, 0, 0) と (1, 1, 1) に設定します．これは，PyVistaとPythonのもう一つの素晴らしい点です．VTKライブラリのすべてを前もって知っておく必要はなく，非常に簡単に始めることができます．慣れてきて，もっと複雑なことをする必要が出てきたら，もっと深く掘り下げることができます．例えば，原点と間隔の変更は以下のように簡単にできます．
```
grid.origin = (10, 20, 10)
grid.spacing = (2, 3, 5)
```
point_array <pyvista.point_array> の名前は，辞書形式で直接与えられます．また，VTKはデータをヒープ(RAMの線形セグメント．C++の概念)に保存するため，データをフラット化してFortranの順序(多次元データが物理的に1次元のメモリにどのようにレイアウトされるかを制御する．numpyはデフォルトで "C "スタイルのメモリレイアウトを使用する)にする必要があります．先ほどの例で， SetValue() の最初の引数が x*300 + y と書かれていたのはこのためです．ここでは，numpyがこれをうまく処理してくれるので，Pythonのベストプラクティスである "Explicit is better than implicit" に従って，コードの中でより明示的にしています．

最後に，PyVistaでは，各ジオメトリクラスにメソッドが用意されているので，プロットの設定をしなくても，すぐにメッシュをプロットすることができます．例えば，VTKでは次のようにする必要があります．

actor = vtk.vtkImageActor()
actor.GetMapper().SetInputData(image_data)
ren = vtk.vtkRenderer()
renWin = vtk.vtkRenderWindow()
renWin.AddRenderer(ren)
renWin.SetWindowName('ReadSTL')
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)
ren.AddActor(actor)
iren.Initialize()
renWin.Render()
iren.Start()

しかし，PyVistaでは，必要なのは以下だけです:

grid.plot(cpos='xy', show_scalar_bar=False, cmap='coolwarm')

スタティックなシーン

インタラクティブなシーン

ポイントセット構築#

PyVista は，VTK の C 配列を効率的に割り当て，アクセスするために NumPy に大きく依存しています．例えば，VTKで点の配列を作るには，通常，リストのすべての点をループして，それを vtkPoints クラスに供給します．例えば，以下のようになります．

import vtk
vtk_array = vtk.vtkDoubleArray()
vtk_array.SetNumberOfComponents(3)
vtk_array.SetNumberOfValues(9)
vtk_array.SetValue(0, 0)
vtk_array.SetValue(1, 0)
vtk_array.SetValue(2, 0)
vtk_array.SetValue(3, 1)
vtk_array.SetValue(4, 0)
vtk_array.SetValue(5, 0)
vtk_array.SetValue(6, 0.5)
vtk_array.SetValue(7, 0.667)
vtk_array.SetValue(8, 0)
vtk_points = vtk.vtkPoints()
vtk_points.SetData(vtk_array)
vtk_points

PyVistaで同じことをするには，単にNumPyの配列を作成する必要があります．

import numpy as np
np_points = np.array([[0, 0, 0], [1, 0, 0], [0.5, 0.667, 0]])

注釈

pyvista.vtk_points() を使って vtkPoints オブジェクトを構築することもできますが，ほとんどの状況では必要ありません．

最終的な目的は， pyvista.DataSet <pyvista.core.dataset.DataSet> を構築することなので， pyvista.PolyData のコンストラクタに np_points の配列を渡すだけです．

import pyvista as pv
poly_data = pv.PolyData(np_points)

VTKではそうする必要があります．

vtk_poly_data = vtk.vtkPolyData()
vtk_poly_data.SetPoints(vtk_points)

面やセルの接続性/トポロジーを割り当てる場合も同様です． VTKでは通常， InsertNextCell() と InsertCellPoint() を使ってループする必要があります．例えば，一つのセル(三角形)を作成して，それを vtkPolyData に割り当てる場合:

cell_arr = vtk.vtkCellArray()
cell_arr.InsertNextCell(3)
cell_arr.InsertCellPoint(0)
cell_arr.InsertCellPoint(1)
cell_arr.InsertCellPoint(2)
vtk_poly_data.SetPolys(cell_arr)

PyVistaでは，コンストラクタでこれを直接割り当て， faces <pyvista.PolyData.faces> 属性からアクセス（変更）することができます．

faces = np.array([3, 0, 1, 2])
poly_data = pv.PolyData(np_points, faces)
poly_data.faces

PyVistaのトレードオフ#

ほとんどの機能は可能ですが，PyVistaでは機能や性能を落とさずにすべてを簡略化できるわけではありません．

collision <pyvista.PolyDataFilters.collision> フィルターでは，2つのメッシュ間のコリジョンを計算する方法を示します．例えば，以下のようになります．

import pyvista as pv

# create a default sphere and a shifted sphere
mesh_a = pv.Sphere()
mesh_b = pv.Sphere(center=(-0.4, 0, 0))
out, n_coll = mesh_a.collision(mesh_b, generate_scalars=True, contact_mode=2)

pl = pv.Plotter()
pl.add_mesh(out)
pl.add_mesh(mesh_b, style='wireframe', color='k')
pl.camera_position = 'xy'
pl.show()

スタティックなシーン

インタラクティブなシーン

フードの下では，コリジョンフィルタは OBB (oriented bounding box) ツリーを使ってメッシュの衝突を検出します．しかし， Collision Example の例のように，同じメッシュで複数の衝突を計算する場合には，OBBツリーが各メッシュに対して一度ずつ計算されるため， vtkCollisionDetectionFilter を使用した方が効率的です．ほとんどのデータサイエンスでは，純粋なPyVistaで十分ですが，VTKクラスを直接使用したい場合もあります．

PyVistaとVTKの連携#

PyVistaは，使いやすいプロットAPIでよく知られており，Matplotlibのようなライブラリですでに経験を積んだほとんどのPythonユーザーにとってなじみ深いものです．多くの人が，データパイプラインにVTKのPythonバインディングの力を，3DレンダリングにPyVistaの柔軟性とシンプルさを組み合わせることで利益を得ています．次のセクションでは，この使用シナリオを説明します．

Tip

メッシュとは? のセクションで明確にされていないかもしれませんが，PyVistaのメッシュクラスはVTKの対応するクラスのサブクラスであり，これはPyVistaとVTKのワークフローを混在させることができることを意味しています．

VTKのクラスを使って，出力をPyVistaでラップして，効率的にプロットすることを止めるものはありません．例えば

import vtk
import pyvista as pv

# Create a circle using vtk
polygonSource = vtk.vtkRegularPolygonSource()
polygonSource.GeneratePolygonOff()
polygonSource.SetNumberOfSides(50)
polygonSource.SetRadius(5.0)
polygonSource.SetCenter(0.0, 0.0, 0.0)
polygonSource.Update()

# wrap and plot using pyvista
mesh = pv.wrap(polygonSource.GetOutput())
mesh.plot(line_width=3, cpos='xy', color='k')

スタティックなシーン

インタラクティブなシーン

このようにして， PyVistaの柔軟性とVTKの原動力の両方を必要とする "best of both worlds" を得ることができるのです．

注釈

上記のVTKコードを1行で置き換えるために， pyvista.Polygon() を使用することができます．

VTKアルゴリズム#

おそらく，PyVistaで（まだ）公開されていないVTKアルゴリズムで，使いたいものがあるのでしょう．PyVistaのオブジェクトはVTKのオブジェクトなので，これは十分に簡単に扱えます．PyVistaのメッシュをVTKアルゴリズムに渡して，出力をラップしてプロットしたり，さらにフィルタリングしたり，何でもできるのです．

import pyvista as pv
from pyvista import examples
import vtk

mesh = examples.download_bunny_coarse()

# Initialize VTK algorithm
splatter = vtk.vtkGaussianSplatter()

# Pass PyVista object as input to VTK
splatter.SetInputData(mesh)

# Set parameters
n = 200
splatter.SetSampleDimensions(n, n, n)
splatter.SetRadius(.02)
splatter.SetExponentFactor(-10)
splatter.SetEccentricity(2)
splatter.Update()

# Retrieve output and wrap with PyVista
vol = pv.wrap(splatter.GetOutput())

# Use PyVista to produce contours
cntrs = vol.contour([.95 * splatter.GetRadius()])

# Use PyVista to plot
p = pv.Plotter()
p.add_mesh(mesh, style='wireframe')
p.add_mesh(cntrs, color=True)
p.show()