グレイスケール画像の解析¶
ここではグレイスケール画像をパーシステントホモロジーで解析します。 ここで解説する内容は以下の通りです。白黒画像との共通点は非常に多いです。
- 画像からパーシステント図を計算する
- その図を可視化する
- テキストデータにbirth-death pairを出力する
- 基本的な逆解析(birth pixle, death pixelの出力)を行う
ここでは入力データとしてはテキストデータを使いますが、 普通の画像でも基本は同じです。
データの基本的な情報を調べる¶
ここで解析するデータは同じディレクトリにある grayscale.txt というファイルです。
数値が2次元に200x200で並んでいます。これを解析してみましょう。
import numpy as np # numpy ライブラリの読み込み
import matplotlib.pyplot as plt # 可視化のため matplitlib の読み込み
# 画像データの読み込み
pict = np.loadtxt("grayscale.txt")
pict.shape
(200, 200)
pict
array([[-1.73075897, -1.71385852, -1.69163311, ..., -1.82802151,
-1.81327553, -1.79317707],
[-1.73451332, -1.71938114, -1.69941846, ..., -1.78830493,
-1.7756238 , -1.75860343],
[-1.7414541 , -1.72962206, -1.71387809, ..., -1.72829107,
-1.72211076, -1.71168612],
...,
[-0.46998849, -0.46657936, -0.45652418, ..., 0.9402477 ,
0.95085668, 0.96084144],
[-0.46945536, -0.46409205, -0.44845549, ..., 0.89788258,
0.90681179, 0.9185174 ],
[-0.46899705, -0.4622184 , -0.44248941, ..., 0.86943194,
0.87788647, 0.8911879 ]], shape=(200, 200))
見ての通り、200x200の2次元配列です。
画像を実際に表示してみましょう。PILライブラリのPIL.Image.fromarrayを使って配列を画像に変換して表示します。jupyter notebookとipythonの組み合わせはこうして作った画像を簡単に表示できます。このあたり詳しくはそれぞれのライブラリのマニュアルを読んでください。
plt.imshow(pict, "gray")
<matplotlib.image.AxesImage at 0x7b5ab01523c0>
画素値の分布を調べましょう。
np.min(pict), np.max(pict), np.mean(pict), np.std(pict)
(np.float64(-1.8280215091732186), np.float64(3.72602997298875), np.float64(-5.684341886080802e-17), np.float64(1.0))
画素値の範囲は -1.8〜3.7、平均値はほぼ0、標準偏差は1ということがわかります。次にヒストグラムを表示しましょう。
plt.hist(pict.ravel(), range=(-4, 4), bins=128); None
複数のピークが見えています。
パーシステント図の計算¶
パーシステント図の計算は以下のようにします。
import homcloud.interface as hc
hc.PDList.from_bitmap_levelset(pict, "superlevel", save_to="grayscale.pdgm")
PDList(path=grayscale.pdgm)
hc.PDList.from_bitmap_levelsetでレベルセットパーシステントを計算します。第一引数に配列を指定します。第二引数に"superlevel"を指定することでスーパーレベルフィルトレーションを使うことを指示します。"sublevel"を代わりに与えることでサブレベルフィルトレーションを代わりに使います。save_to=...で出力ファイルを指定します。
パーシステント図の可視化¶
次に計算結果の0次のパーシステント図、つまり連結成分、島構造、を可視化しましょう。 スーパーレベルを使っているので、値の高い部分、つまり白色のピークの構造などが捉えられます。
pd = hc.PDList("grayscale.pdgm").dth_diagram(0)
pd.histogram().plot(colorbar={"type": "log"})
ちょっとヒストグラムのグリッドが細かすぎてわかりにくいようです(デフォルトは 128x128です)。グリッドを荒くしましょう。
pd.histogram(x_bins=64).plot(colorbar={"type": "log"})
この図ではbirth-death pairは図の右下のほうに来ています。 通常は birth < death なので左上のほうに現れるのですが、 スーパーレベルフィルトレーションを使うということは閾値を大きいほうから 小さいほうに下げていく過程での島の生成と消滅を見ているので birth のほうが閾値では大きくなるのです。
練習問題: birth-death pairは[0,1.5]x[0,1.5]のあたりに多く分布しているようです。
このあたりを拡大して表示しましょう。また、その画像をファイルに保存しましょう。
Birth time, death timeの表示¶
これはどのような入力でも同じです。ここではちょっと工夫して birth timeとdeath timeのペアの配列を表示してみましょう。numpyの機能を使います。
np.vstack([pd.births, pd.deaths]).transpose()
array([[ 3.72049835, 3.7197502 ],
[ 3.7163519 , 3.71581448],
[ 3.62199957, 3.61398173],
[ 1.17646489, 1.17569729],
[ 1.17565024, 1.17560317],
[ 1.17677567, 1.17489635],
[ 1.20250959, 1.17155953],
[ 1.16070165, 1.1514168 ],
[ 1.1448862 , 1.13935918],
[ 1.10316815, 1.10012108],
[ 1.09772847, 1.09706687],
[ 1.08018771, 1.07677108],
[ 1.06793501, 1.06744707],
[ 1.10463065, 1.04357627],
[ 0.88167134, 0.87714672],
[ 0.87779258, 0.87464922],
[ 0.87221495, 0.87149484],
[ 0.87378846, 0.87081415],
[ 0.84911391, 0.79822282],
[ 0.78615028, 0.78444829],
[ 0.81883079, 0.7472538 ],
[ 1.21522446, 0.74060972],
[ 0.71483203, 0.71049393],
[ 0.76913211, 0.70963995],
[ 0.71656008, 0.7091577 ],
[ 0.68458576, 0.67847277],
[ 0.58538702, 0.584013 ],
[ 0.57038128, 0.5610367 ],
[ 0.56187367, 0.5578117 ],
[ 0.57912716, 0.5437553 ],
[ 0.51183726, 0.5114067 ],
[ 0.51330706, 0.50602736],
[ 0.52295205, 0.50019607],
[ 0.43017107, 0.40577379],
[ 0.40833399, 0.36977646],
[ 0.39775193, 0.35347578],
[ 0.93587969, 0.34877654],
[ 0.33259517, 0.33257242],
[ 0.34220866, 0.32492189],
[ 0.32863671, 0.31583634],
[ 0.32598309, 0.30275411],
[ 0.33318458, 0.30085797],
[ 1.29888048, 0.23855504],
[ 0.22924011, 0.21786746],
[ 0.24130412, 0.2082839 ],
[ 0.23668959, 0.19994758],
[ 0.21192509, 0.19917353],
[ 0.29015713, 0.19807269],
[ 0.15068774, 0.06964889],
[ 0.06627242, 0.06472002],
[ 0.06179495, 0.06084787],
[ 1.03477437, 0.04809407],
[ 0.03062501, 0.02347062],
[ 0.03864399, -0.0397299 ],
[-0.10149762, -0.1128864 ],
[-0.1240269 , -0.12764439],
[-0.16523646, -0.18459246],
[-0.18987744, -0.19232931],
[-0.18800893, -0.19267836],
[-0.14841084, -0.19376773],
[-0.09791811, -0.19904108],
[-0.23042979, -0.24591754],
[-0.32730138, -0.33949227],
[-0.37450332, -0.37519826],
[-0.39194886, -0.39529467],
[-0.34620351, -0.40272859],
[-0.40412379, -0.43159528],
[-0.41042166, -0.44578582],
[-0.46482568, -0.46497538],
[-0.29606023, -0.46951701],
[-0.41233545, -0.47659197],
[-0.46199888, -0.47779584],
[-0.47777149, -0.48532429],
[-0.51266144, -0.51710601],
[-0.51813031, -0.51852311],
[-0.51736249, -0.52223673],
[-0.59607204, -0.63550974],
[-0.56677649, -0.65532262],
[-0.80571 , -0.83893901],
[-0.74883863, -0.84006845],
[-1.13211692, -1.1524443 ],
[-1.1729987 , -1.19613943],
[-1.214334 , -1.21494912],
[-1.2109894 , -1.22049237],
[-1.19800118, -1.23013554],
[-1.2360744 , -1.25939254],
[-1.21375775, -1.26177874],
[-1.29852304, -1.2998627 ],
[-1.29742081, -1.31664078],
[-1.31813159, -1.36734127],
[-1.35644958, -1.38143293],
[-1.39827581, -1.41087179],
[-1.44162259, -1.46342071],
[-1.45418752, -1.50165077],
[-1.52926759, -1.54008366],
[-1.54135963, -1.58319969],
[-1.58842435, -1.59080461],
[-1.59300271, -1.59844884],
[-1.56682931, -1.60183314],
[-1.52328 , -1.62831348],
[-1.60949457, -1.63901185],
[-1.59958114, -1.64762312],
[-1.4799324 , -1.65144922],
[-1.6797181 , -1.70790571]])
1列目がbirth time、2列目がdeath timeです。 このデータからも birth > death となっていることがわかります。
pd.essential_births
array([3.72602997])
essential_birthsで得られるものは、 death time が -∞ な birth-death pairの birth timeです。これは0次のパーシステントホモロジー特有のもので、最初に生まれて
最後まで生き延びる島を表現しています。上から閾値を変えるので
「最後まで生き延びる=-∞」となります。
普通は0次の所に一個だけこのようなものがありますが、例えばトーラス上のフィルトレーションなどを考えると他の次数でも1個以上ありえます。
対角線から離れた birth-death pairが重要な構造を表現しているわけなので、上で 見たパーシステント図から、death - birth < -0.3 となるような点の 由来を元データに戻って調べてみることにしましょう。
逆解析の手法も基本的には白黒画像の場合と同じです。 HomCloudのbirth pixel、death pixel機能を使います。島(連結成分)が生まれた/死んだときのピクセルの位置を出力します。
まずは death - birth, つまり lifetime が -0.3 より小さい pair を取り出します。
pairs = [pair for pair in pd.pairs() if pair.lifetime() < -0.3]
pairs
[Pair(1.2152244583757872, 0.7406097154584625), Pair(0.9358796930939498, 0.3487765394031042), Pair(1.2988804792640827, 0.2385550375924738), Pair(1.03477437479703, 0.04809406503892041)]
最後にこの点をplotします。以下のようにします。
hc.draw_birthdeath_pixels_2d(
pairs, pict, draw_birth=True, draw_death=True, draw_line=True
)
赤い点がbirth pixel、青い点がdeath pixel、緑の直線はbirth-deathの対応を示しています。赤い点が白い部分のピーク、青い点は2つの山の峠の位置、となっていることがわかると思います。
では一番のピークの 位置はどうなっているのでしょう？実はこれはdeath timeが-∞となるものと対応しています。 そしてdeath timeが-∞のものに関しては今回は無視しているので こうなっています。
高度な逆解析¶
二値画像の場合には BitmapPHTrees を使って高度な解析を行いました。グレイスケール画像の場合でも同様の解析を
行ってみましょう。
BitmapPHTrees.for_bitmap_levelset を使います。最初の2つの引数はPDList.from_bitmap_levelsetと同じです。
save_to=...で保存するファイルの名前を指定します。拡張子はPDのものと同じ.pdgmです．
hc.BitmapPHTrees.for_bitmap_levelset(pict, "superlevel", save_to="grayscale-tree.pdgm")
PDList(path=grayscale-tree.pdgm)
以下のようにしてファイルを読み込んでbitmap_phtreesで0次のPHの情報を取り出します。
phtrees = hc.PDList("grayscale-tree.pdgm").bitmap_phtrees(0)
実は内部的にはここで計算した結果はツリー構造でツリーのノードが各birth-death pairに対応します。 また、n次元の画像データを解析した場合、0次のPHと(n-1)次のPHに関する情報がツリー構造として計算されます。
1次の情報が欲しければbitmap_phtrees(1)とします．
lifetime < -0.3 という birth-death pair に対するノードを取り出します。
[node for node in phtrees.nodes if node.lifetime() < -0.3]
[BitmapPHTrees.Node(1.2988804792640827, 0.2385550375924738), BitmapPHTrees.Node(1.2152244583757872, 0.7406097154584625), BitmapPHTrees.Node(0.9358796930939498, 0.3487765394031042), BitmapPHTrees.Node(3.72602997298875, -inf), BitmapPHTrees.Node(1.03477437479703, 0.04809406503892041)]
death が -∞ のノードも取得できています。実際のところこれに対応する領域は画像全体になるのでこれを可視化してもあまりおもしろくありません。 そこで death が -∞ なノードは取り除きます。
nodes = [node for node in phtrees.nodes if node.lifetime() < -0.3 and node.death_time() != -np.inf]
最後にdraw_volumes_on_2d_image関数で可視化します。これは二値画像の場合と同じです。
hc.draw_volumes_on_2d_image(nodes, pict, color=(255, 0, 0), alpha=0.2, birth_position=(255,0,0), marker_size=3)
以上でこのチュートリアルは終わりです。