3次元二値画像解析¶

このチュートリアルでは、3次元二値画像解析を行います。基本的な考え方は2次元の二値画像解析と同じなので、まずそちらのチュートリアルを終えてからこのチュートリアルに進んでください。

ここでは、50x50x50ピクセルの3次元画像を解析します。

3次元画像の表現方法はいくつかありますが、このチュートリアルでは各層の断面画像を使います。つまり、50x50ピクセルの画像を50枚用意する形です。

これらのファイルはdata/ディレクトリに格納されており、ファイル名には断面の順序が反映されています。

In [2]:
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ピクセルの白黒（二値）配列を作成します。

In [3]:
# 必要なライブラリの読み込み
import imageio
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

3次元配列を構築する前に、まずは0000.pngを読み込んで、どのような画像か確認してみましょう。

In [4]:
image_0000 = imageio.v3.imread("data/0000.png")
In [5]:
image_0000.shape, image_0000.dtype, np.max(image_0000), np.min(image_0000)
Out[5]:
((50, 50), dtype('uint8'), np.uint8(255), np.uint8(0))

50x50の配列で、0 から 255 の値を持っています。つまりこの画像はグレイスケール画像です。 中身を見ましょう。

In [6]:
image_0000
Out[6]:
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)
In [7]:
plt.imshow(image_0000, "gray")
Out[7]:
<matplotlib.image.AxesImage at 0x7a7dc28e6a80>
No description has been provided for this image

50個の画像をファイル名の順で読み込んで積み重ねます。

In [8]:
pict = np.stack([
    imageio.v3.imread("data/{:04d}.png".format(n)) > 128
    for n in range(50) 
], axis=0)
In [9]:
pict.shape, pict.dtype
Out[9]:
((50, 50, 50), dtype('bool'))

50x50x50の配列が完成しました。これを3次元で可視化します。 plotly.graph_objectsとhomcloud.plotly_3d モジュールを使います。 詳しくはAPIドキュメントを見てください。

コメントアウトしている行を有効化すると不透明度を変えることができます．

In [10]:
# 必要なモジュールを読み込む
import homcloud.plotly_3d as p3d
import plotly.graph_objects as go
import plotly.io
plotly.io.renderers.default = 'notebook'  # plotly Setting, 'jupyterlab' may be better.
In [11]:
fig = go.Figure(data=[p3d.Bitmap3d(pict)], layout=dict(scene=p3d.SimpleScene()))
# fig.update_traces(opacity=0.5, selector=dict(type="mesh3d"))
fig.show()

パーシステンス図の計算¶

パーシステンス図を計算します。 hc.PDList.from_bitmap_levelsetとhc.distance_transformを使います。 2Dの場合と同様に計算します。

In [12]:
import homcloud.interface as hc
In [13]:
hc.PDList.from_bitmap_levelset(hc.distance_transform(pict, signed=True), save_to="binary-3d.pdgm")
Out[13]:
PDList(path=binary-3d.pdgm)

ファイルを読み込みます。

In [14]:
pdlist = hc.PDList("binary-3d.pdgm")

パーシステンス図の解析(1次元)¶

1次元のパーシステンス図を調べていきましょう。

In [15]:
pd1 = pdlist.dth_diagram(1)
In [16]:
pd1.histogram((-15.5, 10.5), 26).plot(colorbar={"type": "log"})
No description has been provided for this image

(-4, 4) の所に2つbirth-death pairがあるようです。これを確認しましょう。

In [17]:
pairs = pd1.pairs_in_rectangle(-4, -4, 4, 4)
pairs
Out[17]:
[Pair(-4.0, 4.0), Pair(-4.0, 4.0)]

確かに2つあります。

逆解析¶

これらのペアがどのような穴を表現しているのかを確認しましょう。1次元パーシステントホモロジーを扱っているので、何らかのリング構造を特定できるはずです。これを計算します。

ここではOptimal 1-cycleという機能を利用します。この機能は最短経路アルゴリズムを使って、最小の長さのリングを探索します。Pair.optimal_1_cycleで計算できます。

In [18]:
optimal_1_cycles = [pair.optimal_1_cycle() for pair in pairs]

まず、これを可視化しましょう。元のボクセルデータを重ねて表示します。

In [19]:
fig = go.Figure(data=[
    cycle.to_plotly3d() for cycle in optimal_1_cycles
] + [
    p3d.Bitmap3d(pict, color="black", name="bitmap"),
], layout=dict(scene=p3d.SimpleScene()))
fig.update_traces(opacity=0.1, selector=dict(name="bitmap"))
fig.show()

拡大や回転をさせたり，不透明度や色を変えて観察しましょう。 抽出された2つのリングはボクセルデータ上のトンネル構造を囲んでいます。

しかしその囲み方は多少適切ではないように見えます。最短経路による計算は近似的なものであるためです。 この問題を改善することはHomCloudの今後の課題です。

次に、このリング構造がどのような経路を通っているのかを確認します。Optimal1CycleForBitmap.pathを使用すると、通過するピクセルの座標を取得できます。

In [20]:
optimal_1_cycles[0].path()
Out[20]:
[(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)]
In [21]:
optimal_1_cycles[1].path()
Out[21]:
[(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)]

0次元パーシステンス図の解析とその逆解析¶

次に、0次元のパーシステンス図を解析していきましょう。

0次元のパーシステンス図をプロットします。

In [22]:
pd0 = pdlist.dth_diagram(0)
pd0.histogram((-20.5, 20.5), 41).plot(colorbar={"type": "log"})
No description has been provided for this image

(-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 というメソッドを使用します。

In [23]:
hc.BitmapPHTrees.for_bitmap_levelset(hc.distance_transform(pict, signed=True), 
                                     save_to="binary-3d-tree.pdgm")
Out[23]:
PDList(path=binary-3d-tree.pdgm)

木構造のデータを読み込みます。

In [24]:
pdlist_with_tree = hc.PDList("binary-3d-tree.pdgm")

0次元パーシステントホモロジーが持つ木構造を取り出すには bitmap_phtrees(0) とします。

In [25]:
phtree_0 = pdlist_with_tree.bitmap_phtrees(0)

(-16.0, -10.0) と (-15.0, -10.0) という2つのペアの情報を取得するには、pair_nodes_in_rectangle を使用します。これは、指定した長方形の範囲内にあるペア（に対応する木構造のノード）をすべて取り出します。

In [26]:
nodes_0 = phtree_0.pair_nodes_in_rectangle(-16, -15, -10, -10)
nodes_0
Out[26]:
[BitmapPHTrees.Node(-15.0, -10.0), BitmapPHTrees.Node(-16.0, -10.0)]

取り出したデータは BitmapPHTrees.Node クラスのインスタンスです。これらの情報を可視化しましょう。

In [27]:
fig = go.Figure(data=[
    n.to_plotly3d() for n in nodes_0
] + [
    p3d.Bitmap3d(pict, color="black", name="bitmap")
], layout=dict(scene=p3d.SimpleScene()))
fig.update_traces(opacity=0.5, selector=dict(type="mesh3d"))
fig.update_traces(opacity=0.2, selector=dict(name="bitmap"))
fig.show()

不透明度 (opacity) を調整するなどして、可視化された情報をよく観察してください。これは0次元のパーシステントホモロジーなので、ボクセル内部の領域に対応しています。

実は、これら2つのペアは木構造上で親子関係にあり、一方がもう一方の部分集合になっています。

赤色の部分が不自然な形に見えるかもしれません。

これは、代わりに stable volume を使うことで改善される可能性があります。stable_volume(k) と指定すると、二値画像の場合は領域を $k$ ピクセル縮小して主要部分だけを取り出す処理を行います（グレースケール画像の場合は異なる処理です）。経験的に、$k=1$ とするとうまくいくことが多いです。

In [28]:
fig = go.Figure(data=[
    n.stable_volume(1).to_plotly3d() for n in nodes_0
] + [
    p3d.Bitmap3d(pict, color="black", name="bitmap")
], layout=dict(scene=p3d.SimpleScene()))
fig.update_traces(opacity=0.5, selector=dict(type="mesh3d"))
fig.update_traces(opacity=0.2, selector=dict(name="bitmap"))
fig.show()

それぞれの塊がよりそれらしく見えていると思います。

また，こうすると2つの領域が分離されます。stable volumeはノイズに強い部分を取り出す仕組みなので、これは「2つのペアの親子関係はノイズで崩れやすい」ということを意味しています。

ボクセルデータで0次元や2次元の逆解析を使う際は、このようにstable_volume(1)としておくことをお勧めします。

2次元パーシステンス図の解析とその逆解析¶

最後に2次元のパーシステンス図を解析していきましょう。

2次元パーシステンス図をプロットします。

In [29]:
pd2= pdlist.dth_diagram(2)
pd2.histogram((-20.5, 20.5), 41).plot(colorbar={"type": "log"})
No description has been provided for this image

(4, 6)というペアを調べます。 bitmap_ph_trees(2) として2次元パーシステンス図上の木構造を取り出します。

In [30]:
phtree_2 = pdlist_with_tree.bitmap_phtrees(2)

pair_nodes_in_rectangle で (4, 6) というペア (に対応した木構造のノード) をすべて取り出します。

In [31]:
nodes_2 = phtree_2.pair_nodes_in_rectangle(4, 4, 6, 6)
nodes_2
Out[31]:
[BitmapPHTrees.Node(4.0, 6.0), BitmapPHTrees.Node(4.0, 6.0)]

同じ座標のペアが2個あります。これを以下のセルで可視化します。空洞の中央に対応していることがわかります。

In [32]:
fig = go.Figure(data=[
    n.to_plotly3d() for n in nodes_2
] + [
    p3d.Bitmap3d(pict, color="black", name="bitmap")
], layout=dict(scene=p3d.SimpleScene()))
fig.update_traces(opacity=0.2, selector=dict(name="bitmap"))
fig.show()

こちらも stable volume を使うことでコンパクトにした構造を計算できます。

In [33]:
fig = go.Figure(data=[
    n.stable_volume(1).to_plotly3d() for n in nodes_2
] + [
    p3d.Bitmap3d(pict, color="black", name="bitmap")
], layout=dict(scene=p3d.SimpleScene()))
fig.update_traces(opacity=0.2, selector=dict(name="bitmap"))
fig.show()

volumeメソッドを使い，この領域に含まれる各ボクセルの座標を得ることができます。

In [34]:
nodes_2[0].volume()
Out[34]:
[[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]]

以上でこのチュートリアルは終わりです。