3次元点集合データ(ポイントクラウド)の解析

この文章で解説するのは

  1. 3次元の点データからパーシステント図を計算する
  2. その図をパラメータを変えながら可視化する
  3. 各birth/death timeを抽出する
  4. 逆解析を行う

の4つです。ここまではおよそどのようなポイントクラウドデータでも 共通ですので、ここまでできるようになると良いでしょう。

パーシステント図の計算

対称となるデータは pointcloud.txt というファイルです。これを解析してみましょう。

まず最初に、必要なライブラリをインポートします。

データを np.loadtxt で読みこみます。

まずは pointcloud の情報を調べてみましょう.

データは 1000x3 の配列です.これは3次元の点が1000個あることを意味します. またX,Y,Z座標の最大最小を見るとデータは$[0, 1]\times[0, 1]\times[0, 1]$という立方体に分布していることがわかります. この点を可視化してみます.

以下のようにポイントクラウドが可視化されます.

データからパーシステント図を計算します。hc.PDList.from_alpha_filtrationという静的メソッドを使います。 save_boundary_map は後で optimal volume というものを計算するのに必要です。

すると pointcloud.pdgm というファイルが生成されます。これが パーシステント図の情報を収めたファイルです。

注意

HomCloud 3.0 以降では .pdgm という拡張子を推奨しています.

パーシステント図の可視化

次に、ここからパーシステント図をプロットしてみましょう。

まず、ファイルを読み込みます。

このオブジェクトはすべての次数のパーシステント図を保持している(0次から2次まで)ので、1次のパーシステント図だけ取り出します。

すべての次数を保持しているので、このクラス自体の名前はPDListという名前です。

このメソッドの返り値は PD class のインスタンスです。

ヒストグラムを構築し、パーシステント図をプロットします。 histogramメソッドでヒストグラムを構築し、plotでそれをプロットします。

左下のほうに何かちょっと見えてあとはかなりぼんやりしています。というのは このデータは左下のほうにデータ(birth-death pair)が偏っているからです。 そこで色付けを log scale にしましょう。

基本的にはパーシステント図の点は対角線から離れるほど「意味のある」リング構造と対応し、 またY軸が大きい値になるほど大きなリング構造を表しています。 つまり(0.0025, 0.008)のあたりにある点がリング構造のなかで最もちゃんとしたリング構造を持っているものと対応していると言えそうです。

ここで注意しておくと、パーシステント図のX、Y軸は点に貼り付ける球の半径と対応していると よくある教科書には書かれています。HomCloudでは実はこれは半径の2乗が使われています。 つまり$\sqrt{0.0025}=0.5$と$\sqrt{0.008} \approx 0.09$ が実際の半径になります。これは主には内部のアルゴリズムの都合に よるものですが、各点に重みづけをしたときにはこのほうが自然に見えるという事情もあり、 HomCloudでは2乗の値が使われます。これを止めたいときはパーシステント図をhc.PDList.from_alpha_filtration で計算するときに no_squared=True という引数を付けると 半径パラメータそのものがX、Y軸に現れます。

さて、次にbirth-death pairが集中している 左下の 0.0〜0.01 のあたりを拡大して調べてみましょう。デフォルトでは ヒストグラムを構築するときにすべてのbirth-death pairが含まれるような 範囲を自動的に切り取ります。これを変更するには histogram メソッドのx_range引数(第1引数)を指定します。

ヒストグラムの解像度を変更するには histogram メソッドの x_bins引数(第2引数)を指定します。デフォルトでは 128x128でヒストグラムを描きますが、もっと細かく256x256にしてみましょう。

これらの図は matplotlib の機能を使って描画されているため、 matplotlib.pyplot.savefigで保存することができます。以下のようにすると、pointcloud-pd1.pngに図が保存されます。

2次元のPDも同様に可視化しましょう.

こちらは death (Y軸) が 0.0125 より大きいあたりにぱらぱらと birth-death pair が分布しているのが特徴的です.

数値データの読み出し

pd1pd2 が指すオブジェクトの birthsdeaths属性を調べることで、birth time、death timeがわかります。

そこで、この2つの差を見ることでlifetimeが計算できます。

lifetime のヒストグラムを表示してみましょう。

行末の ; としているのは plt.hist の返り値を無視するためのトリックです。大半の lifetime は非常に小さいことがわかります。実際対角線のそばに多くの birth-death pair がノイズ的な情報として分布していることはよくあります.

逆解析について

パーシステント図の個々の点は何らかの意味で元のポイントクラウドのリング構造、空隙構造と対応しているはずです。しかしそれがどのようなものであるのかを特定するのは実は そんなに簡単ではないです。このような解析を逆解析と呼びます。

ここでは HomCloud の強力な逆解析ツールである、optimal volumeを使いましょう。とりあえず (0.0025, 0.008) 付近にある birth-death pair の optimal volume を調べます。 optimal volume について詳しくは、大林の論文を参考にしてください。

pd.nearest_pair_to で (0.0025, 0.008) に一番近い birth-death pair を検索することができます。

このオブジェクトは Pair クラスのオブジェクトで、birth timeやdeath timeといった情報を保持しています。

以下のようにすると optimal volume が計算できます。

boundary_points メソッドで、対応するリング構造に含まれる点の座標が得られます。

このリング構造は39点からなっていることがわかります。boundaryメソッドを使うとどの点とどの点がつながっているか等もわかります。

このoptimal volumeを3次元可視化してみましょう。

Paraviewが起動し、リングが表示されます。緑色がリングを表しています。それ以外の情報も表示されていますが、それについては上に挙げた論文が参考になるでしょう。

ポイントクラウドを重ねて表示したい場合には ParaView の機能を使います.

  1. File→Open で pointcloud.txt を開く
  2. Properties Panelが表示されるので,
    • Field Delimiter Characters を適当に変更する(ここではファイルは空白区切りなので空白を入力)     - Have Headers のチェックを外す     - Apply ボタンを押す
  3. Filters→Alphabetical→Table to Points
  4. Properties Panelで
    • X ColumnをField 0     - Y ColumnをField 1     - Z ColumnをField 2    
    • Apply ボタンを押す

これで点が表示されます.

また,HomCloudのparaviewインターフェースを使って表示することも可能です.

以下のように可視化されます.3次元データを回転させて可視化してみましょう.

こんどは複数の点のoptimal volumeを同時に可視化してみます.birth < 0.004, death > 0.006 くらいの領域の birth-death pair を 以下のようにして取り出してきます.

取り出してきた birth-death pair は以下のようにして数えることができます.

では optimal volume を計算して,可視化してみましょう.54個もあるので計算には結構時間がかかります.そこで cutoff_radiusというのを指定します.これは計算ターゲットが存在する領域の大きさを指定することで計算時間を 削減するための仕組みです.だいたいこの球の内側にあるだろう,という領域の半径を指定します. この半径を予想するのは簡単ではありませんがここではデータ全体の座標が0〜1の間なので適当にその半分の 0.5くらいでいいかな,と予想してやりましょう.

2次元のPDの逆解析について

空洞のほうも同様に可視化してみましょう.birth > 0.0125 かつ death - birth > 0.0005 などがおもしろそうです.

結構大きな空洞があるようです.