3次元二値画像解析¶
このチュートリアルでは3次元の画像解析を行ないます。基本的な考え方は2Dの二値画像解析と同じなので、まずそちらのチュートリアルを終わらせてからこちらにトライしてください。
ここでは50x50x50のサイズの画像を解析します。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 の画像を読み込んで 50x50x50 の白黒配列を作ります。
# 必要なライブラリの読み込み
import imageio
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
3次元配列を作る前に 0000.png を読み込んでどんな画像か確認しましょう。
image_0000 = imageio.v3.imread("data/0000.png")
image_0000.shape, image_0000.dtype, np.max(image_0000), np.min(image_0000)
((50, 50), dtype('uint8'), 255, 0)
これは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]], dtype=uint8)
plt.imshow(image_0000, "gray")
<matplotlib.image.AxesImage at 0x7f194bbdb610>
では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 # モジュールの読み込み
import homcloud.pyvistahelper as pvhelper # HomCloud側のヘルパー関数を導入
pv.set_jupyter_backend('trame') # Jupyter notebook のバックエンドを trame にする
pl = pv.Plotter()
pl.add_mesh(pvhelper.Bitmap3D(pict).threshold(0.5), show_scalar_bar=False, opacity=0.7)
pl.show()
Widget(value="<iframe src='http://localhost:35959/index.html?ui=P_0x7f19c4472b60_0&reconnect=auto' style='widt…
パーシステント図の計算¶
パーシステント図を計算します。hc.PDList.from_bitmap_distance_functionを使います。
2Dの場合と同じです。
import homcloud.interface as hc
pdlist = hc.PDList.from_bitmap_levelset(hc.distance_transform(pict, signed=True), save_to="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次)¶
この2つのペアがどのような穴を表現しているのか見てみましょう。1次のパーシステントホモロジーを見ているわけですから、 何らかのリング構造を見ているはずです。それを計算しましょう。
ここでは optimal cycle というものを使います。ポイントクラウドの解析ではoptimal volumeやvolume-optimal 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()
Widget(value="<iframe src='http://localhost:35959/index.html?ui=P_0x7f18d84f6110_1&reconnect=auto' style='widt…
次にこのリング構造がどういうパスを通っているのかを見ます。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のbirth-death pairを詳しく調べていきましょう。
ボクセルデータの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つのbirth-death pairの情報を取り出すのには
pair_nodes_in_rectangleを使いましょう。これは指定した範囲内にある birth-death pair (に対応する木構造のnode)をすべて取り出します。
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()
Widget(value="<iframe src='http://localhost:35959/index.html?ui=P_0x7f18d84aee90_2&reconnect=auto' style='widt…
opacityを変えるなどして，可視化された情報をうまく観察してください。0次のPHなので、ボクセルの内側の領域に 対応しています。実はこの2つのbirth-death pairは木構造上の親子関係にあり、一方が他方の部分集合になっています。
緑の部分の形がどうも「ぎざぎざ」で不自然だな，というようにも見えるのは stable volume というのを代わりに使うと解決される可能性があります．次のように stable_volume(1) として「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()
Widget(value="<iframe src='http://localhost:35959/index.html?ui=P_0x7f18d84c81f0_3&reconnect=auto' style='widt…
それぞれの塊がより良い感じに見えていると思います．こうすると2つのbirth-death pairの対応する領域が分離されています．stable_volumeはノイズに強い部分を取り出す仕組みなので，これは「2つのbirth-death pairの親子関係はノイズで崩れやすい」ということを意味しています．
逆解析を使うときは，このようにstable_volume(1)としておくことをお勧めしておきます．
pd2 = pdlist.dth_diagram(2)
pd2.histogram((-20.5, 20.5), 41).plot(colorbar={"type": "log"})
(4, 6)というbirth-death pairを調べます。逆解析には
bitmap_ph_trees(2)を使います。2が次数です．
phtree_2 = pdlist_with_tree.bitmap_phtrees(2)
(4, 6) という birth-death pairを調べましょう。まず pair_nodes_in_rectangle で (4, 6) という birth-death pair (に対応した tree の node) をすべて取り出します。
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個の 以下が可視化した結果です。くりぬかれた領域の奥のほうがこのbirth-death pairに対応していることがわかると思います。
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()
Widget(value="<iframe src='http://localhost:35959/index.html?ui=P_0x7f18d825afe0_4&reconnect=auto' style='widt…
こちらも 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()
Widget(value="<iframe src='http://localhost:35959/index.html?ui=P_0x7f18d84adc00_5&reconnect=auto' style='widt…
こちらはstable volumeを使ってもあまり見た目は変わらないようです．
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]]
以上でこのチュートリアルは終わりです。