3次元二値画像解析¶

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

ここでは50x50x50のサイズの画像を解析します。3次元画像の表現方法は色々ありますが、ここでは各層ごとの断面画像を用います。 つまり50x50の画像を50個用意します。 data/ ディレクトリにファイルが置かれています。ファイル名は断面の順に付けておきます。

In [1]:
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 の白黒配列を作ります。

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

3次元配列を作る前に 0000.png を読み込んでどんな画像か確認しましょう。

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

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

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

では50個の画像を読み込んで積み重ねます。

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

50x50x50の画像が完成しました。ではこれを3次元可視化しましょう。plotly.graph_objectsと homcloud.plotly_3d モジュールを使います。 詳しくはAPIドキュメントを見てください。 Bitmap3d関数をここでは使います。コメントアウトの所を有効化すると不透明度を変えることができます．

In [9]:
# 必要なモジュールを読み込む
import homcloud.plotly_3d as p3d
import plotly.graph_objects as go
In [10]:
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_distance_functionを使います。 2Dの場合と同じです。

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

以下のようにしてファイルを読み込むこともできます。

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

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

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

In [14]:
pd1 = pdlist.dth_diagram(1)
In [15]:
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 [16]:
pairs = pd1.pairs_in_rectangle(-4, -4, 4, 4)
pairs
Out[16]:
[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で計算します。

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

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

In [18]:
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()