点群データの一部を相対ホモロジーでマスクするチュートリアル

このチュートリアルでは、点群データの一部を相対ホモロジーでマスクする解析について解説します。

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

ライブラリのインポート

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

%load_ext autoreload
%autoreload 2

import homcloud.interface as hc
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches
from scipy.spatial import Voronoi, voronoi_plot_2d

点群データの準備

ここでは、$[-1, 1]^2$上のランダムな点群の中央の曲がった帯状部分をマスクします。 ここでは2種類のランダムな2次元点群を用います。3次元でも計算は可能です。

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

pc = rng.uniform(-1, 1, size=(200, 2))
mask_cond = np.abs(0.5 * np.sin(pc[:, 0] * np.pi) - pc[:, 1]) < 0.35

P = pc[~mask_cond, :]
Q = pc[mask_cond, :]

# 点群のプロット
fig, ax = plt.subplots()
ax.scatter(P[:, 0], P[:, 1], label="P")
ax.scatter(Q[:, 0], Q[:, 1], label="Q")
ax.set_aspect("equal")
ax.legend()

<matplotlib.legend.Legend at 0x7cc0687d1010>

マスクされたPD

$(H_*([\cup_{p \in P} B_r(p)]\cup M, M)_{r \geq 0}$という相対ホモロジーのPHを計算します。

ただし、$\{V_x\}_{x \in P \cup Q}$は$P \cup Q$のボロノイ図で、$M = \cup_{q \in Q} V_q$です。

まずは領域$M$を可視化しましょう。

# ダミー点（Bounding Box）の追加
# データの範囲の外側に4つのダミー点を配置する
# これらのダミー点のおかげで、Qに対応するボロノイセルは有限領域になる。
min_dist = -10  # 十分に外側
max_dist = 10
dummy_points = np.array([
    [min_dist, min_dist],
    [min_dist, max_dist],
    [max_dist, min_dist],
    [max_dist, max_dist]
])

# ボロノイ図を計算
vor = Voronoi(np.vstack([P, Q, dummy_points]))

# 可視化の設定
fig, ax = plt.subplots(figsize=(8, 8))

# ボロノイ図の基本プロット
# show_vertices=False で頂点を非表示、line_colorsで境界線の色を指定
voronoi_plot_2d(vor, ax=ax, show_points=False, show_vertices=False, line_colors='black', line_width=1, line_alpha=0.5)

q_start_idx = len(P)
q_end_idx = len(P) + len(Q)

for i in range(q_start_idx, q_end_idx):
    # 各点に対応する領域（region）のインデックスを取得
    region_index = vor.point_region[i]
    region = vor.regions[region_index]

    polygon = [vor.vertices[idx] for idx in region]
    ax.fill(*zip(*polygon), color='salmon', alpha=0.5)

# 点群のプロット
ax.scatter(P[:, 0], P[:, 1], c='blue', s=20)
ax.scatter(Q[:, 0], Q[:, 1], c='red', s=20)

plt.gca().set_aspect('equal', adjustable='box')
plt.xlim(-1.1, 1.1)
plt.ylim(-1.1, 1.1)
plt.show()

さらに$P$の上に半径0.1の円盤を置いた図も描画しましょう。

fig, ax = plt.subplots(figsize=(8, 8))


for p in P:
    ax.add_artist(matplotlib.patches.Circle(p, 0.1, color="g"))

voronoi_plot_2d(vor, ax=ax, show_points=False, show_vertices=False, line_colors='black', line_width=1, line_alpha=0.5)

for i in range(q_start_idx, q_end_idx):
    region_index = vor.point_region[i]
    region = vor.regions[region_index]
    polygon = [vor.vertices[idx] for idx in region]
    ax.fill(*zip(*polygon), color='salmon', alpha=0.5)

ax.scatter(P[:, 0], P[:, 1], c='blue', s=20)
ax.scatter(Q[:, 0], Q[:, 1], c='red', s=20)

plt.gca().set_aspect('equal', adjustable='box')
plt.xlim(-1.1, 1.1)
plt.ylim(-1.1, 1.1)
plt.show()

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

hc.PDList.from_alpha_voronoi_relative_mask_filtration(P, Q, save_boundary_map=True, save_to="relative_mask.pdgm")
pdlist = hc.PDList("relative_mask.pdgm")

相対PDの可視化

可視化は通常のPDと同じです。1次元PDのみ可視化します。

pdlist[1].histogram((-0.01, 0.05), 64).plot(colorbar={"type": "log"})

マスクなしのPDも計算しましょう。マスクありと一緒に横に並べます。

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

pdlist_normal = hc.PDList.from_alpha_filtration(pc)

fig, axes = plt.subplots(1, 2, figsize=(10, 6))

pdlist[1].histogram((-0.01, 0.05), 64).plot(colorbar={"type": "log"}, ax=axes[0], title="masked")
pdlist_normal[1].histogram((-0.01, 0.05), 64).plot(colorbar={"type": "log"}, ax=axes[1], title="P and Q")
plt.tight_layout()

かなり違いがあります。

ここで、「マスク」しているということは、$P\cup Q$のPDから$Q$のPDを取り除けば似たような効果になるのではと考えられます。 そこで、そのような計算をしてみましょう。

pdlist_Q = hc.PDList.from_alpha_filtration(Q)
birthdeath_Q = set((p.birth_time(), p.death_time()) for p in pdlist_Q[1].pairs())
birthdeath_PQ_minus_Q = [
    (p.birth_time(), p.death_time()) 
    for p in pdlist_normal[1].pairs() 
    if (p.birth_time(), p.death_time()) not in birthdeath_Q
]

fig, axes = plt.subplots(1, 3, figsize=(12, 6))

pdlist[1].histogram((-0.01, 0.05), 64).plot(colorbar={"type": "log"}, ax=axes[0], title="masked")
pdlist_normal[1].histogram((-0.01, 0.05), 64).plot(colorbar={"type": "log"}, ax=axes[2], title="P and Q")


births_PQ_minus_Q = np.array([b for (b, d) in birthdeath_PQ_minus_Q])
deaths_PQ_minus_Q = np.array([d for (b, d) in birthdeath_PQ_minus_Q])
pd1_PQ_minus_Q = hc.PD.from_birth_death(1, births_PQ_minus_Q, deaths_PQ_minus_Q)
pd1_PQ_minus_Q.histogram((-0.01, 0.05), 64).plot(colorbar={"type": "log"}, ax=axes[1], title="PQ\\Q")

plt.tight_layout()

ちょっと近付いたかな、という感じです。

逆解析の適用

HomCloudは相対PHに関する逆解析機能(optimal volume と stable volume)を持っています。ただし、この機能は試験的(experimental)な機能です。

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

pair = pdlist[1].nearest_pair_to(0.005, 0.042)
pair

Pair(0.005016878293821473, 0.041934428640154967)

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

pair.lifetime() / 10

np.float64(0.0036917550346333497)

stable_volume = pair.stable_volume(pair.lifetime() / 10)

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

stable_volume.boundary_points()

[[-0.2267978771111614, -0.8996228507904607],
 [0.08308522618766645, -0.4912833027176391],
 [-0.2975939732909938, -0.9886573640621745],
 [0.19665423119327063, -0.30150668295485694],
 [-0.5326449376597135, -0.7834406313459001],
 [0.035342249990133245, -0.16110218719718516],
 [0.22014223877168804, -0.17109390853228912],
 [0.014750923198010968, -0.5684214133107457],
 [-0.5040852062409231, -0.9197913306766388],
 [-0.10052026869704278, -0.834475290401949],
 [-0.4016271392371318, -0.9622302023465905],
 [-0.08835728636230478, -0.6608371749239166],
 [0.17366582406959963, -0.3939103699823414],
 [0.20232082351324787, -0.26294381064273],
 [-0.1368139569819018, -0.8434039346331299]]

stable_volume.boundary()

[[[0.08308522618766645, -0.4912833027176391],
  [0.17366582406959963, -0.3939103699823414]],
 [[-0.2975939732909938, -0.9886573640621745],
  [-0.4016271392371318, -0.9622302023465905]],
 [[0.19665423119327063, -0.30150668295485694],
  [0.17366582406959963, -0.3939103699823414]],
 [[0.19665423119327063, -0.30150668295485694],
  [0.20232082351324787, -0.26294381064273]],
 [[-0.5040852062409231, -0.9197913306766388],
  [-0.5326449376597135, -0.7834406313459001]],
 [[0.08308522618766645, -0.4912833027176391],
  [0.014750923198010968, -0.5684214133107457]],
 [[-0.1368139569819018, -0.8434039346331299],
  [-0.10052026869704278, -0.834475290401949]],
 [[-0.08835728636230478, -0.6608371749239166],
  [0.014750923198010968, -0.5684214133107457]],
 [[0.22014223877168804, -0.17109390853228912],
  [0.035342249990133245, -0.16110218719718516]],
 [[-0.2975939732909938, -0.9886573640621745],
  [-0.2267978771111614, -0.8996228507904607]],
 [[-0.10052026869704278, -0.834475290401949],
  [-0.08835728636230478, -0.6608371749239166]],
 [[-0.4016271392371318, -0.9622302023465905],
  [-0.5040852062409231, -0.9197913306766388]],
 [[-0.1368139569819018, -0.8434039346331299],
  [-0.2267978771111614, -0.8996228507904607]],
 [[0.20232082351324787, -0.26294381064273],
  [0.22014223877168804, -0.17109390853228912]]]

stable_volume.cells()

[[[0.20232082351324787, -0.26294381064273],
  [0.22014223877168804, -0.17109390853228912],
  [0.035342249990133245, -0.16110218719718516]],
 [[-0.1368139569819018, -0.8434039346331299],
  [-0.10052026869704278, -0.834475290401949],
  [-0.08835728636230478, -0.6608371749239166]],
 [[-0.2975939732909938, -0.9886573640621745],
  [-0.4016271392371318, -0.9622302023465905],
  [-0.2267978771111614, -0.8996228507904607]],
 [[0.19665423119327063, -0.30150668295485694],
  [0.20232082351324787, -0.26294381064273],
  [0.035342249990133245, -0.16110218719718516]],
 [[-0.4016271392371318, -0.9622302023465905],
  [-0.5040852062409231, -0.9197913306766388],
  [-0.5326449376597135, -0.7834406313459001]],
 [[-0.4016271392371318, -0.9622302023465905],
  [-0.5326449376597135, -0.7834406313459001],
  [-0.42677828106990856, -0.6925087544206239]],
 [[-0.4016271392371318, -0.9622302023465905],
  [-0.2267978771111614, -0.8996228507904607],
  [-0.42677828106990856, -0.6925087544206239]],
 [[0.19665423119327063, -0.30150668295485694],
  [0.17366582406959963, -0.3939103699823414],
  [0.035342249990133245, -0.16110218719718516]],
 [[-0.1368139569819018, -0.8434039346331299],
  [-0.2267978771111614, -0.8996228507904607],
  [-0.42677828106990856, -0.6925087544206239]],
 [[-0.08835728636230478, -0.6608371749239166],
  [-0.3239018854854909, -0.5499753882518639],
  [-0.42677828106990856, -0.6925087544206239]],
 [[0.08308522618766645, -0.4912833027176391],
  [0.17366582406959963, -0.3939103699823414],
  [0.035342249990133245, -0.16110218719718516]],
 [[-0.1368139569819018, -0.8434039346331299],
  [-0.08835728636230478, -0.6608371749239166],
  [-0.42677828106990856, -0.6925087544206239]],
 [[-0.08835728636230478, -0.6608371749239166],
  [-0.3239018854854909, -0.5499753882518639],
  [-0.30375450566389306, -0.38705099291249345]],
 [[-0.08835728636230478, -0.6608371749239166],
  [0.014750923198010968, -0.5684214133107457],
  [-0.30375450566389306, -0.38705099291249345]],
 [[0.08308522618766645, -0.4912833027176391],
  [0.035342249990133245, -0.16110218719718516],
  [-0.1341482191418022, -0.19326774627625887]],
 [[0.08308522618766645, -0.4912833027176391],
  [-0.1341482191418022, -0.19326774627625887],
  [-0.25474721045536586, -0.26484394966771174]],
 [[0.08308522618766645, -0.4912833027176391],
  [0.014750923198010968, -0.5684214133107457],
  [-0.30375450566389306, -0.38705099291249345]],
 [[0.08308522618766645, -0.4912833027176391],
  [-0.25474721045536586, -0.26484394966771174],
  [-0.30375450566389306, -0.38705099291249345]]]

可視化しましょう

from matplotlib.patches import Polygon
from matplotlib.lines import Line2D

fig, ax = plt.subplots(figsize=(8, 8))

# ボロノイ図のプロット
voronoi_plot_2d(vor, ax=ax, show_points=False, show_vertices=False, line_colors='black', line_width=1, line_alpha=0.2)

# 点群のプロット
ax.scatter(P[:, 0], P[:, 1])
ax.scatter(Q[:, 0], Q[:, 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")

ax.set_xlim(-1.1, 1.1)
ax.set_ylim(-1.1, 1.1)

(-1.1, 1.1)

マスクされた領域が作る境界と点群とで構成されるループになっています。 理論上、マスクされた領域が作る境界はある意味でPDに対応するループの構成要素にはならないのでその部分は黒く描画されません。

3次元の例

3次元の場合もやってみます。

rng = np.random.Generator(np.random.PCG64(seed=1784594745233))
P = rng.uniform(-1, 1, size=(200, 3))
Q = np.vstack([
    np.hstack([rng.uniform(-0.6, 0.6, size=(100, 2)), rng.uniform(1.01, 1.30, size=(100, 1))]),
    np.hstack([rng.uniform(-0.6, 0.6, size=(100, 2)), rng.uniform(-1.30, -1.01, size=(100, 1))]),
])

fig = plt.figure()
ax = fig.add_subplot(projection='3d')

ax.scatter(P[:, 0], P[:, 1], P[:, 2])
ax.scatter(Q[:, 0], Q[:, 1], Q[:, 2])

ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
plt.show()

plt.scatter(P[:, 0], P[:, 2])
plt.scatter(Q[:, 0], Q[:, 2])


plt.xlabel("X")
plt.ylabel("Z")
plt.gca().set_aspect('equal', adjustable='box')

このように上下にマスクを指定します。

以下PDを計算してプロットします。

pdlist = hc.PDList.from_alpha_voronoi_relative_mask_filtration(P, Q)

pdlist.dth_diagram(1).histogram((-0.01, 0.15), 64).plot(colorbar={"type": "log"}, plot_ess=True)

pdlist.dth_diagram(1).essential_births

array([0.03220152])

すると、消滅時刻が+∞になる生成消滅対が現われます。 マスク領域がこのように2つ以上に分かれている場合は、1次元以上のPDに消滅時刻が+∞になる生成消滅対が現われるようになります。 どの次元にどのように現われるかは、マスク領域の形状に依存します。