3次元点群データの周期境界条件下における解析¶

このチュートリアルでは、3次元点群データに対して周期境界条件を適用し、パーシステントホモロジー解析を行います。

周期境界条件を用いない場合の解析については既に習得済みであることを前提としています。まだそちらを学習していない場合は、先にそちらのチュートリアルを完了してください。

パーシステンス図の計算¶

解析対象のデータは pointcloud.txt というファイルです。これを用いて、パーシステンス図を計算してみましょう。

In [2]:
# ライブラリの読み込み
import homcloud.interface as hc  # HomCloudのインターフェス
import numpy as np
import matplotlib.pyplot as plt

import pyvista as pv  # PyVista (3D可視化)
In [3]:
pv.set_jupyter_backend('trame')
Out[3]:
<MagicMock name='mock()' id='130311050306736'>

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

In [4]:
pointcloud = np.loadtxt("pointcloud.txt")

点群の情報を表示します。

In [5]:
pointcloud.shape
Out[5]:
(1000, 3)
In [6]:
np.min(pointcloud, axis=0), np.max(pointcloud, axis=0)
Out[6]:
(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]$という立方体での周期境界で計算します．

In [7]:
hc.PDList.from_alpha_filtration(pointcloud, 
                                save_to="pointcloud-periodic.pdgm",
                                periodicity=[(0, 1), (0, 1), (0, 1)],
                                save_boundary_map=True)
Out[7]:
PDList(path=pointcloud-periodic.pdgm)

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

注意¶

  • ポイントクラウドデータが小さすぎる場合には周期境界がうまく動作しません。その場合には周期的なコピーをX，Y，Z方向にそれぞれ数倍複製すると良いでしょう．
  • 今のところ周期境界の単位胞としては直方体しかサポートしていません。一般の平行六面体は利用できません。
  • 周期境界でない場合と同様に，点群が以下の条件を満たす場合には計算がうまくいきません。例えば結晶データの場合はこの問題が発生します。その場合各点の座標に微小なノイズを加えてこの条件を満たさないようにしてください。
    • 3点が同一直線上にある
    • 4点が同一平面上にある
    • 5点が同一球面上にある

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

パーシステンス図をプロットします。

ファイルを読み込み，1次元パーシステンス図を取り出し，プロットします。 周期境界条件なしの場合とまったく同じです。

In [8]:
pdlist = hc.PDList("pointcloud-periodic.pdgm")
In [9]:
pd1 = pdlist.dth_diagram(1)
In [10]:
pd1.histogram().plot(colorbar={"type": "log"})
No description has been provided for this image

0 から 0.01 の部分を拡大して表示します．

In [11]:
pd1.histogram((0, 0.01)).plot(colorbar={"type": "log"})
No description has been provided for this image

数値データの読み出し¶

これも周期境界条件なしの場合とまったく同じです。

In [12]:
pd1.births
Out[12]:
array([0.00103596, 0.00100508, 0.00116762, ..., 0.01369992, 0.01460543,
       0.01508903], shape=(3956,))
In [13]:
pd1.deaths
Out[13]:
array([0.00106352, 0.00108161, 0.00119129, ..., 0.01416671, 0.0148616 ,
       0.0153579 ], shape=(3956,))

逆解析について¶

逆解析についても周期境界なしの場合と良く似ていますが，以下の点で異なります。

  • optimal / stable volume は利用可能だが，PH trees は利用できない
  • 可視化の時に工夫をしないと境界部分で不自然に見える

本チュートリアルでは (0.004, 0.007) に一番近いペアに逆解析を適用してみましょう。

In [14]:
pair = pd1.nearest_pair_to(0.0022, 0.0057)
In [15]:
pair
Out[15]:
Pair(0.0021688455616267486, 0.005663084335456451)

optimal volume を計算します。

In [16]:
optimal_volume = pair.optimal_volume()

このoptimal volumeをPyVistaを用いて3次元可視化します。

In [17]:
pl = pv.Plotter()
pl.add_mesh(optimal_volume.to_pyvista_boundary_mesh())
pl.show()
Out[17]:
<MagicMock name='mock()' id='130311046818960'>

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

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

In [18]:
pl = pv.Plotter()
pl.add_mesh(optimal_volume.to_pyvista_boundary_mesh(adjust_periodic_boundary=(0.3, None)))  # 周期境界を考慮した描画をする
pl.add_mesh(pv.Box(bounds=(0, 1, 0, 1, 0, 1), level=2), style='wireframe', color="red")  # 単位胞を描画する
pl.show()
Out[18]:
<MagicMock name='mock()' id='130311046818960'>

ループが周期境界を通り抜けている様子がわかります。 さらに見やすくするため，単位胞の外側でループを複製して表示しましょう。

adjust_periodic_boundary=(0.3, 0.3)とします．2つめの0.3は周期境界の「端」から±0.3の割合の所にある線や頂点を複製します．

In [19]:
pl = pv.Plotter()
pl.add_mesh(optimal_volume.to_pyvista_boundary_mesh(adjust_periodic_boundary=(0.3, 0.3)))
pl.add_mesh(pv.Box(bounds=(0, 1, 0, 1, 0, 1), level=2), style='wireframe', color="red")  # 単位胞を描画する
pl.show()
Out[19]:
<MagicMock name='mock()' id='130311046818960'>

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

In [20]:
optimal_volume.boundary_points()
Out[20]:
[[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の範囲内で出力しています． この逆解析結果を使ってさらなる解析をしたい場合には適当に座標を調整する必要があるでしょう．

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