VTKからPyVistaへの移行

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

たとえば，VTK Pythonのバインディングを使って， vtk.vtkImageData データ構造のポイントをハードコーディングするには，次のように書きます．

>>> import vtk
>>> from math import cos, sin

Create (x, y) points for a 300x300 image dataset

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

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

Create the image structure

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

Assign the points to the image

>>> image_data.GetPointData().SetScalars(points)

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

>>> import pyvista
>>> import numpy as np

Use the meshgrid function to create 2D "grids" of the x and y values.
This section effectively replaces the 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))

Create the grid. Note how the values must use Fortran ordering.

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

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

1. PyVistaでは，データの次元（ numpy.ndarray の形をしている）とデータの値を1行にまとめています．VTKでは，データの形状(空間上の位置)を表すために "タプル" を使用し，データの種類(1=スカラー/スカラーフィールド，2=ベクトル/ベクトルフィールド，n=テンソル/テンソルフィールド)を表すために "コンポーネント" を使用しています．ここでは，1つの変数に形状と値が具体的に格納されています．

2. pyvista.UniformGrid は vtk.vtkImageData をラップしたものです．名前が違うだけで，どちらも等間隔の点のコンテナです． vtk.vtkImageData で使用するデータは "画像" である必要はありません．

   さらに，コンテナが等間隔のデータ用であることがわかっているので， pyvista はデフォルトで原点と間隔を (0, 0, 0) と (1, 1, 1) に設定します．これは，PyVistaとPythonのもう一つの素晴らしい点です．VTKライブラリのすべてを前もって知っておく必要はなく，非常に簡単に始めることができます．慣れてきて，もっと複雑なことをする必要が出てきたら，もっと深く掘り下げることができます．例えば，原点と間隔の変更は以下のように簡単にできます．

   >>> grid.origin = (10, 20, 10)
   >>> grid.spacing = (2, 3, 5)
   

3. 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)
>>> print(vtk_points)

vtkPoints (0x561077401ae0)
  Debug: Off
  Modified Time: 15598
  Reference Count: 1
  Registered Events: (none)
  Data: 0x561077276120
  Data Array Name: Points
  Number Of Points: 3
  Bounds: 
    Xmin,Xmax: (0, 1)
    Ymin,Ymax: (0, 0.667)
    Zmin,Zmax: (0, 0)

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.PolyData のコンストラクタに np_points の配列を渡すだけです．

>>> import pyvista
>>> poly_data = pyvista.PolyData(np_points)

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

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

面やセルの接続性/トポロジーを割り当てる場合も同様です． VTKでは通常， InsertNextCell と InsertCellPoint を使ってループする必要があります． 例えば，一つのセル(3角形)を作成して，それを 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 属性からアクセス（変更）することができます．

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

array([3, 0, 1, 2])

オブジェクトの表現

VTKもPyVistaもオブジェクトの表現を提供しています．

VTK は print() でアクセスできるデータ型の詳細な表現を提供しています(デバッギングに便利)．これは __repr__ が ( __str__ とは異なり) 各オブジェクトに関する最小限の情報しか提供していないためです．

>>> print(vtk_poly_data)

vtkPolyData (0x5610774a18c0)
  Debug: Off
  Modified Time: 15659
  Reference Count: 1
  Registered Events: (none)
  Information: 0x5610774a1210
  Data Released: False
  Global Release Data: Off
  UpdateTime: 0
  Field Data:
    Debug: Off
    Modified Time: 15640
    Reference Count: 1
    Registered Events: (none)
    Number Of Arrays: 0
    Number Of Components: 0
    Number Of Tuples: 0
  Number Of Points: 3
  Number Of Cells: 1
  Cell Data:
    Debug: Off
    Modified Time: 15643
    Reference Count: 1
    Registered Events: 
      Registered Observers:
        vtkObserver (0x5610774a1f80)
          Event: 33
          EventName: ModifiedEvent
          Command: 0x5610774a11c0
          Priority: 0
          Tag: 1
    Number Of Arrays: 0
    Number Of Components: 0
    Number Of Tuples: 0
    Copy Tuple Flags: ( 1 1 1 1 1 0 1 1 1 1 1 1 )
    Interpolate Flags: ( 1 1 1 1 1 0 0 1 1 1 1 0 )
    Pass Through Flags: ( 1 1 1 1 1 1 1 1 1 1 1 1 )
    Scalars: (none)
    Vectors: (none)
    Normals: (none)
    TCoords: (none)
    Tensors: (none)
    GlobalIds: (none)
    PedigreeIds: (none)
    EdgeFlag: (none)
    Tangents: (none)
    RationalWeights: (none)
    HigherOrderDegrees: (none)
    ProcessIds: (none)
  Point Data:
    Debug: Off
    Modified Time: 15642
    Reference Count: 1
    Registered Events: 
      Registered Observers:
        vtkObserver (0x5610774a1190)
          Event: 33
          EventName: ModifiedEvent
          Command: 0x5610774a11c0
          Priority: 0
          Tag: 1
    Number Of Arrays: 0
    Number Of Components: 0
    Number Of Tuples: 0
    Copy Tuple Flags: ( 1 1 1 1 1 0 1 1 1 1 1 1 )
    Interpolate Flags: ( 1 1 1 1 1 0 0 1 1 1 1 0 )
    Pass Through Flags: ( 1 1 1 1 1 1 1 1 1 1 1 1 )
    Scalars: (none)
    Vectors: (none)
    Normals: (none)
    TCoords: (none)
    Tensors: (none)
    GlobalIds: (none)
    PedigreeIds: (none)
    EdgeFlag: (none)
    Tangents: (none)
    RationalWeights: (none)
    HigherOrderDegrees: (none)
    ProcessIds: (none)
  Bounds: 
    Xmin,Xmax: (0, 1)
    Ymin,Ymax: (0, 0.667)
    Zmin,Zmax: (0, 0)
  Compute Time: 15711
  Editable: false
  Number Of Points: 3
  Point Coordinates: 0x561077401ae0
  PointLocator: 0
  CellLocator: 0
  Number Of Vertices: 0
  Number Of Lines: 0
  Number Of Polygons: 1
  Number Of Triangle Strips: 0
  Number Of Pieces: 1
  Piece: -1
  Ghost Level: 0
  CellsBounds: 
    Xmin,Xmax: (0, 1)
    Ymin,Ymax: (0, 0.667)
    Zmin,Zmax: (0, 0)
  CellsBounds Time: 15712

PyVista は repr() で最小限のデータを表示することにしており，アトリビュートの大部分はメッシュ上の明示的なアトリビュートアクセスを優先しています．例えば，以下のようになります．

>>> poly_data

PolyData	Information
N Cells	1
N Points	3
N Strips	0
X Bounds	0.000e+00, 1.000e+00
Y Bounds	0.000e+00, 6.670e-01
Z Bounds	0.000e+00, 0.000e+00
N Arrays	0

この表現では，次のようになります．

• セルの数 n_cells

• ポイントの数 n_points

• メッシュの境界 bounds

• データ配列の数 n_arrays

lines， point_data， cell_data などの他のすべての属性は，オブジェクトから直接アクセスすることができます． このアプローチは，ユーザーを圧倒することなく， DataSet の重要な部分を示す簡単な要約を可能にするために選ばれました．

トレードオフ

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

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

import pyvista

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

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

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

VTKクラスを使用して，その出力をPyVistaでラッピングすることを妨げるものは何もないことに注意してください． 例えば，以下のようになります．

import vtk
import pyvista

# 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 = pyvista.wrap(polygonSource.GetOutput())
mesh.plot(line_width=3, cpos='xy', color='k')

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

注釈

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