6. NumPy を使用してデータを処理する
5 章 では、NumPyを使って配列にアクセスし、配列に対する操作を実行することに大きく依存していたデータ処理用のPythonスクリプトを作成するためのいくつかのレシピについて説明しました。この章では、2つのシステム間でデータ表現に大きな違いがあるにもかかわらず、VTKとNumPyを一緒に使用できるようにするVTK-NumPy統合レイヤーについて詳しく説明します。
6.2. dataset_adapter モジュールを理解する
本節では、 dataset_adapter モジュールについて詳しく説明します。このモジュールは、PythonからVTKデータセットおよびアレイに簡単にアクセスできるように設計されており、NumPyスタイルのインターフェイスを提供します。
前節の例に進みましょう。このスクリプトは、 Sphere に接続された Programmable Filter の Script に置かれ、その後に Elevation フィルタパイプラインが続きます。
from vtk.numpy_interface import dataset_adapter as dsa
...
print(data)
print(isinstance(data, dsa.VTKObjectWrapper))
これは下記のように出力されます。:
<paraview.vtk.numpy_interface.dataset_adapter.PolyData object at 0x14b7caa50>
True
VTKObjectメンバを使用して、基礎となるVTKオブジェクトにアクセスできます。
print(type(data.VTKObject))
は下記のように出力されます。:
<type 'vtkCommonDataModelPython.vtkPolyData'>
Programmable Filter で$inputs$として取得されるものは、実際にはVTKデータ・オブジェクト自体をラップするPythonオブジェクトです。 Programmable Filter は、VTKデータオブジェクトの vtk.numpy_interface.dataset_adapter モジュールから WrapDataObject 関数を手動で呼び出すことでこれを行います。WrapDataObject 関数は、すべての vtkDataSet サブクラス、vtkTable サブクラス、およびすべての vtkCompositeData サブクラスの適切なラッパークラスを返します。 他の vtkDataObject サブクラスは現在サポートされていません。
VTKObjectWrapper は、VTKメソッドをVTKObjectに転送するので、以下のようにVTK APIに直接アクセスすることができます。
print(data.GetNumberOfCells())
96L
ただし、VTKObjectWrapper を引数としてVTKメソッドに直接渡すことはできません。
from paraview.vtk.vtkFiltersGeneral import vtkShrinkPolyData
s = vtkShrinkPolyData()
s.SetInputData(data)
データを設定しようとすると、エラーメッセージが表示されます。
TypeError: SetInputData argument 1: method requires a VTK object
代わりに、次のようにVTKオブジェクトをVTKフィルタに渡す必要があります。
s.SetInputData(data.VTKObject)
上記の例で注意すべき重要な点は、Pythonクラス vtkShrinkPolyData がスクリプトで使用するためにどのように読み込まれたかです。VTKでは、クラスは関連する機能の異なるグループに編成され、これらのグループはPythonモジュールとしてインポートすることができます。クラスを使用するには、まずそのクラスが存在するモジュールを特定します。これは、[KitwareInc] クラスのDoxygenドキュメントから確認できます。クラスのDoxygenページに移動し、ページの一番下にある dox/<first directory>/<second directory>/<class name> という形式のドキュメントが生成されたファイルのパスを探します。モジュールは vtk <最初のディレクトリ> <セカンドディレクトリ> として導出されます。例として、 vtkShrinkPolyData のドキュメントは dox/Filters/General/vtkShrinkPolyData から生成されているので、モジュールは vtkFiltersGeneral です。次に、フォームのステートメントを使用してクラスをインポートできます。
6.2.1. Dataset属性
これまでのところ、VTKデータオブジェクト用のラッパーがあります。 はVTKデータオブジェクトのように部分的に動作します。このデータセットに含まれるフィールド (配列) にアクセスする方法を調べ始めると、これが少し興味深いものになります。
簡単にするために、スクリプトによって生成された出力をコード自体に埋め込み、>>> プレフィックスを使用してコードと出力を区別します。
>>> print(data.PointData)
<vtk.numpy_interface.dataset_adapter.DataSetAttributes at 0x110f5b750>
>>> print(data.PointData.keys())
['Normals', 'Elevation']
>>> print(data.CellData.keys())
[]
>>> print(data.PointData['Elevation'])
VTKArray([ 0.5 , 0. , 0.45048442, 0.3117449 , 0.11126047,
0. , 0. , 0. , 0.45048442, 0.3117449 ,
0.11126047, 0. , 0. , 0. , 0.45048442,
...,
0.11126047, 0. , 0. , 0. , 0.45048442,
0.3117449 , 0.11126047, 0. , 0. , 0. ], dtype=float32)
>>> elevation = data.PointData['Elevation']
>>> print(elevation[:5])
VTKArray([0.5, 0., 0.45048442, 0.3117449, 0.11126047], dtype=float32)
# Note that this works with composite datasets as well:
>>> mb = vtk.vtkMultiBlockDataSet()
>>> mb.SetNumberOfBlocks(2)
>>> mb.SetBlock(0, data.VTKObject)
>>> mb.SetBlock(1, data.VTKObject)
>>> mbw = dsa.WrapDataObject(mb)
>>> print(mbw.PointData)
<vtk.numpy_interface.dataset_adapter.CompositeDataSetAttributes instance at 0x11109f758>
>>> print(mbw.PointData.keys())
['Normals', 'Elevation']
>>> print(mbw.PointData['Elevation'])
<vtk.numpy_interface.dataset_adapter.VTKCompositeDataArray at 0x1110a32d0>
この方法で、PointData 、CellData、FieldData、Points (vtkPointSet のみのサブクラス)、Polygons (vtkPolyData のみ)にアクセスすることができます。このAPIを使用して、さらに多くのタイプの配列にアクセサを追加します。
6.3. 配列の操作
本節では、テストパイプラインを、 Programmable Filter に接続された Wavelet ソースで構成するように変更します。
Script では、次のように RTData ポイントデータ配列にアクセスします。
# Code for 'Script'
from paraview.vtk.vtkFiltersGeneral import vtkDataSetTriangleFilter
image = inputs[0]
rtdata = image.PointData['RTData']
# Let's transform this data as well, using another VTK filter.
tets = vtkDataSetTriangleFilter()
tets.SetInputDataObject(image.VTKObject)
tets.Update()
# Here, now we need to explicitly wrap the output dataset to get a
# VTKObjectWrapper instance.
ugrid = dsa.WrapDataObject(tets.GetOutput())
rtdata2 = ugrid.PointData['RTData']
ここでは、画像データ( vtkImageData )と非構造格子( vtkUnstructuredGrid )の2つのデータセットを作成しました。基本的には同じデータを表しますが、非構造格子は画像データを四面体化することによって作成されます。したがって、非構造格子は、同じ点を持ちますが、より多くのセル(四面体)を持つことが予想されます。
from paraview.vtk.vtkFiltersGeneral import vtkDataSetTriangleFilter
6.3.1. 配列API
numpy_interface 配列オブジェクトの動作はNumPy配列と非常に似ています。実際、vtkDataSet サブクラスの配列はVTKArrayのインスタンスです。numpy.ndarray のサブクラスです。vtkCompositeDataSet とサブクラスの配列はNumPy配列ではありませんが、非常によく似た動作をします。相違点については、別の項で概説します。基本から始めましょう。次の作業はすべて予想どおりに行われました。
これまでと同様、単純化するために、スクリプトによって生成された出力をコード自体に埋め込み、>>> プレフィックスを使用してコードと出力を区別します。
>>> print(rtdata[0])
60.763466
>>> print(rtdata[-1])
57.113735
>>> print(repr(rtdata[0:10:3]))
VTKArray([ 60.76346588, 95.53707886, 94.97672272, 108.49817657], dtype=float32)
>>> print(repr(rtdata + 1))
VTKArray([ 61.76346588, 86.87795258, 73.80931091, ..., 68.51051331,
44.34006882, 58.1137352 ], dtype=float32)
>>> print(repr(rtdata < 70))
VTKArray([ True , False, False, ..., True, True, True])
# We will cover algorithms later. This is to generate a vector field.
>>> avector = algs.gradient(rtdata)
# To demonstrate that avector is really a vector
>>> print(algs.shape(rtdata))
(9261,)
>>> print(algs.shape(avector))
(9261, 3)
>>> print(repr(avector[:, 0]))
VTKArray([ 25.69367027, 6.59600449, 5.38400745, ..., -6.58120966,
-5.77147198, 13.19447994])
この例では、次の点に注意してください。
単一コンポーネントの配列は、常に次のような形をしています: (n-tuples,) であり、(n-tuples, 1)ではありません。
複数のコンポーネント配列の形状は次のとおりです:(n-tuples, n-components)
テンソル配列の形状:(n-tuples, 3, 3)
上記は、画像や他の構造データにも当てはまります。すべての配列は1つの次元(1成分配列)、2つの次元(多成分配列)または3つの次元(テンソル配列)を持ちます。
もう1つの素晴らしい点は、配列のインデックスを作成するためにboolean配列を使用できることです。したがって、次のようにすると非常に効果的です。
>>> print(repr(rtdata[rtdata < 70]))
VTKArray([ 60.76346588, 66.75043488, 69.19681549, 50.62128448,
64.8801651 , 57.72655106, 49.75050354, 65.05570221,
57.38450241, 69.51113129, 64.24596405, 67.54656982,
...,
61.18143463, 66.61872864, 55.39360428, 67.51051331,
43.34006882, 57.1137352 ], dtype=float32)
>>> print(repr(avector[avector[:,0] > 10]))
VTKArray([[ 25.69367027, 9.01253319, 7.51076698],
[ 13.1944809 , 9.01253128, 7.51076508],
[ 25.98717642, -4.49800825, 7.80427408],
...,
[ 12.9009738 , -16.86548471, -7.80427504],
[ 25.69366837, -3.48665428, -7.51076889],
[ 13.19447994, -3.48665524, -7.51076794]])
6.3.2. アルゴリズム
配列APIを使えば簡単にいろいろなことができますが、numpy_interface.algorithms モジュールを使い始めるともっと面白くなります。前の例で簡単に紹介しました。ここではさらに詳しく説明します。アルゴリズムの完全なリストについては、help(algs) を使用してください。以下に簡単に説明できる例をいくつか示します。
>>> import paraview.vtk.numpy_interface.algorithms as algs
>>> print(repr(algs.sin(rtdata)))
VTKArray([-0.87873501, -0.86987603, -0.52497 , ..., -0.99943125,
-0.59898132, 0.53547275], dtype=float32)
>>> print(repr(algs.min(rtdata)))
VTKArray(37.35310363769531)
>>> print(repr(algs.max(avector)))
VTKArray(34.781060218811035)
>>> print(repr(algs.max(avector, axis=0)))
VTKArray([ 34.78106022, 29.01940918, 18.34743023])
>>> print(repr(algs.max(avector, axis=1)))
VTKArray([ 25.69367027, 9.30603981, 9.88350773, ..., -4.35762835,
-3.78016186, 13.19447994])
軸引数を使用したことがない場合は、非常に簡単です。軸の値を渡さない場合、関数は次元性を考慮せずに配列のすべての値に適用されます。axis=0 の場合、関数は配列の各要素に独立に適用されます。axis=1 の場合、関数は各タプルに独立に適用されます。これがはっきりしない場合は、実験してみてください。 このように動作する関数には、 sum 、 min 、 max 、 std 、 var があります。
もう1つの興味深い便利な関数は、特定の条件が発生した場合に配列のインデックスが返される関数です。
>>> print(repr(algs.where(rtdata < 40)))
(array([ 420, 9240]),)
# For vectors, this will also return the component index if an axis is not
# defined.
>>> print(repr(algs.where(avector < -29.7)))
(VTKArray([4357, 4797, 4798, 4799, 5239]), VTKArray([1, 1, 1, 1, 1]))
これまで説明したすべての関数は、NumPyによって直接提供されています。 NumPyの関数の多くは、algorithmsモジュールに含まれています。これらはすべて単一配列と複合データ配列で動作します。アルゴリズムには、NumPyのアルゴリズムとは少し異なる動作をする関数もあります。これらの関数には、cross、dot、inverse、determinant、eigenvalue、eigenvectorなどがあります。 5.9.3 章 の非網羅的なリストを参照してください。 これらの関数はすべて、配列/行列全体ではなく各タプルに適用されます。例えば:
>>> amatrix = algs.gradient(avector)
>>> print(repr(algs.determinant(amatrix)))
VTKArray([-1221.2732624 , -648.48272183, -3.55133937, ..., 28.2577152 ,
-629.28507693, -1205.81370163])
上記のすべてはタプルごとの情報のみを利用しており、メッシュには依存していないことに注意してください。VTKの最大の強みの1つは、そのデータ・モデルが多種多様なメッシュをサポートする一方で、そのアルゴリズムがこれらすべてのメッシュ・タイプで汎用的に機能することです。algorithmsモジュールは、この機能の一部を公開しています。その他の機能は、既存のVTKフィルタを利用することで容易に実装できます。以前に勾配を用いてベクトルと行列を生成しました。これでまた行います:
>>> avector = algs.gradient(rtdata)
>>> amatrix = algs.gradient(avector)
このような機能を使用するには、配列と関連するメッシュを含むデータセットにアクセスする必要があります。これは、dataset_adapterで ndarray のサブクラスを使用する理由の1つです。
>>> print(repr(rtdata.DataSet))
<paraview.vtk.numpy_interface.dataset_adapter.DataSet at 0x11b61e9d0>
各配列は、それを含むデータセットを指します。gradientなどの関数は、メッシュと配列を一緒に使用します。NumPyにはgradient関数も提供されています。
>>> print(repr(algs.gradient(rtdata2)))
VTKArray([[ 25.46767712, 8.78654003, 7.28477383],
[ 6.02292252, 8.99845123, 7.49668884],
[ 5.23528767, 9.80230141, 8.3005352 ],
...,
[ -6.43249083, -4.27642155, -8.30053616],
[ -5.19838905, -3.47257614, -7.49668884],
[ 13.42047501, -3.26066017, -7.28477287]])
>>> print(rtdata2.DataSet.GetClassName())
vtkUnstructuredGrid
メッシュへのアクセスを必要とする勾配とアルゴリズムは、メッシュが一様な格子、曲線格子、またはVTKのデータモデルのおかげで非構造格子のいずれであっても機能します。アルゴリズムモジュールのさまざまな関数を見てください。残りのセクションでは、これらのモジュールを使用して特定の問題を解決する方法について説明します。
6.4. 複合データセットの処理
本節では、複合データセットについて詳しく説明します。この例では、パイプラインは Sphere ソースであり、 Cone ソースは Programmable Filter への2つの入力として設定されています。
次のように、 Programmable Filter の Script にマルチブロックデータセットを作成できます。
# Let's assume inputs[0] is the output from Sphere and
# inputs[1] is the output from Cone.
mb = vtk.vtkMultiBlockDataSet()
mb.SetBlock(0, inputs[0].VTKObject)
mb.SetBlock(1, inputs[1].VTKObject)
VTKのアルゴリズムの多くは、変更することなく複合データセットで動作します。例えば:
e = vtk.vtkElevationFilter()
e.SetInputData(mb)
e.Update()
mbe = e.GetOutputDataObject(0)
print(mbe.GetClassName())
これは vtkMultiBlockDataSet を出力します。
スカラーを持つ複合データセットができたので、numpy_interface を使うことができます。 これまでと同様、単純化するために、スクリプトによって生成された出力をコード自体に埋め込み、>>> プレフィックスを使用してコードと出力を区別します。
>>> from paraview.vtk.numpy_interface import dataset_adapter as dsa
>>> mbw = dsa.WrapDataObject(mbe)
>>> print(repr(mbw.PointData.keys()))
['Normals', 'Elevation']
>>> elev = mbw.PointData['Elevation']
>>> print(repr(elev))
<paraview.vtk.numpy_interface.dataset_adapter.VTKCompositeDataArray at 0x1189ee410>
配列型がこれまで見てきたもの (`VTKArray`) と異なることに注意してください。ただし、同じように動作します。
>>> from paraview.vtk.numpy_interface import algorithms as algs
>>> print(algs.max(elev))
0.5
>>> print(algs.max(elev + 1))
1.5
各ブロックの配列には、次の方法で個別にアクセスできます。
>>> print(repr(elev.Arrays[0]))
VTKArray([ 0.5 , 0. , 0.45048442, 0.3117449 , 0.11126047,
0. , 0. , 0. , 0.45048442, 0.3117449 ,
0.11126047, 0. , 0. , 0. , 0.45048442,
0.3117449 , 0.11126047, 0. , 0. , 0. ,
0.45048442, 0.3117449 , 0.11126047, 0. , 0. ,
0. , 0.45048442, 0.3117449 , 0.11126047, 0. ,
0. , 0. , 0.45048442, 0.3117449 , 0.11126047,
0. , 0. , 0. , 0.45048442, 0.3117449 ,
0.11126047, 0. , 0. , 0. , 0.45048442,
0.3117449 , 0.11126047, 0. , 0. , 0. ], dtype=float32)
インデックス作成は少し異なることに注意してください。
>>> print(elev[0:3])
[VTKArray([ 0.5, 0., 0.45048442], dtype=float32),
VTKArray([ 0., 0., 0.43301269], dtype=float32)]
戻り値は、2つのVTKArrayで構成される複合配列です。[] 演算子は単に各配列の最初の4つの値を返しました。一般に、すべてのインデックス作成操作は、複合配列コレクション内の各VTKArrayに適用されます。アルゴリズムの場合も同様です。
>>> print(algs.where(elev < 0.5))
[(array([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,
18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34,
35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49]),),
(array([0, 1, 2, 3, 4, 5, 6]),)]
次に、Normalsと呼ばれるもう1つの配列を見てみましょう。
>>> normals = mbw.PointData['Normals']
>>> print(repr(normals.Arrays[0]))
VTKArray([[ 0.00000000e+00, 0.00000000e+00, 1.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, -1.00000000e+00],
[ 4.33883727e-01, 0.00000000e+00, 9.00968850e-01],
[ 7.81831503e-01, 0.00000000e+00, 6.23489797e-01],
[ 9.74927902e-01, 0.00000000e+00, 2.22520933e-01],
...
[ 6.89378142e-01, -6.89378142e-01, 2.22520933e-01],
[ 6.89378142e-01, -6.89378142e-01, -2.22520933e-01],
[ 5.52838326e-01, -5.52838326e-01, -6.23489797e-01],
[ 3.06802124e-01, -3.06802124e-01, -9.00968850e-01]], dtype=float32)
>>> print(repr(normals.Arrays[1]))
<paraview.vtk.numpy_interface.dataset_adapter.VTKNoneArray at 0x1189e7790>
2番目の配列が VTKNoneArray であることに注目してください。これは、vtkConeSource がNormalsを生成しないためです。配列が存在しない場合は、VTKNoneArrayをプレースホルダとして使用します。これにより、複合データセットのデータセットとVTKCompositeDataArray内のアレイ間の1対1のマッピングを維持することができます。また、特別なコードをたくさん使わなくても、アルゴリズムを並列に動かすことができます。
多くのアルゴリズムは、コレクション内の各配列に個別に適用されますが、一部のアルゴリズムはグローバルです。例えば、上で説明したように、min と max を使用します。 ブロック単位の回答を得ることが役立つ場合があります。これには _per_block アルゴリズムを使うことができます。
>>> print(algs.max_per_block(elev))
[0.5, 0.4330127]
これらは他のオペレーションと非常にうまく連携して動作します。たとえば、各ブロックの標高値を正規化する方法を次に示します。
>>> _min = algs.min_per_block(elev)
>>> _max = algs.max_per_block(elev)
>>> _norm = (elev - _min) / (_max - _min)
>>> print(algs.min(_norm))
0.0
>>> print(algs.max(_norm))
1.0
これらの機能を理解すると、単一の配列と非常によく似た複合配列を使用できるようになります。
複合データセットに関する最後の注意: numpy_interface.dataset_adapter によって提供される複合データラッパーには、複合データセットをトラバースするための便利な関数がいくつか用意されています。以下に簡単な例を示します。
>>> for ds in mbw:
>>> print(type(ds))
<class 'paraview.vtk.numpy_interface.dataset_adapter.PolyData'>
<class 'paraview.vtk.numpy_interface.dataset_adapter.PolyData'>