PyVistaのデータセット#
PyVistaは、Pythonでの3D可視化とメッシュ操作を容易にする強力なライブラリです。データを適切な構造で表現することは、可視化の第一歩となります。本章では、PyVistaのデータセットの基本概念を紹介し、そのデータ構造を理解することで、3Dデータの操作や解析結果の視覚的表現をより効果的に行う方法を解説します。
from helper.python import print_subclasses
import helper.magics
import pyvista as pv
import numpy as np
from helper.pyvista import plotter_to_iframe, plot_point_indices
PyVistaは、さまざまな種類のデータを表現するために複数のデータセット（Dataset）を提供しています。データセットは、点（Point）とデータ（Data）の2つの要素で構成されます。点は互いに接続されている場合とされていない場合があり、複数の関連する点が集まることでセル（Cell）を形成します。点同士の接続は明示的または暗黙的に定義されます。データはスカラーまたはベクトルとして表現され、点またはセルに付随することができます。
本章では、いくつかの例を通じてデータセットの構造を段階的に理解していきます。PyVistaのデータセットクラスはすべてpv.DataSetを基底クラスとして継承しており、次にその継承関係を示します。
Tip
読者がデータセットの構造をより具体的に理解できるよう、PyVistaのPlotterを使用してデータセットの構造を3D図として描画します。学習の過程でこれらのプログラムを実行し、データセットへの理解を深めてください。
print_subclasses(pv.DataSet)
└──DataSet
├──Grid
│ ├──RectilinearGrid
│ └──ImageData
└──_PointSet
├──PointSet
├──PolyData
└──PointGrid
├──UnstructuredGrid
├──StructuredGrid
└──ExplicitStructuredGrid
ImageData#
最も理解しやすいデータセットはImageDataで、2次元または3次元の画像を表現するデータ構造です。単純に2次元または3次元の配列として考えることができます。配列にはデータが格納されており、点は直交する等間隔のグリッド上にあるため、点の座標を指定する必要はありません。点間の接続関係も配列内の位置によって決定されるため、接続も暗黙的です。
以下のプログラムはImageDataオブジェクトを作成し、spacing、origin、dimensions属性を設定しています。
img = pv.ImageData(
spacing=(0.1, 0.1, 0.1), origin=(0.1, 0.2, 0.3), dimensions=(3, 4, 5)
)
img
| ImageData | Information |
|---|---|
| N Cells | 24 |
| N Points | 60 |
| X Bounds | 1.000e-01, 3.000e-01 |
| Y Bounds | 2.000e-01, 5.000e-01 |
| Z Bounds | 3.000e-01, 7.000e-01 |
| Dimensions | 3, 4, 5 |
| Spacing | 1.000e-01, 1.000e-01, 1.000e-01 |
| N Arrays | 0 |
origin属性は3Dグリッドデータの開始座標を、spacing属性はX、Y、Z軸方向のグリッド間隔を、dimensions属性はX、Y、Z軸方向のグリッド数をそれぞれ示します。img.pointsを使用すると、グリッド内の各点の座標値が格納された配列を取得できます。この配列の座標は、X、Y、Z軸に沿って順に増加していきます。以下のプログラムは、最初の6点の座標を出力します。
img.points[:6]
array([[0.1, 0.2, 0.3],
[0.2, 0.2, 0.3],
[0.3, 0.2, 0.3],
[0.1, 0.3, 0.3],
[0.2, 0.3, 0.3],
[0.3, 0.3, 0.3]])
各点に対応するデータは、point_data属性に保存されており、これはDataSetAttributesオブジェクトです。
img.point_data
pyvista DataSetAttributes
Association : POINT
Active Scalars : None
Active Vectors : None
Active Texture : None
Active Normals : None
Contains arrays : None
DataSetAttributesオブジェクトには、複数のデータセットを保存できます。active_scalars属性は、各点に対応するスカラー値を格納します。NumPy配列を代入すると、自動的にPyVistaの対応する配列オブジェクトが作成され、その内容が保存されます。
例えば、以下のプログラムを実行すると、active_scalars属性はNoneからpyvista_ndarray配列に変わります。このプログラムは、各点にインデックスと同じスカラー値を追加します：
print(img.point_data.active_scalars) # データなし
img.point_data.set_scalars(np.arange(0.0, img.n_points))
print(type(img.point_data.active_scalars))
img.point_data.active_scalars
None
<class 'pyvista.core.pyvista_ndarray.pyvista_ndarray'>
pyvista_ndarray([ 0., 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., 50., 51., 52., 53., 54.,
55., 56., 57., 58., 59.])
pyvista_ndarrayは、NumPyのndarrayを継承したクラスです。そのため、NumPy配列がサポートするすべての操作を利用できます。
img.point_data.active_scalars[:2] = 10, 11
print(img.point_data.active_scalars[:5])
[10. 11. 2. 3. 4.]
pyvista_ndarrayは、VTKの配列メソッドもサポートしています。次の例では、GetName()メソッドを使って、この配列の名前を取得します。
img.point_data.active_scalars.GetName()
'scalars'
DataSetAttributesオブジェクトは、複数の配列を保持できます。また、辞書のようなインターフェースを提供しているため、他のデータを簡単に追加できます。
img.point_data.keys()
['scalars']
img.point_data["scalars"]
pyvista_ndarray([10., 11., 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., 50., 51., 52., 53., 54.,
55., 56., 57., 58., 59.])
img.point_data["zerodata"] = np.zeros(img.n_points)
img.point_data.keys()
['scalars', 'zerodata']
delキーワードを使用することで、配列名を指定して配列を削除することができます：
del img.point_data["zerodata"]
img.point_data.keys()
['scalars']
各点に対応する値は、スカラーだけでなくベクトルも指定できます。ベクトル配列は、set_vectors()メソッドを使って追加できます。以下のプログラムでは、形状が (N, 3) の NumPy 配列を vectors という名前で追加しています。ここで N は点の数を表し、ImageData オブジェクトの各点が3次元空間内のベクトルに対応します：
img.point_data.set_vectors(np.arange(0.0, img.n_points * 3).reshape(-1, 3), "vectors")
img.point_data
pyvista DataSetAttributes
Association : POINT
Active Scalars : scalars
Active Vectors : vectors
Active Texture : None
Active Normals : None
Contains arrays :
scalars float64 (60,) SCALARS
vectors float64 (60, 3) VECTORS
作成された配列は、active_vectorsとして扱われます。これは依然としてNumPy配列として操作することができます。
img.point_data.active_vectors.shape
(60, 3)
ImageDataオブジェクトでは、点の座標はspacing、origin、dimensionsなどの属性によって暗黙的に定義されます。同様に、点とセルの関係も暗黙的に定義されます。セルと点の関係は、次のグラフに示されています。セルは8つの隣接する点から構成される立方体で、図では第0セルが半透明の灰色の立方体として示されています。
plotter = pv.Plotter()
plotter.add_mesh(img.extract_all_edges(), color="black")
plotter.add_mesh(
img.glyph(
scale=False,
orient=False,
factor=0.02,
geom=pv.Sphere(theta_resolution=6, phi_resolution=6),
)
)
first_cell = img.extract_cells([0]).extract_surface(pass_pointid=False)
plotter.add_mesh(first_cell, color="gray", opacity=0.5)
plot_point_indices(plotter, first_cell, scale=0.02)
plotter.add_axes()
plotter_to_iframe(plotter)
get_cell()メソッドを使用すると、VOXEL(ボクセル)タイプのCellオブジェクトを取得できます。n_points、n_edges、n_facesなどの属性は、それぞれボクセルオブジェクトを構成する点の数、辺の数、面の数を表します。
cell = img.get_cell(0)
print(cell)
%C cell.n_points; cell.n_edges; cell.n_faces
Cell (0x1e943db1840)
Type: <CellType.VOXEL: 11>
Linear: True
Dimension: 3
N Points: 8
N Faces: 6
N Edges: 12
X Bounds: 1.000e-01, 2.000e-01
Y Bounds: 2.000e-01, 3.000e-01
Z Bounds: 3.000e-01, 4.000e-01
cell.n_points cell.n_edges cell.n_faces
------------- ------------ ------------
8 12 6
point_ids属性は、ボクセルを構成する点のインデックスリストを取得し、points属性はボクセルを構成する点の座標を取得します：
print(cell.point_ids)
cell.points
[0, 1, 3, 4, 12, 13, 15, 16]
array([[0.1, 0.2, 0.3],
[0.2, 0.2, 0.3],
[0.1, 0.3, 0.3],
[0.2, 0.3, 0.3],
[0.1, 0.2, 0.4],
[0.2, 0.2, 0.4],
[0.1, 0.3, 0.4],
[0.2, 0.3, 0.4]])
ImageDataオブジェクトのn_cells属性は、データセット内のセル数を表します。つまり、上のグラフに示された小さな立方体の数です。
img.n_cells
24
データセットは、点とセル間の関係を取得するためのさまざまなメソッドを提供しています。例えば、point_cell_ids()メソッドは特定の点を含むすべてのセルのインデックスを取得し、cell_neighbors()メソッドは特定のセルと隣接するセルを取得します。connections引数で隣接の定義を指定できます。facesの場合、指定されたセルと共通の面を持つセルが取得されます。
print(f"cells of point 3: {img.point_cell_ids(3)}")
print(f"neighbors of cell 0: {img.cell_neighbors(0, connections='faces')}")
cells of point 3: [2, 0]
neighbors of cell 0: [1, 2, 6]
各セルに対応するデータはcell_data属性に保存されており、その使用方法はpoint_dataと似ています。ここでは、その詳細について省略します。
img.cell_data
pyvista DataSetAttributes
Association : CELL
Active Scalars : None
Active Vectors : None
Active Texture : None
Active Normals : None
Contains arrays :
vtkOriginalCellIds int64 (24,)
RectilinearGrid#
ImageDataは最も単純なデータセットで、すべての点が等間隔の3Dグリッド上に配置されています。そのため、開始座標、グリッドサイズ、グリッド間隔などの情報を使って、グリッド上のすべての点の座標を計算することができます。不均一な間隔のグリッドを表現するには、次のRectilinearGridデータセットを使用します。
x = np.array([0, 3, 9, 15])
y = np.array([0, 1, 5])
z = np.array([0, 2, 3])
rgrid = pv.RectilinearGrid(x, y, z) #❶
print(f'{rgrid.dimensions = }') #❷
rgrid.point_data["scalars"] = np.arange(0.0, rgrid.n_points, dtype=np.float64) #❸
rgrid.dimensions = (4, 3, 3)
RectilinearGridはImageDataと同様に、すべての点が直交するグリッド上に配置されていますが、グリッドの分布は不均一です。そのため、X、Y、Z軸の各グリッド平面の位置を設定する必要があります。❶RectilinearGridを作成する際、x、y、zなどの3つの配列で、X軸、Y軸、Z軸に垂直な平面の位置を設定します。RectilinearGridオブジェクト内の点は、これらの平面の交点になります。❷RectilinearGridオブジェクトのdimensions属性は、これらの配列の長さに基づいています。❸最後に、point_data属性を使用して、各点に対応するデータを設定します。
次のグラフは、作成したRectilinearGridオブジェクトを表しています。半透明の長方体は、0番のセルを示しています。
plotter = pv.Plotter()
plotter.add_mesh(rgrid.extract_all_edges(), color="black")
plotter.add_mesh(
rgrid.glyph(
scale=False,
orient=False,
factor=0.5,
geom=pv.Sphere(theta_resolution=6, phi_resolution=6),
)
)
first_cell = rgrid.extract_cells([0]).extract_surface(pass_pointid=False)
plotter.add_mesh(first_cell, color="gray", opacity=0.5)
plot_point_indices(plotter, first_cell, scale=0.5)
plotter.add_axes()
plotter_to_iframe(plotter)
ImageDataオブジェクトと同様に、点のインデックスはX、Y、Z軸に沿って順次増加します。例えば：
rgrid.points[:10]
array([[ 0., 0., 0.],
[ 3., 0., 0.],
[ 9., 0., 0.],
[15., 0., 0.],
[ 0., 1., 0.],
[ 3., 1., 0.],
[ 9., 1., 0.],
[15., 1., 0.],
[ 0., 5., 0.],
[ 3., 5., 0.]])
セルと点の関係は、ImageDataと同様です：
c = rgrid.get_cell(1)
print("points of cell 1:", repr(c.point_ids))
print(c.points)
points of cell 1: [1, 2, 5, 6, 13, 14, 17, 18]
[[3. 0. 0.]
[9. 0. 0.]
[3. 1. 0.]
[9. 1. 0.]
[3. 0. 2.]
[9. 0. 2.]
[3. 1. 2.]
[9. 1. 2.]]
StructuredGrid#
RectilinearGrid よりも柔軟な構造を持つ StructuredGrid では、すべての点の座標を個別に指定する必要があります。ただし、点とセルの関係は、依然としてグリッド内での点の配置によって自動的に決定されます。以下のコードは、台形状の StructuredGrid を作成する例です。
❶ まず、NumPy の mgrid オブジェクトを使用して、3 つの配列 x、y、z を作成します。それぞれの形状はすべて (4, 5, 3) です。配列 x の値は第 0 軸に沿って、y の値は第 1 軸に沿って、z の値は第 2 軸に沿って変化します。これら 3 つの配列の対応するインデックスが、等間隔のグリッド上の点を表します。
❷ 各点の X 座標および Y 座標に、Z 座標に応じた係数を掛け算します。Z 軸の正方向に進むにつれてこの係数が徐々に小さくなるため、Z 軸に垂直な断面が異なる比率で収縮し、最終的に台形状のグリッドが形成されます。
❸ x、y、z を用いて StructuredGrid オブジェクトを作成します。StructuredGrid 内部では、これらの配列を Fortran 順（.ravel('F')）の一次元配列として扱います。
❹ 各点のインデックスを point_data["scalars"] に格納します。
x, y, z = np.mgrid[:4.0, :5.0, :3.0] #❶
x *= (4 - z) / 3 #❷
y *= (4 - z) / 3
sgrid = pv.StructuredGrid(x, y, z) #❸
sgrid.point_data["scalars"] = np.arange(0, sgrid.n_points, dtype=np.float64)
配列xのshape属性は、StructuredGridのdimensions属性として使用されます。
x.shape, sgrid.dimensions
((4, 5, 3), (4, 5, 3))
配列xは第0軸に沿って値が変化するため、一次元配列に変換すると X 軸のデータが最も早く変化します。
sgrid.points[:10]
pyvista_ndarray([[0. , 0. , 0. ],
[1.33333333, 0. , 0. ],
[2.66666667, 0. , 0. ],
[4. , 0. , 0. ],
[0. , 1.33333333, 0. ],
[1.33333333, 1.33333333, 0. ],
[2.66666667, 1.33333333, 0. ],
[4. , 1.33333333, 0. ],
[0. , 2.66666667, 0. ],
[1.33333333, 2.66666667, 0. ]])
次のグラフは、作成したグリッドの様子を示しています。灰色の半透明部分は、0番のセルを表しています。
plotter = pv.Plotter()
plotter.add_mesh(sgrid.extract_all_edges(), color="black")
plotter.add_mesh(
sgrid.glyph(
scale=False,
orient=False,
factor=0.15,
geom=pv.Sphere(theta_resolution=6, phi_resolution=6),
)
)
first_cell = sgrid.extract_cells([0]).extract_surface(pass_pointid=False)
plotter.add_mesh(first_cell, color="gray", opacity=0.5)
plot_point_indices(plotter, first_cell, scale=0.2)
plotter.add_axes()
plotter_to_iframe(plotter)
セルの形状は直方体ではありませんが、PyVista ではこのセルを Hexahedron（六面体）型として表現します。セルオブジェクトには faces や edges といった属性があり、それぞれセルを構成する面や辺を取得できます。
c = sgrid.get_cell(2)
print("cell type:", repr(c.type))
print("number_of_faces:", len(c.faces)) # セルの面数
f = c.faces[0] # 第0面を取得
print("face type:", repr(f.type)) # 各面はQuadオブジェクトで表される
print("points of face 0:", repr(f.point_ids)) # 第0面を構成する4つの点のインデックス
print("edge count of cell:", len(c.edges)) # セルの辺数
e = c.edges[0] # 第0辺を取得
print("edge type:", repr(e.type))
print("points of edge 0:", repr(e.point_ids)) # 第0辺を構成する2つの点のインデックス
cell type: <CellType.HEXAHEDRON: 12>
number_of_faces: 6
face type: <CellType.QUAD: 9>
points of face 0: [2, 22, 26, 6]
edge count of cell: 12
edge type: <CellType.LINE: 3>
points of edge 0: [2, 3]
StructuredGridを使用することで、任意の形状のグリッドを作成できます。例えば、次のプログラムは、次のグラフに示すように、中空円柱の半分を作成します。プログラムでは、まず円柱座標系で等間隔のグリッドを作成し、その後、各点の座標を直交座標系に変換します。具体的な手順については、読者自身で分析してください。
r, theta, z2 = np.mgrid[2:3:3j, -np.pi / 2 : np.pi / 2 : 6j, 0:4:7j]
x2 = np.cos(theta) * r
y2 = np.sin(theta) * r
sgrid2 = pv.StructuredGrid(x2, y2, z2)
sgrid2.point_data["scalars"] = np.arange(0, sgrid2.n_points, dtype=np.float64)
plotter = pv.Plotter()
plotter.add_mesh(sgrid2.extract_all_edges(), color="black")
plotter.add_mesh(
sgrid2.glyph(
scale=False,
orient=False,
factor=0.15,
geom=pv.Sphere(theta_resolution=6, phi_resolution=6),
)
)
plotter.add_axes()
plotter_to_iframe(plotter)
PolyData#
PolyData データセットは、一連の点、点間の接続線、および点によって構成される多角形の面で構成されています。これらの情報はすべてユーザーが設定する必要があるため、プログラムで PolyData オブジェクトを作成するのは煩雑です。しかし、pyvista の多くの関数は PolyData オブジェクトを出力します。たとえば、円錐データを作成する Cone 関数では次のように使用します：
cone = pv.Cone(resolution=4)
type(cone)
pyvista.core.pointset.PolyData
PolyData オブジェクトの points 属性は、点の座標を保存する配列です。各点の座標を簡単に確認するために、NumPy の array_str() 関数を使用してこの配列の内容を出力し、suppress_small パラメータを使って非常に小さい数を 0 として表示することができます：
print(np.array_str(cone.points, suppress_small=True))
[[ 0.5 0. 0. ]
[-0.5 0.5 0. ]
[-0.5 0. 0.5]
[-0.5 -0.5 0. ]
[-0.5 -0. -0.5]]
Cone 関数の resolution 引数を 4 に設定したため、円錐の底面は正方形になります。各点の座標から、底面が X 軸に垂直で、X = -0.5 の平面上にあり、円錐の頂点が X 軸上の X = 0.5 にあることが容易にわかります。各点間の関係は faces 属性によって決定されます：
print(cone.n_cells) # 円錐には5つの面がある
print(cone.faces)
5
[4 4 3 2 1 3 0 1 2 3 0 2 3 3 0 3 4 3 0 4 1]
faces 属性は、各面と点の関係を保存します。作成した円錐には 5 つの面があるため、PolyData オブジェクトの n_cells 属性は 5 です。faces 属性は 1 次元の整数配列を使用して、各面を構成する点のインデックスを保存します。各面を構成する点の数は異なる場合があるため、その情報も保存する必要があります。この 1 次元配列は次のように理解できます：各行は 1 つの面に対応し、コロンの前の値はその面を構成する点の数、コロンの後の一連の数字は points 属性内の各点のインデックスです。以下のデータから、円錐は 4 つの三角形と 1 つの四角形で構成されていることがわかります。
4 : 4, 3, 2, 1
3 : 0, 1, 2
3 : 0, 2, 3
3 : 0, 3, 4
3 : 0, 4, 1
次に、直接 PolyData オブジェクトを作成する方法を見てみましょう。まず、簡単な角錐の例です。
❶ 形状は (N, 3) の NumPy 配列で N 点の座標データを作成します。
❷ 各面と点の関係を表すリストを作成します。
❸ points と faces で PolyData を作成します。
❹ 各点のデータを scalars として作成します。
points = np.array(
[(1, 1, 0), (1, -1, 0), (-1, -1, 0), (-1, 1, 0), (0, 0, 2)], dtype=np.float64
) # ❶
faces = [4, 0, 1, 2, 3, 3, 4, 0, 1, 3, 4, 1, 2, 3, 4, 2, 3, 3, 4, 3, 0] #❷
p1 = pv.PolyData(points, faces) #❸
p1.point_data["scalars"] = np.arange(0.0, len(p1.points)) #❹
作成した角錐は次のグラフに示されています。図中には各点の番号が表示されています。PolyData オブジェクト内の各面は get_cell() メソッドで取得でき、図中には第 0 面と第 1 面が示されています。
plotter = pv.Plotter()
plotter.add_mesh(p1, show_edges=True)
plotter.add_mesh(
p1.glyph(
scale=False,
orient=False,
factor=0.15,
geom=pv.Sphere(theta_resolution=6, phi_resolution=6),
)
)
plotter.add_mesh(p1.extract_cells([0, 1]), color='gray', opacity=1)
plot_point_indices(plotter, p1, 0.2, index_name='scalars')
plotter.add_axes()
plotter_to_iframe(plotter)
次に、第 0 面と第 1 面の点のインデックスを取得します：
print(p1.get_cell(0).point_ids)
print(p1.get_cell(1).point_ids)
[0, 1, 2, 3]
[4, 0, 1]
以下は、PolyDataを使用して半球面を作成するプログラムです。まず、球座標系のグリッドを作成し、直交座標系の座標に変換します。x、y、zは(10, 10)の形状の２次元配列です。❶次に、三つの配列を形状が(N, 3)の配列に変換します。❷各面の点の数が同じ場合、2次元配列を使用して面と点の関係を表すことができます。ここで、第0軸の長さは面の数、第1軸の長さは各面の点の数+1です。❸各面は4つの点から構成されます。この2次元配列を直接PolyDataに渡すことができます。
N = 10
a, b = np.mgrid[0 : np.pi : N * 1j, 0 : np.pi : N * 1j]
x = np.sin(a) * np.cos(b)
y = np.sin(a) * np.sin(b)
z = np.cos(a)
points = np.c_[x.ravel('F'), y.ravel('F'), z.ravel('F')] #❶
faces = np.zeros(((N - 1) ** 2, 5), np.int32) #❷
t1, t2 = np.mgrid[: (N - 1) * N : N, : N - 1]
faces[:, 0] = 4 #❸
faces[:, 1] = (t1 + t2).ravel()
faces[:, 2] = faces[:, 1] + 1
faces[:, 3] = faces[:, 2] + N
faces[:, 4] = faces[:, 1] + N
p2 = pv.PolyData(points, faces)
p2.point_data["scalars"] = np.arange(0.0, p2.n_points)
次のグラフは、半球面の表示効果です。
plotter = pv.Plotter()
plotter.add_mesh(p2, show_edges=True)
plotter.add_axes()
plotter_to_iframe(plotter)