点群データの相対パーシステントホモロジー解析¶

このチュートリアルでは、点群データの相対パーシステントホモロジー解析ついて解説します。

この機能はまだ試験的なので、インターフェースが変わったり機能自体がなくなる可能性があります。

ライブラリのインポート¶

まずは、解析に必要な Python ライブラリをインポートします。

In [2]:
import homcloud.interface as hc
import numpy as np
import matplotlib.pyplot as plt

点群データの準備¶

ここでは2種類のランダムな2次元点群を用います。3次元でも計算は可能です。 点群 X は円周上、Y は円盤内の一様分布です。

In [3]:
# 乱数アルゴリズムと乱数シードを固定することで、得られるランダムデータを再現可能にする
rng = np.random.Generator(np.random.PCG64(seed=8402))

t = np.linspace(0, 2 * np.pi, 20)
X = np.vstack([np.sin(t), np.cos(t)]).transpose() + rng.uniform(-0.1, 0.1, size=(20, 2))
Y = rng.uniform(-1, 1, size=(50, 2))
mask = np.linalg.norm(Y, axis=1) < 0.93
Y = Y[mask, :]
In [4]:
# 点群のプロット
fig, ax = plt.subplots()
ax.scatter(X[:, 0], X[:, 1])
ax.scatter(Y[:, 0], Y[:, 1])
ax.set_aspect("equal")
No description has been provided for this image

相対PDの計算¶

$B_r(X) = \cup_{x \in X} B_r(x)$という記法の下、$(H_*(B_r(X \cup Y), B_r(X)))_{r \geq 0}$という相対ホモロジー群のパーシステントホモロジーを考えます。

以下のセルで計算します。 後に逆解析を適用するために save_boundary_map=True という引数を追加しておきます。

In [5]:
hc.PDList.from_coupled_alpha_relative_filtration(X, Y, save_boundary_map=True, save_to="relative.pdgm")
pdlist = hc.PDList("relative.pdgm")

相対PDの可視化¶

可視化は通常のPDと同じです。

In [6]:
fig, axes = plt.subplots(1, 3, figsize=(12, 5))
ranges = [(-0.01, 0.05), (0, 0.1), (0, 1)]

for d in range(3):
    pdlist[d].histogram(ranges[d], 64).plot(ax=axes[d])

plt.tight_layout()
No description has been provided for this image

通常の2次元点群では2次元PDは生成消滅対を持ちませんが、相対ホモロジーを考えると2次元PDにも対が現われます。

$X \cup Y$からアルファ複体を使ってPDを計算して、比較してみましょう。

In [7]:
pdlist_normal = hc.PDList.from_alpha_filtration(np.vstack([X, Y]))
fig, axes = plt.subplots(1, 3, figsize=(12, 5))
ranges = [(-0.01, 0.05), (0, 0.1), (0, 1)]
for d in range(3):
    pdlist_normal[d].histogram(ranges[d], 64).plot(ax=axes[d])

plt.tight_layout()
No description has been provided for this image

$X\cup Y$のPDと$(X\cup Y, X)$に関する相対PDが異なっていることがわかります。

逆解析の適用¶

HomCloudは相対PHに関する逆解析機能(optimal volume と stable volume)を持っています。ただし、この機能は試験的(experimental)な機能で、出力結果に対する数学的な正当化は現状では存在しません。

$(0.02, 0.05)$付近のbirth-death pairを対象とします。

In [8]:
pair = pdlist[1].nearest_pair_to(0.02, 0.05)
pair
Out[8]:
Pair(0.021722880343497643, 0.052110525961362716)

ここでは stable volume を計算しましょう。epsilonはこのpairのlifetimeの1/10とします。

In [9]:
pair.lifetime() / 10
Out[9]:
np.float64(0.0030387645617865073)
In [10]:
stable_volume = pair.stable_volume(pair.lifetime() / 10)

情報の抽出法は通常の2D/3D点群の場合と同じです。

In [11]:
stable_volume.boundary_points()
Out[11]:
[[-0.1482899537807676, 0.10552173871366266],
 [-0.05582540313078943, 0.2910552409135023],
 [-0.7822657186432491, 0.5677363493557619],
 [-0.6223111815017949, 0.7310348965581595],
 [-0.5295934569092886, 0.3857465058637355],
 [-0.2236996789980048, 0.5619884121483119],
 [-0.0851144909483299, 0.4141882232605183],
 [-0.358987724057769, 0.1267950246294185],
 [-0.3120621450718857, 0.7752386089859276]]
In [12]:
stable_volume.boundary()
Out[12]:
[[[-0.7822657186432491, 0.5677363493557619],
  [-0.5295934569092886, 0.3857465058637355]],
 [[-0.6223111815017949, 0.7310348965581595],
  [-0.3120621450718857, 0.7752386089859276]],
 [[-0.358987724057769, 0.1267950246294185],
  [-0.5295934569092886, 0.3857465058637355]],
 [[-0.0851144909483299, 0.4141882232605183],
  [-0.2236996789980048, 0.5619884121483119]],
 [[-0.1482899537807676, 0.10552173871366266],
  [-0.05582540313078943, 0.2910552409135023]],
 [[-0.3120621450718857, 0.7752386089859276],
  [-0.2236996789980048, 0.5619884121483119]],
 [[-0.0851144909483299, 0.4141882232605183],
  [-0.05582540313078943, 0.2910552409135023]],
 [[-0.358987724057769, 0.1267950246294185],
  [-0.1482899537807676, 0.10552173871366266]]]
In [13]:
stable_volume.cells()
Out[13]:
[[[-0.7822657186432491, 0.5677363493557619],
  [-0.6223111815017949, 0.7310348965581595],
  [-0.5295934569092886, 0.3857465058637355]],
 [[-0.358987724057769, 0.1267950246294185],
  [-0.1482899537807676, 0.10552173871366266],
  [-0.05582540313078943, 0.2910552409135023]],
 [[-0.0851144909483299, 0.4141882232605183],
  [-0.358987724057769, 0.1267950246294185],
  [-0.05582540313078943, 0.2910552409135023]],
 [[-0.6223111815017949, 0.7310348965581595],
  [-0.3120621450718857, 0.7752386089859276],
  [-0.2236996789980048, 0.5619884121483119]],
 [[-0.6223111815017949, 0.7310348965581595],
  [-0.5295934569092886, 0.3857465058637355],
  [-0.2236996789980048, 0.5619884121483119]],
 [[-0.0851144909483299, 0.4141882232605183],
  [-0.358987724057769, 0.1267950246294185],
  [-0.2236996789980048, 0.5619884121483119]],
 [[-0.358987724057769, 0.1267950246294185],
  [-0.5295934569092886, 0.3857465058637355],
  [-0.2236996789980048, 0.5619884121483119]]]

可視化しましょう

In [14]:
from matplotlib.patches import Polygon
from matplotlib.lines import Line2D

# 点群のプロット
fig, ax = plt.subplots()
ax.scatter(X[:, 0], X[:, 1])
ax.scatter(Y[:, 0], Y[:, 1])

# 内部領域のプロット (薄い青色)
for triangle in stable_volume.cells():
    ax.add_artist(Polygon(triangle, color="b", alpha=0.2))

# 境界部のプロット (黒色)
for edge in stable_volume.boundary():
    ax.add_artist(Line2D([edge[0][0], edge[1][0]], [edge[0][1], edge[1][1]], color="k"))

# 縦横比を一致させる
ax.set_aspect("equal")
No description has been provided for this image