3次元二値画像解析¶
このチュートリアルでは、3次元二値画像解析を行います。基本的な考え方は2次元の二値画像解析と同じなので、まずそちらのチュートリアルを終えてからこのチュートリアルに進んでください。
ここでは、50x50x50ピクセルの3次元画像を解析します。
3次元画像の表現方法はいくつかありますが、このチュートリアルでは各層の断面画像を使います。つまり、50x50ピクセルの画像を50枚用意する形です。
これらのファイルはdata/ディレクトリに格納されており、ファイル名には断面の順序が反映されています。
ls data/
0000.png 0007.png 0014.png 0021.png 0028.png 0035.png 0042.png 0049.png 0001.png 0008.png 0015.png 0022.png 0029.png 0036.png 0043.png 0002.png 0009.png 0016.png 0023.png 0030.png 0037.png 0044.png 0003.png 0010.png 0017.png 0024.png 0031.png 0038.png 0045.png 0004.png 0011.png 0018.png 0025.png 0032.png 0039.png 0046.png 0005.png 0012.png 0019.png 0026.png 0033.png 0040.png 0047.png 0006.png 0013.png 0020.png 0027.png 0034.png 0041.png 0048.png
データファイルの読み込み¶
これらの画像を読み込み、そのデータから50x50x50ピクセルの白黒（二値）配列を作成します。
# ライブラリの読み込み
import imageio
import numpy as np
import matplotlib.pyplot as plt
3次元配列を構築する前に、まずは0000.pngを読み込んで、どのような画像か確認してみましょう。
image_0000 = imageio.v3.imread("data/0000.png")
image_0000.shape, image_0000.dtype, np.min(image_0000), np.max(image_0000)
((50, 50), dtype('uint8'), np.uint8(0), np.uint8(255))
50x50の配列で、0 から255の値を持っています。つまりこの画像はグレイスケール画像です。 中身を見ましょう。
image_0000
array([[255, 255, 255, ..., 255, 255, 255],
[255, 255, 255, ..., 255, 255, 255],
[255, 255, 255, ..., 255, 255, 255],
...,
[255, 255, 255, ..., 0, 0, 0],
[255, 255, 255, ..., 0, 0, 0],
[255, 255, 255, ..., 0, 0, 0]], shape=(50, 50), dtype=uint8)
plt.imshow(image_0000, "gray")
<matplotlib.image.AxesImage at 0x7ffb267c7ce0>
50個の画像をファイル名の順で読み込んで積み重ねます。
pict = np.stack([
imageio.v3.imread("data/{:04d}.png".format(n)) > 128
for n in range(50)
], axis=0)
pict.shape, pict.dtype
((50, 50, 50), dtype('bool'))
50x50x50の配列が完成しました。これを3次元で可視化します。
可視化にはPyVistaを使用し、opacity=0.7を設定して半透明にします。
import pyvista as pv # PyVistaの読み込み
import homcloud.pyvistahelper as pvhelper # HomCloud側のPyVistaヘルパー関数を導入
pv.set_jupyter_backend('trame') # PyVistaの Jupyter notebook のバックエンドを trame にする
<MagicMock name='mock()' id='140716659031216'>
pl = pv.Plotter()
pl.add_mesh(pvhelper.Bitmap3D(pict).threshold(0.5), show_scalar_bar=False, opacity=0.7)
pl.show()
<MagicMock name='mock()' id='140716145419824'>
パーシステンス図の計算¶
パーシステンス図を計算します。
hc.PDList.from_bitmap_levelsetとhc.distance_transformを使います。
2Dの場合と同様に計算します。
import homcloud.interface as hc
hc.PDList.from_bitmap_levelset(hc.distance_transform(pict, signed=True), save_to="binary-3d.pdgm")
PDList(path=binary-3d.pdgm)
ファイルを読み込みます。
pdlist = hc.PDList("binary-3d.pdgm")
パーシステンス図の解析(1次元)¶
1次元のパーシステンス図を調べていきましょう。
pd1 = pdlist.dth_diagram(1)
pd1.histogram((-15.5, 10.5), 26).plot(colorbar={"type": "log"})
(-4, 4)の位置に2つのbirth-death pairがあるようです。これを確認しましょう。
pairs = pd1.pairs_in_rectangle(-4, -4, 4, 4)
pairs
[Pair(-4.0, 4.0), Pair(-4.0, 4.0)]
確かに2つあります。
逆解析¶
これらのペアがどのような穴を表現しているのかを確認しましょう。1次元パーシステントホモロジーを扱っているので、何らかのリング構造を特定できるはずです。これを計算します。
ここではOptimal 1-cycleという機能を利用します。この機能は最短経路アルゴリズムを使って、最小の長さのリングを探索します。Pair.optimal_1_cycleで計算できます。
optimal_1_cycles = [pair.optimal_1_cycle() for pair in pairs]
まず、これを可視化しましょう。元のボクセルデータを重ねて表示します。
pl = pv.Plotter()
pl.add_mesh(pvhelper.Bitmap3D(pict).threshold(0.5), show_scalar_bar=False, opacity=0.7)
for cycle in optimal_1_cycles:
pl.add_mesh(cycle.to_pyvista_mesh(), color="red")
pl.show()
<MagicMock name='mock()' id='140716145419824'>
拡大や回転をさせたり，不透明度や色を変えて観察しましょう。 抽出された2つのリングはボクセルデータ上のトンネル構造を囲んでいます。
しかしその囲み方は多少適切ではないように見えます。最短経路による計算は近似的なものであるためです。
次に、このリング構造がどのような経路を通っているのかを確認します。Optimal1CycleForBitmap.pathを使用すると、通過するピクセルの座標を取得できます。
optimal_1_cycles[0].path()
[(49, 4, 44), (49, 4, 43), (49, 4, 42), (49, 4, 41), (49, 4, 40), (49, 4, 39), (49, 4, 38), (48, 4, 38), (48, 4, 37), (48, 4, 36), (47, 4, 36), (47, 4, 35), (47, 4, 34), (46, 4, 34), (45, 4, 34), (44, 4, 34), (43, 4, 34), (43, 4, 33), (42, 4, 33), (41, 4, 33), (40, 4, 33), (39, 4, 33), (38, 4, 33), (37, 4, 33), (36, 4, 33), (36, 5, 33), (35, 5, 33), (34, 5, 33), (34, 6, 33), (34, 7, 33), (34, 8, 33), (34, 9, 33), (34, 10, 33), (34, 11, 33), (34, 11, 34), (34, 11, 35), (34, 11, 36), (34, 11, 37), (34, 11, 38), (34, 11, 39), (34, 11, 40), (34, 11, 41), (34, 11, 42), (34, 11, 43), (34, 11, 44), (34, 11, 45), (35, 11, 45), (35, 11, 46), (36, 11, 46), (37, 11, 46), (37, 11, 47), (37, 11, 48), (38, 11, 48), (39, 11, 48), (40, 11, 48), (40, 10, 48), (41, 10, 48), (42, 10, 48), (43, 10, 48), (43, 9, 48), (44, 9, 48), (45, 9, 48), (45, 8, 48), (45, 8, 47), (46, 8, 47), (46, 7, 47), (46, 7, 46), (47, 7, 46), (47, 7, 45), (48, 7, 45), (48, 7, 44), (48, 6, 44), (48, 5, 44), (49, 5, 44), (49, 4, 44)]
optimal_1_cycles[1].path()
[(16, 42, 45), (16, 41, 45), (16, 41, 46), (15, 41, 46), (15, 40, 46), (15, 40, 47), (14, 40, 47), (14, 39, 47), (14, 39, 48), (13, 39, 48), (13, 39, 49), (12, 39, 49), (11, 39, 49), (11, 38, 49), (11, 37, 49), (10, 37, 49), (10, 36, 49), (10, 35, 49), (10, 34, 49), (9, 34, 49), (9, 33, 49), (9, 32, 49), (8, 32, 49), (8, 31, 49), (8, 30, 49), (8, 29, 49), (8, 29, 48), (8, 29, 47), (8, 29, 46), (8, 29, 45), (8, 29, 44), (8, 29, 43), (8, 29, 42), (8, 29, 41), (9, 29, 41), (9, 29, 40), (9, 29, 39), (10, 29, 39), (10, 29, 38), (10, 29, 37), (10, 29, 36), (10, 29, 35), (10, 30, 35), (10, 31, 35), (10, 32, 35), (10, 33, 35), (11, 33, 35), (11, 34, 35), (12, 34, 35), (13, 34, 35), (13, 35, 35), (14, 35, 35), (14, 36, 35), (14, 37, 35), (14, 38, 35), (14, 39, 35), (14, 40, 35), (14, 41, 35), (14, 41, 36), (15, 41, 36), (15, 41, 37), (15, 42, 37), (15, 42, 38), (16, 42, 38), (16, 42, 39), (16, 42, 40), (17, 42, 40), (18, 42, 40), (18, 42, 41), (18, 42, 42), (18, 42, 43), (18, 42, 44), (18, 42, 45), (17, 42, 45), (16, 42, 45)]
pd0 = pdlist.dth_diagram(0)
pd0.histogram((-20.5, 20.5), 41).plot(colorbar={"type": "log"})
(-16.0, -10.0) と (-15.0, -10.0) という2つのペアについて詳しく調べていきましょう。
ボクセルデータの0次元と2次元のパーシステンス図の逆解析には BitmapPHTrees というクラスを使用します（1次元には使えないため、上で説明した optimal_1_cycle を代わりに使いました）。
BitmapPHTrees.for_bitmap_levelset は PDList.from_bitmap_levelset と似たインターフェースで逆解析用の情報を計算します。
この PHTrees という名前は、ここで計算されるものが木構造であることを意味しています。n次元のビットマップデータに対して0次とn-1次の2つの情報を同時に計算し、2つの木構造を持っています。この木構造を取り出すためにはPDList.bitmap_phtrees というメソッドを使用します。
hc.BitmapPHTrees.for_bitmap_levelset(hc.distance_transform(pict, signed=True),
save_to="binary-3d-tree.pdgm")
PDList(path=binary-3d-tree.pdgm)
木構造のデータを読み込みます。
pdlist_with_tree = hc.PDList("binary-3d-tree.pdgm")
0次元パーシステントホモロジーが持つ木構造を取り出すには bitmap_phtrees(0) とします。
phtree_0 = pdlist_with_tree.bitmap_phtrees(0)
(-16.0, -10.0) と (-15.0, -10.0) という2つのペアの情報を取得するには、pair_nodes_in_rectangle を使用します。これは、指定した長方形の範囲内にあるペア（に対応する木構造のノード）をすべて取り出します。
nodes_0 = phtree_0.pair_nodes_in_rectangle(-16, -15, -10, -10)
nodes_0
[BitmapPHTrees.Node(-15.0, -10.0), BitmapPHTrees.Node(-16.0, -10.0)]
取り出したデータは BitmapPHTrees.Node クラスのインスタンスです。これらの情報を可視化しましょう。
pl = pv.Plotter()
pl.add_mesh(pvhelper.Bitmap3D(pict).threshold(0.5), show_scalar_bar=False, opacity=0.2)
pl.add_mesh(nodes_0[0].to_pyvista_mesh(), color="green")
pl.add_mesh(nodes_0[1].to_pyvista_mesh(), color="red", opacity=0.4)
pl.show()
<MagicMock name='mock()' id='140716145419824'>
不透明度 (opacity) を調整するなどして、可視化された情報をよく観察してください。これは0次元のパーシステントホモロジーなので、ボクセル内部の領域に対応しています。
実は、これら2つのペアは木構造上で親子関係にあり、一方がもう一方の部分集合になっています。
赤色の部分が不自然な形に見えるかもしれません。
これは、代わりに stable volume を使うことで改善される可能性があります。stable_volume(k) と指定すると、二値画像の場合は領域を $k$ ピクセル縮小して主要部分だけを取り出す処理を行います（グレースケール画像の場合は異なる処理です）。経験的に、$k=1$ とするとうまくいくことが多いです。
pl = pv.Plotter()
pl.add_mesh(pvhelper.Bitmap3D(pict).threshold(0.5), show_scalar_bar=False, opacity=0.2)
pl.add_mesh(nodes_0[0].stable_volume(1).to_pyvista_mesh(), color="green")
pl.add_mesh(nodes_0[1].stable_volume(1).to_pyvista_mesh(), color="red", opacity=0.4)
pl.show()
<MagicMock name='mock()' id='140716145419824'>
それぞれの塊がよりそれらしく見えていると思います。
また，こうすると2つの領域が分離されます。stable volumeはノイズに強い部分を取り出す仕組みなので、これは「2つのペアの親子関係はノイズで崩れやすい」ということを意味しています。
ボクセルデータで0次元や2次元の逆解析を使う際は、このようにstable_volume(1)としておくことをお勧めします。
pd2 = pdlist.dth_diagram(2)
pd2.histogram((-20.5, 20.5), 41).plot(colorbar={"type": "log"})
(4, 6)というペアを調べます。
bitmap_ph_trees(2) として2次元パーシステンス図上の木構造を取り出します。
phtree_2 = pdlist_with_tree.bitmap_phtrees(2)
pair_nodes_in_rectangle で (4, 6) というペア (に対応した木構造のノード) をすべて取り出します。
nodes_2 = phtree_2.pair_nodes_in_rectangle(4, 4, 6, 6)
nodes_2
[BitmapPHTrees.Node(4.0, 6.0), BitmapPHTrees.Node(4.0, 6.0)]
同じ座標のペアが2個あります。これを以下のセルで可視化します。空洞の中央に対応していることがわかります。
pl = pv.Plotter()
pl.add_mesh(pvhelper.Bitmap3D(pict).threshold(0.5), show_scalar_bar=False, opacity=0.2)
pl.add_mesh(nodes_2[0].to_pyvista_mesh())
pl.add_mesh(nodes_2[1].to_pyvista_mesh())
pl.show()
<MagicMock name='mock()' id='140716145419824'>
こちらも stable volume を使うことでコンパクトにした構造を計算できます。
pl = pv.Plotter()
pl.add_mesh(pvhelper.Bitmap3D(pict).threshold(0.5), show_scalar_bar=False, opacity=0.2)
pl.add_mesh(nodes_2[0].stable_volume(1).to_pyvista_mesh())
pl.add_mesh(nodes_2[1].stable_volume(1).to_pyvista_mesh())
pl.show()
<MagicMock name='mock()' id='140716145419824'>
volumeメソッドを使い，この領域に含まれる各ボクセルの座標を得ることができます。
nodes_2[0].volume()
[[27, 24, 25], [26, 24, 25], [27, 23, 23], [26, 24, 23], [27, 24, 23], [26, 23, 23], [26, 23, 24], [27, 24, 24], [26, 24, 24], [27, 23, 24], [27, 23, 25], [26, 24, 22], [27, 23, 26], [26, 24, 26], [27, 24, 26], [28, 25, 25], [28, 24, 25], [25, 25, 25], [26, 25, 25], [27, 25, 25], [28, 23, 25], [25, 23, 25], [26, 23, 25], [25, 24, 25], [27, 22, 25], [25, 22, 25], [26, 22, 25], [26, 25, 22], [27, 25, 22], [26, 22, 23], [25, 22, 23], [27, 22, 23], [25, 23, 23], [28, 24, 23], [25, 24, 23], [28, 23, 23], [28, 25, 23], [27, 25, 23], [26, 25, 23], [25, 25, 23], [27, 24, 21], [26, 24, 21], [27, 23, 21], [26, 23, 21], [27, 23, 22], [25, 24, 22], [26, 23, 22], [27, 24, 22], [28, 25, 24], [27, 25, 24], [26, 25, 24], [28, 23, 24], [25, 25, 24], [28, 24, 24], [25, 24, 24], [28, 22, 24], [25, 23, 24], [27, 22, 24], [26, 22, 24], [25, 22, 24], [29, 24, 25], [29, 25, 25], [24, 25, 25], [25, 26, 25], [26, 26, 25], [27, 26, 25], [28, 26, 25], [26, 24, 20], [26, 23, 20], [27, 23, 20], [27, 24, 20], [27, 22, 21], [28, 22, 21], [26, 22, 21], [25, 22, 21], [24, 25, 24], [29, 24, 24], [29, 25, 24], [25, 26, 24], [26, 26, 24], [27, 26, 24], [28, 26, 24], [28, 23, 21], [28, 25, 21], [25, 23, 21], [28, 24, 21], [25, 24, 21], [25, 25, 21], [26, 25, 21], [27, 25, 21], [29, 22, 24], [24, 23, 24], [27, 21, 24], [26, 21, 24], [25, 21, 24], [24, 22, 24], [28, 21, 24], [24, 24, 24], [26, 22, 19], [29, 23, 24], [26, 23, 19]]
以上でこのチュートリアルは終わりです。