3次元二値画像解析

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

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

ls data/

データファイルの読み込み

まずはこの 50x50x50 の画像を読み込んで 50x50x50 の白黒配列を作ります。

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

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

image_0000 = imageio.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 0x7f5384f1adc0>

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

pict = np.stack([
    imageio.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次元可視化しましょう。homcloud.paraview_interface モジュールを使います。 詳しくはAPIドキュメントを見てください。 VoxelData クラスを主に使います。「Apply」ボタンを押すと3次元的に表示されます。Opacityを変えるなどして中まで見てください。

import homcloud.paraview_interface as pv  # モジュールの読み込み

pv.show([pv.VoxelData(pict).threshold("value", (0.5, 1.0))], wait=False)

パーシステント図の計算

パーシステント図を計算します。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]

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

pv.show([
    pv.VoxelData(pict, gui_name="voxel data").threshold("value", (0.5, 1.0)),
    hc.BitmapOptimal1Cycle.to_paraview_node_for_1cycles(optimal_1_cycles, gui_name="optimal 1-cycle"),
], wait=False)

Paraviewを色々操作して見てください。ボクセルデータのOpacityを色々変えてみる、optimal 1-cyclesの線の太さを太くする、などして 見やすい設定をさぐってください。ボクセルデータの場合には1次のbirth-death pairに対応するのはトンネル構造である、というのが この可視化を見るとわかると思います。

次にこのリング構造がどういうパスを通っているのかを見ます。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)]

パーシステント図の解析(0次)とその逆解析

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

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

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_distance_functionはPDList.from_bitmap_distance_functionと似たインターフェースで 逆解析用の情報を計算します。

この 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 というクラスのインスタンスです。次にこれらの情報を可視化しましょう。

pv.show([
    pv.VoxelData(pict).threshold("value", (0.5, 1.0)).set_opacity(0.5),
    nodes_0[0].to_pvnode().set_color((1, 0, 0)),
    nodes_0[1].to_pvnode().set_color((0, 1, 0)),
], wait=False)

Paraviewをうまく使って(Opacityを変えるなど)可視化された情報をうまく観察してください。0次のPHなので、ボクセルの内側の領域に 対応しています。実はこの2つのbirth-death pairは木構造上の親子関係にあり、一方が他方の部分集合になっています。

パーシステント図の解析(2次)とその逆解析

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

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

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に対応していることがわかると思います。

pv.show([
    pv.VoxelData(pict).threshold("value", (0.5, 1.0)).set_opacity(0.5),
    hc.BitmapPHTrees.to_pvnode_from_nodes(nodes_2).set_color((1, 0, 0)),
], wait=False)

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]]