ボリューム分析#

データセットの体積や面積などのマスプロパティを計算します

import numpy as np

from pyvista import examples

PyVistaでデータセットの体積や面積などの質量特性を計算するには,すべてのPyVistaメッシュで pyvista.DataSetFilters.compute_cell_sizes() フィルタと pyvista.DataSet.volume プロパティを使用すると非常に簡単です.

簡単なグリッドメッシュから始めましょう.

# Load a simple example mesh
dataset = examples.load_uniform()
dataset.set_active_scalars("Spatial Cell Data")
(<FieldAssociation.CELL: 1>, pyvista_ndarray([  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
                   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
                   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
                   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
                   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
                   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
                   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
                   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
                   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
                   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
                   0.,   1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.,
                   0.,   2.,   4.,   6.,   8.,  10.,  12.,  14.,  16.,
                   0.,   3.,   6.,   9.,  12.,  15.,  18.,  21.,  24.,
                   0.,   4.,   8.,  12.,  16.,  20.,  24.,  28.,  32.,
                   0.,   5.,  10.,  15.,  20.,  25.,  30.,  35.,  40.,
                   0.,   6.,  12.,  18.,  24.,  30.,  36.,  42.,  48.,
                   0.,   7.,  14.,  21.,  28.,  35.,  42.,  49.,  56.,
                   0.,   8.,  16.,  24.,  32.,  40.,  48.,  56.,  64.,
                   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
                   0.,   2.,   4.,   6.,   8.,  10.,  12.,  14.,  16.,
                   0.,   4.,   8.,  12.,  16.,  20.,  24.,  28.,  32.,
                   0.,   6.,  12.,  18.,  24.,  30.,  36.,  42.,  48.,
                   0.,   8.,  16.,  24.,  32.,  40.,  48.,  56.,  64.,
                   0.,  10.,  20.,  30.,  40.,  50.,  60.,  70.,  80.,
                   0.,  12.,  24.,  36.,  48.,  60.,  72.,  84.,  96.,
                   0.,  14.,  28.,  42.,  56.,  70.,  84.,  98., 112.,
                   0.,  16.,  32.,  48.,  64.,  80.,  96., 112., 128.,
                   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
                   0.,   3.,   6.,   9.,  12.,  15.,  18.,  21.,  24.,
                   0.,   6.,  12.,  18.,  24.,  30.,  36.,  42.,  48.,
                   0.,   9.,  18.,  27.,  36.,  45.,  54.,  63.,  72.,
                   0.,  12.,  24.,  36.,  48.,  60.,  72.,  84.,  96.,
                   0.,  15.,  30.,  45.,  60.,  75.,  90., 105., 120.,
                   0.,  18.,  36.,  54.,  72.,  90., 108., 126., 144.,
                   0.,  21.,  42.,  63.,  84., 105., 126., 147., 168.,
                   0.,  24.,  48.,  72.,  96., 120., 144., 168., 192.,
                   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
                   0.,   4.,   8.,  12.,  16.,  20.,  24.,  28.,  32.,
                   0.,   8.,  16.,  24.,  32.,  40.,  48.,  56.,  64.,
                   0.,  12.,  24.,  36.,  48.,  60.,  72.,  84.,  96.,
                   0.,  16.,  32.,  48.,  64.,  80.,  96., 112., 128.,
                   0.,  20.,  40.,  60.,  80., 100., 120., 140., 160.,
                   0.,  24.,  48.,  72.,  96., 120., 144., 168., 192.,
                   0.,  28.,  56.,  84., 112., 140., 168., 196., 224.,
                   0.,  32.,  64.,  96., 128., 160., 192., 224., 256.,
                   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
                   0.,   5.,  10.,  15.,  20.,  25.,  30.,  35.,  40.,
                   0.,  10.,  20.,  30.,  40.,  50.,  60.,  70.,  80.,
                   0.,  15.,  30.,  45.,  60.,  75.,  90., 105., 120.,
                   0.,  20.,  40.,  60.,  80., 100., 120., 140., 160.,
                   0.,  25.,  50.,  75., 100., 125., 150., 175., 200.,
                   0.,  30.,  60.,  90., 120., 150., 180., 210., 240.,
                   0.,  35.,  70., 105., 140., 175., 210., 245., 280.,
                   0.,  40.,  80., 120., 160., 200., 240., 280., 320.,
                   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
                   0.,   6.,  12.,  18.,  24.,  30.,  36.,  42.,  48.,
                   0.,  12.,  24.,  36.,  48.,  60.,  72.,  84.,  96.,
                   0.,  18.,  36.,  54.,  72.,  90., 108., 126., 144.,
                   0.,  24.,  48.,  72.,  96., 120., 144., 168., 192.,
                   0.,  30.,  60.,  90., 120., 150., 180., 210., 240.,
                   0.,  36.,  72., 108., 144., 180., 216., 252., 288.,
                   0.,  42.,  84., 126., 168., 210., 252., 294., 336.,
                   0.,  48.,  96., 144., 192., 240., 288., 336., 384.,
                   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
                   0.,   7.,  14.,  21.,  28.,  35.,  42.,  49.,  56.,
                   0.,  14.,  28.,  42.,  56.,  70.,  84.,  98., 112.,
                   0.,  21.,  42.,  63.,  84., 105., 126., 147., 168.,
                   0.,  28.,  56.,  84., 112., 140., 168., 196., 224.,
                   0.,  35.,  70., 105., 140., 175., 210., 245., 280.,
                   0.,  42.,  84., 126., 168., 210., 252., 294., 336.,
                   0.,  49.,  98., 147., 196., 245., 294., 343., 392.,
                   0.,  56., 112., 168., 224., 280., 336., 392., 448.,
                   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
                   0.,   8.,  16.,  24.,  32.,  40.,  48.,  56.,  64.,
                   0.,  16.,  32.,  48.,  64.,  80.,  96., 112., 128.,
                   0.,  24.,  48.,  72.,  96., 120., 144., 168., 192.,
                   0.,  32.,  64.,  96., 128., 160., 192., 224., 256.,
                   0.,  40.,  80., 120., 160., 200., 240., 280., 320.,
                   0.,  48.,  96., 144., 192., 240., 288., 336., 384.,
                   0.,  56., 112., 168., 224., 280., 336., 392., 448.,
                   0.,  64., 128., 192., 256., 320., 384., 448., 512.]))

次に, .compute_cell_sizes フィルタを使用して,配列内のすべてのセルの体積を計算できます. .compute_cell_sizes フィルタは,メッシュコアのセルデータに配列を追加します.既定では,体積と面積です.

# Compute volumes and areas
sized = dataset.compute_cell_sizes()

# Grab volumes for all cells in the mesh
cell_volumes = sized.cell_data["Volume"]

また, .volume プロパティを使用してメッシュの総体積を計算することもできます.

# Compute the total volume of the mesh
volume = dataset.volume

しかし,1つのデータセットに2つのボリュームボディを残したまま閾値を設定したデータセットがあるとしたらどうでしょうか.たとえば,次のようになります.

threshed = dataset.threshold_percent([0.15, 0.50], invert=True)
threshed.plot(show_grid=True, cpos=[-2, 5, 3])
compute volume

次に,2つのボディの分類配列を割り当て,セルサイズを計算し,各ボディの体積を抽出します.以下の ボリュームの分割 でより簡単な実装があることに注意してください.

# Create a classifying array to ID each body
rng = dataset.get_data_range()
cval = ((rng[1] - rng[0]) * 0.20) + rng[0]
classifier = threshed.cell_data["Spatial Cell Data"] > cval

# Compute cell volumes
sizes = threshed.compute_cell_sizes()
volumes = sizes.cell_data["Volume"]

# Split volumes based on classifier and get the volumes
idx = np.argwhere(classifier)
hvol = np.sum(volumes[idx])
idx = np.argwhere(~classifier)
lvol = np.sum(volumes[idx])

print(f"Low grade volume: {lvol}")
print(f"High grade volume: {hvol}")
print(f"Original volume: {dataset.volume}")
Low grade volume: 518.0
High grade volume: 32.0
Original volume: 729.0

あるいは,単純に connectivity'largest' を渡し,目的のスカラー範囲を指定することで,データセットから最大のボリュームを直接抽出することもできます.

# Grab the largest connected volume within a scalar range
scalar_range = [0, 77]  # Range corresponding to bottom 15% of values
largest = threshed.connectivity('largest', scalar_range=scalar_range)

# Get volume as numeric value
large_volume = largest.volume

# Display it
largest.plot(show_grid=True, cpos=[-2, 5, 3])
compute volume

ボリュームの分割#

その代わりに,上記のようなデータセットで,接続されているすべてのボディ/ボリュームを分割したいとしたらどうでしょうか. pyvista.DataSetFilters.split_bodies() フィルタを使用して,データセット内のすべての異なる接続ボリュームを pyvista.MultiBlock データセット内のブロックに抽出することができます.たとえば,上記の例で,閾値設定されたボリュームを分割します.

# Load a simple example mesh
dataset = examples.load_uniform()
dataset.set_active_scalars("Spatial Cell Data")
threshed = dataset.threshold_percent([0.15, 0.50], invert=True)

bodies = threshed.split_bodies()

for i, body in enumerate(bodies):
    print(f"Body {i} volume: {body.volume:.3f}")
Body 0 volume: 518.000
Body 1 volume: 32.000
bodies.plot(show_grid=True, multi_colors=True, cpos=[-2, 5, 3])
compute volume

実際のデータセット#

これは,地下の河道の現実的なトレーニングデータセットです.これにより,データセットのチャネルが閾値に設定され,非常に大きなボディのそれぞれが分離され,それぞれのボリュームが計算されます.

データをロードし,チャンネルを閾値に設定します.

data = examples.load_channels()
channels = data.threshold([0.9, 1.1])

ここで,すべての異なるボディを抽出し,それらのボリュームを計算します.

bodies = channels.split_bodies()
# Now remove all bodies with a small volume
for key in bodies.keys():
    b = bodies[key]
    vol = b.volume
    if vol < 1000.0:
        del bodies[key]
        continue
    # Now lets add a volume array to all blocks
    b.cell_data["TOTAL VOLUME"] = np.full(b.n_cells, vol)

各ボディの体積をプリントします.

for i, body in enumerate(bodies):
    print(f"Body {i:02d} volume: {body.volume:.3f}")
Body 00 volume: 152667.000
Body 01 volume: 152638.000
Body 02 volume: 108024.000
Body 03 volume: 66761.000
Body 04 volume: 32520.000
Body 05 volume: 31866.000
Body 06 volume: 27857.000
Body 07 volume: 18269.000
Body 08 volume: 18238.000
Body 09 volume: 16120.000
Body 10 volume: 12550.000
Body 11 volume: 12490.000
Body 12 volume: 9861.000
Body 13 volume: 8239.000
Body 14 volume: 5166.000
Body 15 volume: 2270.000
Body 16 volume: 2085.000
Body 17 volume: 1889.000
Body 18 volume: 1548.000
Body 19 volume: 1443.000
Body 20 volume: 1150.000

そして,あらゆるボリュームを可視化します.

bodies.plot(scalars="TOTAL VOLUME", cmap="viridis", show_grid=True)
compute volume

Total running time of the script: (0 minutes 4.833 seconds)

Sphinx-Galleryによるギャラリー