In [3]:
display(Image("binary-image.png"))
ここでは白黒の画像をパーシステントホモロジーで解析します。 ここで解説する内容は以下の通りです。
基本的にやっていることは、与えられた白黒画像の黒い領域を 膨らませたり萎ませたりして島や穴の生成と消滅を調べています。 https://arxiv.org/abs/1706.10082 の論文の2.3節、特にFig.2が参考になるでしょう 。
このnotebookがあるディレクトリには binary-image.png というファイルがあります。これを解析してみましょう。まず画像を表示させます。
# 画像表示に必要なライブラリの読み込み
# IPython パッケージは jupyter notebook での各種ライブラリで
# そこから画面表示用のモジュールをインポートしている
from IPython.display import display, Image
%matplotlib inline
display(Image("binary-image.png"))
まずはimageioを使って画像を読み込みます
# imageio ライブラリの読み込み
import imageio
# 画像の読み込み
# as_gray=True とするとグレースケールで読み込まれる
image = imageio.imread("binary-image.png", as_gray=True)
画像の画素値のヒストグラムを確認しておきましょう
image.shape
(128, 128)
import matplotlib.pyplot as plt
import numpy as np
plt.hist(np.ravel(image), range=(0,256), bins=256); None
画素値は0と255の所に集中しているようです。つまり入力の時点でほぼ二値化された画像といってよいでしょう。HomCloudではTrue/Falseの値を持つ配列を二値画像とみなすので、numpyの比較演算子を使ってこのような配列に変換します。閾値は128にしておけばよいのでそうします。また、ここでは黒い領域に注目するので、128より小さい領域を取り出します。
binary_image = image < 128
以下のようにしてパーシステント図を計算します。
# HomCloudライブラリの読み出し
import homcloud.interface as hc
hc.PDList.from_bitmap_levelset(hc.distance_transform(binary_image, signed=True), save_to="binary-image.pdgm")
PDList(path=binary-image.pdgm)
すると binary-image.pdgm というファイルが生成されます。これがパーシステント図の情報を収めたファイルです。save_to=...で保存先のファイルを指定できます。signed=Trueは画像が縮んでいくのと広がっていくのと両方を考慮することを指定しています(これは
白黒画像の解析では常に有効にしておいて良いでしょう)。以下のようにしてファイルを読み出すことができます。
pdlist = hc.PDList("binary-image.pdgm")
実は hc.PDList.from_bitmap_distance_function の返り値とこれで読みこまれるパーシステント図の情報は同じものなので、ファイルを読みこまずに解析を進めることもできます。ただパーシステント図の計算は比較的高コストなので、計算結果はどこかに保存しておいて解析の歳に改めて読み込むのが良いでしょう。
次に計算結果の0次のパーシステント図、つまり連結成分、島構造、を可視化しましょう。 上の計算で黒の領域に注目するようにしたので、黒の島構造に注目します。
pdlist.dth_diagram(0).histogram().plot(colorbar={"type":"log"})
pdlist自体はすべての次数のパーシステント図の情報を持っているので、.dth_diagram(0)で0次のパーシスント図の情報を取り出し、.histogram()で二次元ヒストグラムを計算、.plot(...)でプロットします。colorbar={"type":"log"}とすることで色付けをlog scaleにします。このあたりは pointcloud での場合と同じです。
何か小さい点が図の上にぽうぽつと現れます。実はbirth time、death timeは この画像解析ルールだと整数の値しかとりません。デフォルトのヒストグラムは 128x128で計算しますから、これは細かすぎるのです。またこの図からすべての birth time, death timeは-20〜+7 くらいの範囲にあることがわかりますので、 これらをうまく表示されるようにプロットの領域や解像度を調整します。そこで 次のようにしてみましょう。
pdlist.dth_diagram(0).histogram(x_range=(-20.5, 7.5), x_bins=28).plot(colorbar={"type":"log"})
x_range=(-20.5, 7.5)でプロットする範囲を -20.5 から 7.5 に、x_bins=28ヒストグラムの分割を28x28に指定します。birth time、death timeは整数値なので、
ヒストグラムの各ビンの中心が整数になるように、±0.5の余裕を取るようにしています。
これを画像に保存するには matplotlib の savefig を使います.このあたりもポイントクラウドの場合と同じです。
pdlist.dth_diagram(0).histogram(x_range=(-20.5, 7.5), x_bins=28).plot(colorbar={"type":"log"})
plt.savefig("binary-image-pd0.png")
birth timeとdeath timeはpdlist.dth_diagram(0)で得られるオブジェクトの births, deaths を見るとわかります。
pd = pdlist.dth_diagram(0)
pd.births, pd.deaths
(array([-19., -18., -18., -18., -18., -18., -18., -18., -18., -16., -16.,
-14., -14., -14., -14., -14., -14., -14., -14., -14., -14., -14.,
-13., -11., -9., -9., -8., -8., -8., -7., -7., -6., -6.,
-6., -6., -6., -6., -6., -6., -6., -6., -6., -6., -6.,
-6., -6., -7., -6., -6., -5., -5., -5., -5., -5., -5.,
-5., -6., -5., -11., -5., -5., -5., -5., -5., -5., -5.,
-5., -5., -5., -5., -5., -8., -5., -5., -5., -5., -5.,
-5., -5., -5., -5., -5., -5., -6., -5., -5., -5., -5.,
-5., -5., -5., -5., -5., -5., -5., -5., -5., -5., -5.,
-5., -5., -5., -5., -5., -5., -6., -6., -5., -5., -6.,
-15.]),
array([-18., -17., -17., -17., -17., -17., -17., -17., -17., -15., -15.,
-13., -13., -13., -13., -13., -13., -13., -13., -13., -13., -13.,
-12., -10., -8., -8., -7., -7., -7., -6., -6., -5., -5.,
-5., -5., -5., -5., -5., -5., -5., -5., -5., -5., -5.,
-5., -5., -5., -5., -5., -4., -4., -4., -4., -4., -4.,
-4., -4., -4., -4., -4., -4., -4., -4., -4., -4., -4.,
-4., -4., -4., -4., -4., -4., -4., -4., -4., -4., -4.,
-4., -4., -4., -4., -4., -4., -4., -4., -4., -4., -4.,
-4., -4., -4., -4., -4., -4., -4., -4., -4., -4., -4.,
-4., -4., -4., -4., -4., -4., -4., -4., -4., -4., 5.,
6.]))
さて、ここでそれぞれのbirth-death pairが対応する幾何構造を調べていきましょう。 上のパーシステント図を見ると(-5,-4)の所に何かbirth death pairが 集中しているようです。これが何か調べてみましょう。HomCloudのbirth pixel、 death pixel出力機能を使ってみましょう。これは島(連結成分)が生まれた/死んだ ときのピクセルの位置を出力するものです。birth/death pixelについて詳しくは https://arxiv.org/abs/1706.10082 の論文の2.3節、特にFig.2が参考になる と思います(この論文ではbirth/death positionという名前で呼んでいます)。
pd.birth_positions
[[57, 121], [47, 112], [48, 113], [49, 114], [50, 115], [51, 116], [52, 117], [53, 118], [46, 111], [73, 111], [74, 110], [40, 112], [83, 113], [84, 114], [85, 115], [86, 116], [87, 117], [88, 118], [89, 119], [90, 120], [91, 121], [92, 122], [37, 114], [9, 112], [75, 98], [74, 97], [47, 96], [50, 94], [71, 95], [48, 34], [68, 93], [26, 55], [28, 105], [36, 30], [51, 27], [51, 36], [52, 26], [53, 25], [54, 24], [55, 23], [56, 22], [57, 21], [58, 20], [59, 12], [59, 19], [58, 11], [62, 15], [103, 31], [104, 32], [9, 35], [10, 36], [11, 37], [18, 45], [19, 46], [20, 48], [21, 49], [25, 54], [25, 103], [7, 111], [28, 17], [30, 58], [31, 23], [32, 25], [27, 15], [44, 79], [49, 1], [50, 3], [48, 0], [73, 81], [75, 82], [82, 85], [44, 32], [84, 86], [86, 87], [87, 67], [87, 88], [91, 69], [92, 48], [93, 70], [94, 49], [94, 71], [94, 95], [95, 50], [98, 74], [96, 51], [99, 53], [42, 80], [102, 79], [103, 81], [104, 83], [105, 84], [106, 85], [108, 35], [109, 36], [111, 37], [112, 38], [112, 64], [114, 39], [115, 40], [116, 67], [117, 41], [117, 68], [118, 42], [119, 98], [120, 70], [127, 105], [127, 49], [121, 71], [124, 73], [4, 29], [99, 20]]
これで birth position がわかります。death のほうも pd.death_positionsで同様のことができます。
ここで注意すべきは、座標は(y,x)の順に並んでいるということです。 これはHomCloudの仕様の問題ですが、全体の整合性を保つためどうしても こうなっているので気を付けてください。
ここまでは、birth timeの配列、death timeの配列、birth positionの配列、という形でデータを取り扱ってきましたが、逆に birth time, death time, birth position, death position を一塊のオブジェクトにしたほうが解析しやすい場合も多いでしょう。そういうときは pd.pairs()を使うとよいです。
pairs = pd.pairs()
pairs[0].birth_time(), pairs[0].death_time()
(-19.0, -18.0)
pairs[0].birth_position
[57, 121]
pairs[0].death_position
[56, 121]
さて、ここから(-5,-4)に対応するものだけを取り出しましょう。Pythonのリスト内包表記を使います。
pairs_m5_m4 = [pair for pair in pairs if pair.birth_time() == -5 and pair.death_time() == -4]
pairs_m5_m4
[Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0), Pair(-5.0, -4.0)]
確かに (-5, -4) のペアのみが取り出せているようです。最後にここで取り出した birth-death pair の birth position を入力画像の上に描画しましょう。このために便利なdraw_birthdeath_pixels_2dという関数があります。
birth_pixels_image = hc.draw_birthdeath_pixels_2d(
pairs_m5_m4, "binary-image.png", draw_birth=True
)
display(birth_pixels_image)
最初の引数で描画するbirth-death pairを、2番目の引数で下地の画像を指定します。また、birth positionのみ描画したいのでdraw_birth=Trueと指定します。draw_death=Trueなどと指定するとdeath positionのほうも描画されます。
この結果から、 黒い道路状の部分に赤い点が表示されていることがわかります。これは(-5, -4)というのは 幅の1/2が5ピクセルくらいの道状の形状に対応していることを意味しています。
ここでは0次のパーシステント図を見てきました。この図の上の各点は何らかの意味でデータの島構造を表現しているのですが、 上のbirth pixelではその「中心」とおぼしき一点のみを可視化しています。これはこれで手軽で便利なのですが 島を具体的に見ることができないのでしょうか。実は2次元画像の場合はそれを計算することが可能です。
BitmapPHTreesというクラスを利用します。まず、以下のようにして計算します。
hc.BitmapPHTrees.for_bitmap_levelset(hc.distance_transform(binary_image, signed=True), save_to="binary-image-tree.pdgm")
PDList(path=binary-image-tree.pdgm)
次のようにしてデータを読み出します。
phtrees = hc.PDList("binary-image-tree.pdgm").bitmap_phtrees(0)
pair_nodes_in_rectangleというメソッドを使って(-5, -4)に対応する情報だけ取り出します。
nodes = phtrees.pair_nodes_in_rectangle(-5, -5, -4, -4)
draw_volumes_on_2d_imageで可視化します。
引数は以下の通りです。
hc.draw_volumes_on_2d_image(nodes, "binary-image.png", color=(255, 0, 0), alpha=0.5, birth_position=(0, 255, 0))
問題の領域は薄い赤、birth positionは緑で描画されます。
黒い道路状の部分が対応しているということがよりわかりやすく可視化されています。
以上で二値画像の解析のチュートリアルは終わりです。