3次元周期境界点集合データの解析

この文章を読む前に周期境界でないポイントクラウドデータの解析のチュートリアルを先にやってください．

パーシステント図の計算

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

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

%load_ext autoreload
%autoreload 2
import homcloud.interface as hc  # HomCloudのインターフェス
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import pyvista as pv  # PyVista (3D可視化)
pv.set_jupyter_backend('trame')

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

データを 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 というファイルが生成されます。これが パーシステント図の情報を収めたファイルです。ここから後の解析は普通のポイントクラウドの場合とほぼ同じです．

注意

• ポイントクラウドデータが小さすぎる場合には周期境界がうまく動作しない場合があります．その場合には周期的なコピーをX，Y，Z方向にそれぞれ何倍か複製すると良いでしょう．
• 今のところ周期境界条件は立方体，つまりX,Y,Zの幅が等しい場合にしか対応していません
• 周期境界でない場合でも同様ですが，ポイントクラウドが整いすぎている，例えば3点が同じ直線にのったり4点が同じ平面にのったりする場合には計算がうまくいかない場合があります．結晶データの場合には特この状況が発生します．解決策としては各点の座標に小さなノイズを加えてください．

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

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

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

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次元可視化してみましょう。

pl = pv.Plotter()
pl.add_mesh(optimal_volume.to_pyvista_boundary_mesh())
pl.show()

Widget(value='<iframe src="http://localhost:41573/index.html?ui=P_0x7d9680187c10_3&reconnect=auto" class="pyvi…

Paraviewが起動し、リングが表示されます。一部リングが変な場所に飛んでいるように見えますが，それは周期境界を跨いでいるのでこう見えます

これの見た目を改善するためには，以下のようにします．adjust_peridoic_boundary=(0.3, None)で，周期境界の「端」から±0.3の割合の所にある頂点に周期境界を気にさせるようにします． 周期境界も表示するようにします．

pl = pv.Plotter()
pl.add_mesh(optimal_volume.to_pyvista_boundary_mesh(adjust_periodic_boundary=(0.3, None)))
# bounds=(0, 1, ...) は周期境界条件の範囲．
pl.add_mesh(pv.Box(bounds=(0, 1, 0, 1, 0, 1), level=2), style='wireframe', color="red")
pl.show()

Widget(value='<iframe src="http://localhost:41573/index.html?ui=P_0x7d96801879a0_4&reconnect=auto" class="pyvi…

これでも周期境界を通り抜けてループをなしていることはわかりますが，周期境界で線を複製してもう少し見やすくしてみましょう． adjust_periodic_boundary=(0.3, 0.3)とします．2つめの0.3は周期境界の「端」から±0.3の割合の所にある線や頂点を複製します．

pl = pv.Plotter()
pl.add_mesh(optimal_volume.to_pyvista_boundary_mesh(adjust_periodic_boundary=(0.3, 0.3)))
# bounds=(0, 1, ...) は周期境界条件の範囲．
pl.add_mesh(pv.Box(bounds=(0, 1, 0, 1, 0, 1), level=2), style='wireframe', color="red")
pl.show()

Widget(value='<iframe src="http://localhost:41573/index.html?ui=P_0x7d96801859f0_5&reconnect=auto" class="pyvi…

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

optimal_volume.boundary_points()

[[0.14317939281543823, 0.08641612275154564, 0.2752064138682909],
 [0.025843202045195746, 0.4412407932483413, 0.3131146662301325],
 [0.9682911290945488, 0.3916907506112238, 0.2609586940046116],
 [0.28595062229479495, 0.6349264952707766, 0.5122748612177933],
 [0.12413139370195103, 0.9955678640871992, 0.28471679969247143],
 [0.48209910273061174, 0.6505350114882772, 0.5486480098843861],
 [0.46981160438464675, 0.5219067759595548, 0.4412818861069293],
 [0.4643529555106557, 0.38891409104439234, 0.37758039548286215],
 [0.9989527470705499, 0.25841793345924924, 0.6440558456175349],
 [0.39009567851528515, 0.5793629250289116, 0.5797318835031173],
 [0.13409028719087068, 0.6354942001243541, 0.40887284917547795],
 [0.002104056512585717, 0.29622758604183674, 0.5093623894255657],
 [0.3948113683540698, 0.5084120755128144, 0.4651290775456587],
 [0.9917615011994608, 0.33788429342340687, 0.3191661932756794],
 [0.07478379187666995, 0.1318144315592794, 0.6167927456566142],
 [0.32192148499694273, 0.27568117312184415, 0.27938277535368994],
 [0.37077607976531557, 0.18699903517410632, 0.17464767863151776],
 [0.920248389750528, 0.326293124436118, 0.36752925977801376],
 [0.09398645654180116, 0.5112288078274436, 0.4101723844530256],
 [0.3666635611645276, 0.2424199961546708, 0.24272150619865196],
 [0.4887373110939709, 0.5596049391459755, 0.6878955598748427],
 [0.1779098717306835, 0.16300248669786055, 0.2618915483858949],
 [0.06800441912888222, 0.1305283961003515, 0.5058547294424975],
 [0.4619450486593817, 0.6017627787429797, 0.47949711406702566],
 [0.41898602615433533, 0.44362916692289767, 0.3710002303374377],
 [0.0009620214867345211, 0.32624210024523104, 0.3299785356144932],
 [0.9459593549121383, 0.30824172357662616, 0.45444415865244425],
 [0.41872842706632685, 0.3117355628197983, 0.36357591323139904],
 [0.5624614375978, 0.6782509559928688, 0.5782499596280716],
 [0.21146515573678626, 0.7087412818645273, 0.3981480720640632],
 [0.22896706291545432, 0.19008820307844132, 0.2038858832397007],
 [0.006270880292196468, 0.10979896194917471, 0.5705032405071087],
 [0.3806980910779977, 0.4089441782551778, 0.43539259723846135],
 [0.28890235055155566, 0.19755580527182426, 0.1506399166240383],
 [0.5978877785048142, 0.6208192117413165, 0.6243499650652846],
 [0.4209916164127331, 0.550325389811247, 0.6350627818165291],
 [0.1655314819171636, 0.5632778523877888, 0.42906124661705347],
 [0.5168428169854183, 0.5772551474145238, 0.6098838405615449],
 [0.36385197619445886, 0.6046102044115498, 0.5050128183567079],
 [0.1773422292280451, 0.6978490974425432, 0.4591252721946766],
 [0.25272442784352755, 0.17629982556763446, 0.24425649532600746],
 [0.9636268547804409, 0.24647613378352018, 0.538180598867259],
 [0.02463949624447459, 0.5129939701650063, 0.3685529375160984],
 [0.1180617075890974, 0.03607810750091334, 0.3767880296175038],
 [0.3373945662185728, 0.46415106121781036, 0.4856613355037861],
 [0.03223409787023468, 0.08546385914763066, 0.4515319922985789],
 [0.03570453518233241, 0.17677934680678953, 0.6507549326106109],
 [0.13176352745648068, 0.9611707326966676, 0.36717940449927045],
 [0.1736042250190336, 0.012175214910754906, 0.24208220916104362],
 [0.3659135246377513, 0.33581293452045724, 0.30788962330093506],
 [0.2634758723979528, 0.6406149467861669, 0.4258730692776027],
 [0.05540748469500101, 0.006777502709981453, 0.4234521764308643],
 [0.9634957773822149, 0.19393250568789067, 0.5995595844093459]]

このリング構造はかなり沢山の点からなっていることがわかります。この座標出力も周期境界は考慮せずに0から1の範囲内で出力しています．この逆解析結果を使ってさらなる解析をしたい場合には適当に座標を調整する必要があるでしょう． boundaryメソッドを使うとどの点とどの点がつながっているか等もわかります。