この文章を読む前に周期境界でないポイントクラウドデータの解析のチュートリアルを先にやってください.
import homcloud.interface as hc # HomCloudのインターフェース
import plotly.graph_objects as go # 3次元可視化用
import homcloud.plotly_3d as p3d # 3次元可視化用
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
データを np.loadtxt
で読みこみます。
pointcloud = np.loadtxt("pointcloud.txt")
まずは pointcloud の情報を調べてみましょう.
pointcloud.shape
(1000, 3)
np.min(pointcloud, axis=0), np.max(pointcloud, axis=0)
(array([0.00026274, 0.00100671, 0.00020878]), array([0.99997803, 0.9993682 , 0.99825316]))
データは 1000x3 の配列です.これは3次元の点が1000個あることを意味します. またX,Y,Z座標の最大最小を見るとデータは$[0, 1]\times[0, 1]\times[0, 1]$という立方体に分布していることがわかります.
データからパーシステント図を計算します。hc.PDList.from_alpha_filtration
という静的メソッドを使います。
periodicity=[(0, 1), (0, 1), (0, 1)]
という引数を渡すことで$[0, 1]\times[0, 1]\times[0, 1]$という立方体での周期境界で計算します.
hc.PDList.from_alpha_filtration(pointcloud,
save_to="pointcloud-periodic.pdgm",
periodicity=[(0, 1), (0, 1), (0, 1)],
save_boundary_map=True)
PDList(path=pointcloud-periodic.pdgm)
すると pointcloud-periodic.pdgm
というファイルが生成されます。これが
パーシステント図の情報を収めたファイルです。ここから後の解析は普通のポイントクラウドの場合とほぼ同じです.
pdlist = hc.PDList("pointcloud-periodic.pdgm")
このオブジェクトはすべての次数のパーシステント図を保持している(0次から2次まで)ので、1次のパーシステント図だけ取り出します。
すべての次数を保持しているので、このクラス自体の名前はPDList
という名前です。
pd1 = pdlist.dth_diagram(1)
このメソッドの返り値は PD
class のインスタンスです。
ヒストグラムを構築し、パーシステント図をプロットします。
histogram
メソッドでヒストグラムを構築し、plot
でそれをプロットします。colorbar={"type": "log"}
を指定することでログスケールで色付けしましょう.
pd1.histogram().plot(colorbar={"type": "log"})
0 から 0.01 の部分を拡大して表示します.
pd1.histogram((0, 0.01)).plot(colorbar={"type": "log"})
このあたりの方法も周期条件でない場合と同じです.
pd1
が指すオブジェクトの births
、deaths
属性を調べることで、birth time、death timeがわかります。
pd1.births
array([0.00103596, 0.00100508, 0.00116762, ..., 0.01369992, 0.01460543, 0.01508903])
pd1.deaths
array([0.00106352, 0.00108161, 0.00119129, ..., 0.01416671, 0.0148616 , 0.0153579 ])
パーシステント図の個々の点は何らかの意味で元のポイントクラウドのリング構造、空隙構造と対応しているはずです。しかしそれがどのようなものであるのかを特定するのは実は そんなに簡単ではないです。このような解析を逆解析と呼びます。
ここでは HomCloud の強力な逆解析ツールである、optimal volumeを使いましょう。とりあえず (0.0025, 0.008) 付近にある birth-death pair の optimal volume を調べます。 optimal volume について詳しくは、大林の論文を参考にしてください。
pd.nearest_pair_to
で (0.004, 0.007) に一番近い birth-death pair を検索することができます。
pair = pd1.nearest_pair_to(0.0022, 0.0057)
これは Pair
クラスのオブジェクトで、birth timeやdeath timeといった情報を保持しています。
pair
Pair(0.0021688455616267486, 0.005663084335456451)
以下のようにすると optimal volume が計算できます。
optimal_volume = pair.optimal_volume()
このoptimal volumeを3次元可視化してみましょう。
go.Figure(data=[optimal_volume.to_plotly3d()], layout=dict(scene=p3d.SimpleScene()))