使用GUDHI从二进制掩码PNG生成正确persistence diagram的问题求助

问题描述 #

尝试从分割用的二进制掩码生成persistence diagram，但运行代码后图表始终仅显示两个红点，更换不同PNG或随机生成二进制掩码也一样，不过cofaces_of_persistence_pairs输出看似正常。代码如下：

import numpy as np
import gudhi as gd
import matplotlib.pyplot as plt
from PIL.Image import open as open

image = 1 - np.asarray(open(r'./samples/Image_01L_1stHO.png').convert('L'))
plt.imshow(image, cmap='gray')
plt.title("Original Image")
plt.show()

image_vector = image.flatten()

cubical_complex = gd.CubicalComplex(
    dimensions=[image.shape[0], image.shape[1]],
    top_dimensional_cells=image_vector
)

cubical_complex.compute_persistence(homology_coeff_field=2, min_persistence=0)
print(cubical_complex)

cofaces = cubical_complex.cofaces_of_persistence_pairs()

print("Cofaces of persistence pairs:", cofaces)

persistence_diagram = cubical_complex.persistence()
print("Persistence Diagram:", persistence_diagram)

gd.plot_persistence_diagram(persistence_diagram)
plt.show()

问题根源 #

问题出在立方复形的构造逻辑上：

• Gudhi的CubicalComplex默认按过滤值从小到大构建复形，0对应先出现的低维结构，1对应后出现的高维结构。
• 二进制掩码的0/1值会让所有背景被识别为同一个初始连通分支，所有前景被当作填充该分支的高维结构，最终只会产生0维的两个点（背景的出生、前景填充后连通分支的死亡），这就是图表仅显示两个红点的原因。

修复方案 #

方式1：调整二值图的过滤值逻辑 #

反转过滤值的含义，让前景区域的过滤值更小，保证前景的连通分支被优先识别，背景作为高维“空”结构：

# 替换原图像处理代码，根据掩码实际前景颜色调整判断条件
image = np.asarray(open(r'./samples/Image_01L_1stHO.png').convert('L'))
# 如果前景是黑色（值为0），用这行；如果是白色（值为255），改为image == 255
binary_mask = (image == 0).astype(int)
image_vector = binary_mask.flatten()

# 后续立方复形构造和计算代码不变
cubical_complex = gd.CubicalComplex(
    dimensions=[binary_mask.shape[0], binary_mask.shape[1]],
    top_dimensional_cells=image_vector
)

方式2：用距离变换捕捉拓扑特征（更适合分割掩码） #

对掩码做距离变换，用距离值作为过滤值，能有效捕捉前景的孔洞、分支等拓扑结构：

import scipy.ndimage as ndi

# 先得到二值掩码（前景为1，背景为0）
image = np.asarray(open(r'./samples/Image_01L_1stHO.png').convert('L'))
binary_mask = (image > 127).astype(int)  # 根据实际阈值调整

# 计算前景到背景的欧式距离变换
distance = ndi.distance_transform_edt(binary_mask)
# 反转距离值，让前景中心（距离大）的过滤值更小，构造复形时从中心向外扩展
filtered_image = np.max(distance) - distance

# 构造立方复形
image_vector = filtered_image.flatten()
cubical_complex = gd.CubicalComplex(
    dimensions=[filtered_image.shape[0], filtered_image.shape[1]],
    top_dimensional_cells=image_vector
)

额外优化建议 #

• 不要设置min_persistence=0，这会保留所有噪声拓扑点，建议设置小阈值过滤无效结构：

  cubical_complex.compute_persistence(homology_coeff_field=2, min_persistence=0.1)
  

• 绘图时开启图例，区分不同维度的拓扑特征：

  gd.plot_persistence_diagram(persistence_diagram, legend=True)
  plt.show()
  

内容的提问来源于stack exchange，提问作者Donal Huang