以下のようにポイントクラウドが可視化されます．
この文章で解説するのは
の4つです。ここまではおよそどのようなポイントクラウドデータでも 共通ですので、ここまでできるようになると良いでしょう。
import homcloud.interface as hc # HomCloudのインターフェス
import homcloud.paraview_interface as pv # 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]$という立方体に分布していることがわかります． この点を可視化してみます．
pv.show([pv.PointCloud.from_array(pointcloud)])
以下のようにポイントクラウドが可視化されます．
データからパーシステント図を計算します。hc.PDList.from_alpha_filtrationという静的メソッドを使います。
save_boundary_map は後で optimal volume というものを計算するのに必要です。
hc.PDList.from_alpha_filtration(pointcloud,
save_to="pointcloud.pdgm",
save_boundary_map=True)
PDList(path=pointcloud.pdgm)
すると pointcloud.pdgm というファイルが生成されます。これが
パーシステント図の情報を収めたファイルです。
HomCloud 3.0 以降では .pdgm という拡張子を推奨しています．
pdlist = hc.PDList("pointcloud.pdgm")
このオブジェクトはすべての次数のパーシステント図を保持している(0次から2次まで)ので、1次のパーシステント図だけ取り出します。
すべての次数を保持しているので、このクラス自体の名前はPDListという名前です。
pd1 = pdlist.dth_diagram(1)
このメソッドの返り値は PD class のインスタンスです。
ヒストグラムを構築し、パーシステント図をプロットします。
histogramメソッドでヒストグラムを構築し、plotでそれをプロットします。
pd1.histogram().plot()
左下のほうに何かちょっと見えてあとはかなりぼんやりしています。というのは このデータは左下のほうにデータ(birth-death pair)が偏っているからです。 そこで色付けを log scale にしましょう。
pd1.histogram().plot(colorbar={"type": "log"})
基本的にはパーシステント図の点は対角線から離れるほど「意味のある」リング構造と対応し、 また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引数）を指定します。
pd1.histogram((0, 0.01)).plot(colorbar={"type": "log"})
ヒストグラムの解像度を変更するには histogram メソッドの
x_bins引数（第2引数）を指定します。デフォルトでは
128x128でヒストグラムを描きますが、もっと細かく256x256にしてみましょう。
pd1.histogram((0, 0.01), 256).plot(colorbar={"type": "log"})
これらの図は matplotlib の機能を使って描画されているため、
matplotlib.pyplot.savefigで保存することができます。以下のようにすると、pointcloud-pd1.pngに図が保存されます。
pd1.histogram((0, 0.01), 256).plot(colorbar={"type": "log"})
plt.savefig("pointcloud-pd1.png")
2次元のPDも同様に可視化しましょう．
pd2 = pdlist.dth_diagram(2)
pd2.histogram().plot(colorbar={"type": "log"})
こちらは death (Y軸) が 0.0125 より大きいあたりにぱらぱらと birth-death pair が分布しているのが特徴的です．
pd1 や pd2 が指すオブジェクトの births、deaths属性を調べることで、birth time、death timeがわかります。
pd1.births
array([0.00140906, 0.00139591, 0.00134116, ..., 0.01572219, 0.01709659,
0.01871588])
pd1.deaths
array([0.0015926 , 0.00164399, 0.00170775, ..., 0.0157662 , 0.01717465,
0.01880101])
そこで、この2つの差を見ることでlifetimeが計算できます。
pd1.deaths - pd1.births
array([1.83534435e-04, 2.48078852e-04, 3.66588805e-04, ...,
4.40061409e-05, 7.80638251e-05, 8.51306753e-05])
lifetime のヒストグラムを表示してみましょう。
plt.hist(pd1.deaths - pd1.births, bins=100);
行末の ; としているのは 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 = pd1.nearest_pair_to(0.0025, 0.008)
このオブジェクトは Pair クラスのオブジェクトで、birth timeやdeath timeといった情報を保持しています。
pair
Pair(0.0025623095880009357, 0.008135892057712315)
以下のようにすると optimal volume が計算できます。
optimal_volume = pair.optimal_volume()
boundary_points メソッドで、対応するリング構造に含まれる点の座標が得られます。
optimal_volume.boundary_points()
[[0.24041004740553085, 0.3620263927982651, 0.3938519814000714], [0.017231099351335488, 0.10413818329710312, 0.6827502061744146], [0.07850470821281752, 0.31832538225761675, 0.2726051865809148], [0.26866264289512143, 0.5145750814813962, 0.5502132034822121], [0.1255383953826733, 0.0800900311964785, 0.7190957020829821], [0.1252375238953729, 0.20843145021582987, 0.37381081780757186], [0.27896903483481006, 0.2900038108788685, 0.3853792732744701], [0.11968719444459341, 0.46213675229879714, 0.6086047255929676], [0.10439142042970573, 0.12004955606294554, 0.3999294686167837], [0.07478379187666995, 0.1318144315592794, 0.6167927456566142], [0.241996106473225, 0.4429189262369284, 0.4639170406136188], [0.2422121988843836, 0.4155487622456445, 0.7075748545403898], [0.33836973880327526, 0.5114239253242512, 0.7178588306122617], [0.19443369660602383, 0.2643008696498902, 0.7685550494906768], [0.33213749494421396, 0.4315334142137247, 0.7323967387271755], [0.41898602615433533, 0.44362916692289767, 0.3710002303374377], [0.12195279248515212, 0.003239458152972219, 0.6061209669659768], [0.21608561791517278, 0.2918069725141198, 0.3219258262415402], [0.1894939883612946, 0.03095999913576697, 0.6762565272706147], [0.0708147888171643, 0.05566245889432686, 0.6142169362499716], [0.39171929078549805, 0.4405548898038244, 0.288435385498083], [0.21062788160417378, 0.11633529728096403, 0.5162505592185838], [0.3166872139448015, 0.38348372951946463, 0.33415602643779374], [0.3806980910779977, 0.4089441782551778, 0.43539259723846135], [0.17739556404702383, 0.10003158291956527, 0.43755332223623955], [0.3171296008248521, 0.531306547414789, 0.6355063863352421], [0.1563951908464587, 0.13629542421604657, 0.594303189697284], [0.15814254898557367, 0.44579878013749064, 0.4799628611505681], [0.17887993848973416, 0.42554668220768377, 0.6217943706156428], [0.09845847975583133, 0.26404248130005237, 0.7778474474038211], [0.3373945662185728, 0.46415106121781036, 0.4856613355037861], [0.08346167978435703, 0.26153343284444097, 0.3241366037151565], [0.03570453518233241, 0.17677934680678953, 0.6507549326106109], [0.08561651425017891, 0.42297356026275923, 0.528205350556188], [0.12886238414379725, 0.17512618876757002, 0.7438074952407614], [0.16660019888277078, 0.33716865816694763, 0.30111066188327684], [0.18595842828877163, 0.502526980520708, 0.6007699835223476], [0.32035651808142795, 0.4751033316678963, 0.34029734006837487], [0.18678772331743787, 0.34044366488617717, 0.7226556446102448]]
このリング構造は39点からなっていることがわかります。boundaryメソッドを使うとどの点とどの点がつながっているか等もわかります。
このoptimal volumeを3次元可視化してみましょう。
pv.show([optimal_volume], wait=False)
Paraviewが起動し、リングが表示されます。緑色がリングを表しています。それ以外の情報も表示されていますが、それについては上に挙げた論文が参考になるでしょう。
ポイントクラウドを重ねて表示したい場合には ParaView の機能を使います．
これで点が表示されます．
また，HomCloudのparaviewインターフェースを使って表示することも可能です．
pv.show([
pv.PointCloud.from_array(pointcloud, gui_name="Point Cloud"),
optimal_volume.to_pvnode("Optimal Volume (0.0025, 0.008)"),
], wait=False)
以下のように可視化されます．3次元データを回転させて可視化してみましょう．
こんどは複数の点のoptimal volumeを同時に可視化してみます．birth < 0.004, death > 0.006 くらいの領域の birth-death pair を 以下のようにして取り出してきます．
pairs = [pair for pair in pd1.pairs() if pair.birth_time() < 0.004 and pair.death_time() > 0.006]
取り出してきた birth-death pair は以下のようにして数えることができます．
len(pairs)
54
では optimal volume を計算して，可視化してみましょう．54個もあるので計算には結構時間がかかります．そこで
cutoff_radiusというのを指定します．これは計算ターゲットが存在する領域の大きさを指定することで計算時間を
削減するための仕組みです．だいたいこの球の内側にあるだろう，という領域の半径を指定します．
この半径を予想するのは簡単ではありませんがここではデータ全体の座標が0〜1の間なので適当にその半分の
0.5くらいでいいかな，と予想してやりましょう．
%%time
optimal_volumes = [pair.optimal_volume(cutoff_radius=0.5) for pair in pairs]
CPU times: user 37.2 s, sys: 1.24 s, total: 38.4 s Wall time: 45.4 s
pv.show(
[v.to_pvnode(gui_name="Volume for ({:.5f}, {:.5f})".format(v.birth_time(), v.death_time())) for v in optimal_volumes]
+ [pv.PointCloud.from_array(pointcloud, gui_name="Point cloud")],
wait=False
)
空洞のほうも同様に可視化してみましょう．birth > 0.0125 かつ death - birth > 0.0005 などがおもしろそうです．
pairs = [pair for pair in pd2.pairs() if pair.birth_time() > 0.0125 and pair.lifetime() > 0.0005]
optimal_volumes = [pair.optimal_volume() for pair in pairs]
pv.show(optimal_volumes + [pv.PointCloud.from_array(pointcloud, gui_name="Pointcloud")], wait=False)
結構大きな空洞があるようです．
