示例09 图像 Lower-Star 过滤

"""ex09_lower_star_image.py — 第 7 章 图像 Lower-Star 过滤

对应教程：tutorials/07-lower-star.md。

演示：

• lower_star_img 对灰度图计算 H0 持久图
• 合成「黑底亮环」图，观察 H1 出现（用负图做下星等效）
• 加噪声前后持久性变化

运行：
python ex09_lower_star_image.py
"""

import numpy as np
from ripser import lower_star_img
from persim import plot_diagrams
import matplotlib.pyplot as plt

def make_ring_image(size=64, center=(32, 32), radius=18, width=3):
"""合成黑底亮环图：背景 0，环上像素高值。"""
yy, xx = np.mgrid[0:size, 0:size]
dist = np.sqrt((xx - center[0])**2 + (yy - center[1])**2)
img = np.zeros((size, size), dtype=float)
img[np.abs(dist - radius) <= width / 2] = 1.0
return img

def main():
img = make_ring_image(size=64)
print(f"图像 shape={img.shape}, 像素值范围 [{img.min()}, {img.max()}]")

# ── 1. lower_star_img 默认算 H0（下水平集）──
dgm_h0 = lower_star_img(img)
print(f"\n下水平集 H0: {len(dgm_h0)} 点")

# ── 2. 对负图做下星 → 等效「上水平集」，看亮区结构 ──
dgm_super = lower_star_img(-img)
pers_super = dgm_super[:, 1] - dgm_super[:, 0]
pers_super = pers_super[np.isfinite(pers_super)]
print(f"上水平集 H0: {len(dgm_super)} 点, "
      f"max persistence={pers_super.max() if len(pers_super) else 0:.4f}")
print("  （亮环会形成一个显著的高持久 H0/H1 类特征）")

# ── 3. 加噪声对比 ──
rng = np.random.default_rng(0)
img_noisy = img + rng.uniform(-0.2, 0.2, img.shape)
dgm_noisy = lower_star_img(-img_noisy)
pers_noisy = dgm_noisy[:, 1] - dgm_noisy[:, 0]
pers_noisy = pers_noisy[np.isfinite(pers_noisy)]
print(f"加噪后 max persistence="
      f"{pers_noisy.max() if len(pers_noisy) else 0:.4f}")

# ── 可视化 ──
fig, axes = plt.subplots(2, 2, figsize=(11, 8))
axes[0, 0].imshow(img, cmap='gray')
axes[0, 0].set_title('合成亮环图')

axes[0, 1].imshow(img_noisy, cmap='gray')
axes[0, 1].set_title('加噪图')

plot_diagrams([dgm_super], labels=['上水平集 H0'],
              ax=axes[1, 0], show=False)
axes[1, 0].set_title('原始图 上水平集持久图')

plot_diagrams([dgm_noisy], labels=['上水平集 H0'],
              ax=axes[1, 1], show=False)
axes[1, 1].set_title('加噪图 上水平集持久图')

plt.tight_layout()
plt.show()
​

if name == "main":
main()