3次元二値画像解析

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

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

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

入力データの可視化

解析する前にまずこの画像を3次元で可視化してみましょう。ここでは paraview というソフトウェアを使います。 paraview は VTK フォーマットの入力を受け付けるため、画像をこのフォーマットに変換します。

python3 -m homcloud.pict.pict3d_vtk -o bin3d.vtk -T pictures3d data/*.png

次のコマンドでparaviewを起動します。

paraview bin3d.vtk

最初は何も表示されないと思います。まずは左側ペインにある「Apply」というボタンを押してください。 すると赤青二色で塗られたキューブが表示されると思います。これは画像データが0(青)と255(赤)の二値画像であることから このように表示されます。まずはこれを回転させたり拡大縮小させたりして操作に慣れてください。

ただ、このままだと内部の様子が見れません。そこで paraview の フィルタという機能を使って 255 の値の部分だけを表示しましょう。 メニューの Filters > Common > Threshold というのを選んでください。すると左側のペイン「Pipeline browser」の所の bin3d.vtk の下に Threshold1 というのが表示されます。 この Threshold1 をクリックしてその後その下側の「Properties」ペインの中にある「Minimum」というスクロールバーの値を 大きくします。適当に128くらいにしておけばOKです(すると128より大きい所だけ切り取られるため)。その後 Apply ボタンを押すと 今いじった部分が反映されます。これで中身が見えるはずです。

ここで「Minimum」を0、「Maximum」を128などとしてApplyをすると、反転した領域が表示されます。paraviewは後でも使うので ちょっと慣れておいてください。

試し終わったら paraview を閉じてください。

パーシステント図の計算

パーシステント図を計算します。homcloud.pict.binarize_nd を用います。オプションの意味は基本的に2Dの場合と同じですが、 復習も兼ねて以下に説明します。

• -T pictures3d というので入力データが2次元画像の積み重ねで3次元画像を表現していることを指定します
  • 他にも npy でnumpyのnpyフォーマットが指定できます。ファイル1個で3次元画像を表現できるのでこれも便利でしょう
• -m white-base というので白の領域に注目することを指定
• -t 128 というので白黒の閾値を128に決める(この入力画像では黒=0、白=255なのでどんな値でも あまり問題にはなりません)
• -o bin3d.pdgmで出力ファイルを指定
• -s で画像が縮んでいくのと広がっていくのと両方を考慮することを指定(これは 白黒画像の解析では常に有効にしておいて良いでしょう)

python3 -m homcloud.pict.binarize_nd -t 128 -m white-base -T pictures3d -o bin3d.pdgm -s data/*.png

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

python3 -m homcloud.plot_PD -d 1 -x "[-15.5:10.5]" -X 26 -l --aspect equal -o bin3d-1.png bin3d.pdgm 
display < bin3d-1.png

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

python3 -m homcloud.dump_diagram -d 1 -S yes bin3d.pdgm | grep "^-4.0 4.0"

-4.0 4.0 (49,4,44) (49,14,35)
-4.0 4.0 (16,42,45) (14,36,42)

逆解析

この2つのペアがどのような穴を表現しているのか見てみましょう。そのために、このdeathの場所を可視化します。

1次の場合、deathの場所が穴(3次元上の通路状の形)の場所になるのでそれを調べます。以下のコマンドで paraview で読みこめる VTKファイルを生成します。

オプションの意味などは以下の通りです。

• -d 1 1次のPDを対象にしている
• -f "birth == -4" -f "death == 4" birth がちょうど -4、 death がちょうど 4 の pair を対象にする
• -D death の場所を可視化
• --vtk-output death_pos.vtk 出力ファイルを death_pos.vtk に指定

python3 -m homcloud.view_index_pict3d -d 1 -f "birth == -4" -f "death == 4" -D --vtk-output death_pos.vtk bin3d.pdgm

これを paraview で可視化します。以下のような手順です。

1. paraview を起動する (次のセルを実行) (ただし、MacOSの場合は代わりにopen /Applications/Paraview-5.6.0.appなどとして実行できます。)

2. 入力データを読みこみ表示します。以下の手順です。

   • メニューの File > Open... で 最初のほうで作った bin3d.vtk を読み込み
   • 左下ペイン「Properties」の Apply ボタンを押す
   • メニューの Filters > Common > Threshold を選択
   • 左下ペイン「Properties」の「Minimum」を128くらいに変更
   • 左下ペイン「Properties」の Apply ボタンを押す

3. death_pos.vtk を読み込んで表示

   • メニューの File > Open... で 最初のほうで作った death_pos.vtk を読み込み
   • 左下ペイン「Properties」の Apply ボタンを押す

すると元データと death の場所が表示されます。death の場所は青いボックスで表示されるので データを回転させたりして探してください。 見つからないときは左上ペインの「Threshold1」の隣にある目玉ボタンを押すと元データを 表示っしたり消したりできるのでそれも使ってみてください。 見つかったならば、確かにそれが通路状の構造の所に death のピクセルがある、ということがわかると思います。

paraview

他にも気になる birth-death pair があれば同様に可視化してみましょう。

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