R中持久同伦：识别生成拓扑特征的点

我正在用R语言的TDA包，通过gridDiag()函数完成持久同伦分析，重点关注1维单纯复形（环路）相关的点。设置location = TRUE参数后，能通过cycleLocation把单纯复形绘制回点云，示例代码如下：

# 生成数据
set.seed(2)
x = runif(60, min=0, max=100)
y = runif(60, min=0, max=100)
coords <- cbind(x,y)
plot(coords)

# 计算持久同伦，设置location = TRUE
library(TDA)
Xlim=c(min(coords[,1]), max(coords[,1]))
Ylim=c(min(coords[,2]), max(coords[,2]))
by=1
lim = cbind(Xlim, Ylim)
Diag <- gridDiag(coords, distFct, lim = lim, by = by, sublevel = TRUE,
                 library = "Dionysus", location = TRUE, printProgress = TRUE)

# 绘图
par(mfrow = c(1, 3))
plot(coords, cex = 0.5, pch = 19)
title(main = "数据")
threshold = 1 # 用于绘制拓扑特征的持久值阈值
plot(Diag[["diagram"]], band = 2*threshold)
title(main = "距离函数图")
one <- which(Diag[["diagram"]][, 1] == 1 & sqrt(0.5*(Diag[["diagram"]][, "Death"]-Diag[["diagram"]][, "Birth"]))>threshold)
plot(coords, col = 2, main = "网格点的代表性环路")
for (i in seq(along = one)) {
  points(Diag[["birthLocation"]][one[i], , drop = FALSE], pch = 15, cex = 3,
         col = i)
  points(Diag[["deathLocation"]][one[i], , drop = FALSE], pch = 17, cex = 3,
         col = i)
  for (j in seq_len(dim(Diag[["cycleLocation"]][[one[i]]])[1])) {
    lines(Diag[["cycleLocation"]][[one[i]]][j, , ], pch = 19, cex = 1, col = i)
  }
}

示例代码生成的图包含三个子图：原始数据点云、距离函数持久同伦图、网格点的代表性环路可视化。

目前得到的是增长半径球之间的空空间，我想找简便方法获取启动环路的原始点坐标——也就是识别环路诞生时，因球重叠生成环路的那些原始点。

之前有类似问题被提出，但解决方案用的聚类算法只适用于特定数据集。我的数据集（比如示例里的随机分布点）分布不清晰，希望直接从已完成的计算结果中提取：理想形式是一个列表，每个子列表对应一个阈值化单纯复形，包含coords中生成该复形的顶点索引向量。

解决方案 #

1. 核心思路 #

gridDiag结合Dionysus计算时，Diag对象已存储了环路的诞生位置和诞生半径（以平方距离形式）。我们只需要找到原始数据中，与环路诞生位置的距离小于等于诞生半径的点——这些点就是球重叠形成环路的关键原始点，无需额外聚类。

2. 提取代码实现 #

在你的分析代码后添加以下内容：

# 定义函数：根据环路信息提取诞生时的原始点索引
get_cycle_birth_vertices <- function(diag_obj, coords, cycle_idx) {
  # 获取当前环路的诞生位置与平方诞生值
  birth_loc <- diag_obj[["birthLocation"]][cycle_idx, ]
  birth_sq <- diag_obj[["diagram"]][cycle_idx, "Birth"]
  # 转换为实际半径（gridDiag使用平方欧氏距离，对应子水平集的距离为sqrt(0.5*平方距离)）
  birth_radius <- sqrt(0.5 * birth_sq)
  
  # 计算每个原始点到诞生位置的欧氏距离
  dists <- sqrt(rowSums((coords - birth_loc)^2))
  # 返回距离小于等于诞生半径的点的索引
  which(dists <= birth_radius)
}

# 生成目标结果列表：每个子元素对应一个符合阈值的环路的原始点索引
cycle_vertex_list <- list()
for (idx in seq_along(one)) {
  cycle_vertex_list[[idx]] <- get_cycle_birth_vertices(Diag, coords, one[idx])
}

# 查看结果
print(cycle_vertex_list)

3. 结果验证 #

可以通过绘图确认提取的点是否对应环路核心区域：

# 绘制第一个环路的诞生原始点
plot(coords, cex = 0.5, pch = 19, main = "第一个环路诞生时的原始点")
points(coords[cycle_vertex_list[[1]], ], col = "red", pch = 19, cex = 1.2)

关键说明 #

• gridDiag默认采用平方欧氏距离计算子水平集，因此需要将Birth值转换为实际半径：sqrt(0.5 * Birth)，这是因为子水平集的距离函数对应d(x)^2 ≤ 2r²，转换后得到实际半径r。
• 该方法直接复用已完成的持久同伦计算结果，不需要额外训练或聚类，适配任意分布的数据集。
• 如果需要获取构成环路的1维单纯形（边）对应的顶点，可以进一步解析Dionysus的底层输出，但上述方法已满足“找到导致环路诞生的原始点”的核心需求。

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