R：查找二维点云的 2 个点之间的最短测地线路径

Question

我在 Vincent Zoonekynd 编写的两个函数的帮助下创建了下图（您可以找到它们here）（在帖子末尾找到我的代码）。

为了能够解释Isometric Feature Mapping 使用的邻域图和参数“k”是什么。 “k”指定每个点直接连接到多少个点。它们的距离只是彼此之间的欧几里得距离。任何点与其 (k + 1) 最近点（或任何更远的点）之间的距离称为“测地线”，是到达那里所需的所有边缘长度的最小总和。这有时比欧几里得距离长得多。我图中的 A 点和 B 点就是这种情况。

现在我想添加一条黑线，显示从 A 点到 B 点的测地线距离。我知道命令 segments()，这可能是添加线的最佳选择，而且我知道要找到的一种算法最短路径（测地线距离）是 Dijkstra 算法，它在包igraph 中实现。但是，我既不能让igraph 解释我的图表，也不能自己找出需要传递的点（顶点）（及其坐标）。

顺便说一下，如果 k = 18，即如果每个点都直接连接到最近的 18 个点，那么 A 和 B 之间的测地线距离将只是欧几里得距离。

---

isomap.incidence.matrix <- function (d, eps=NA, k=NA) {
  stopifnot(xor( is.na(eps), is.na(k) ))
  d <- as.matrix(d)
  if(!is.na(eps)) {
    im <- d <= eps
  } else {
    im <- apply(d,1,rank) <= k+1
    diag(im) <- FALSE
  }
  im | t(im)
}

plot.graph <- function (im,x,y=NULL, ...) {
  if(is.null(y)) {
    y <- x[,2]
    x <- x[,1]
  }
  plot(x,y, ...)
  k <- which(  as.vector(im)  )
  i <- as.vector(col(im))[ k ]
  j <- as.vector(row(im))[ k ]
  segments( x[i], y[i], x[j], y[j], col = "grey")
}

z <- seq(1.1,3.7,length=140)*pi

set.seed(4)
zz <- rnorm(1:length(z))+z*sin(z)
zz <- cbind(zz,z*cos(z)*seq(3,1,length=length(z)))

dist.grafik <- dist(zz)

pca.grafik <- princomp(zz)

x11(8, 8)
par(mar=c(0,0,0,0))
plot.graph(isomap.incidence.matrix(dist.grafik, k=3), pca.grafik$scores[,1], pca.grafik$scores[,2],
           xaxt = "n", yaxt = "n", xlab = "", ylab = "", cex = 1.3)
legend("topright", inset = 0.02, legend = "k = 3", col = "grey", lty = 1, cex = 1.3) 
segments(x0 = -8.57, y0 = -1.11, x1 = -10.83, y1 = -5.6, col = "black", lwd = 2, lty = "dashed")
text(x = -8.2, y = -1.4, labels = "A", font = 2, cex = 1.2)
text(x = -11, y = -5.1, labels = "B", font = 2, cex = 1.2)

Answer 1

以下代码可能会对您有所帮助，它使用您的数据创建一个 igraph 对象，在您的情况下，权重是节点之间的欧几里德距离。 然后找到sp$vpath[[1]] 返回的加权最短路径。在以下示例中，它是节点 5 和 66 之间的最短路径。 我用mattu的解决方案编辑了代码

isomap.incidence.matrix <- function (d, eps=NA, k=NA) {
  stopifnot(xor( is.na(eps), is.na(k) ))
  d <- as.matrix(d)
  if(!is.na(eps)) {
    im <- d <= eps
  } else {
    im <- apply(d,1,rank) <= k+1
    diag(im) <- FALSE
  }
  im | t(im)
}

plot.graph <- function (im,x,y=NULL, ...) {
  if(is.null(y)) {
    y <- x[,2]
    x <- x[,1]
  }
  plot(x,y, ...)
  k <- which(  as.vector(im)  )
  i <- as.vector(col(im))[ k ]
  j <- as.vector(row(im))[ k ]
  segments( x[i], y[i], x[j], y[j], col = "grey")
}

z <- seq(1.1,3.7,length=100)*pi

set.seed(4)
zz <- rnorm(1:length(z))+z*sin(z)
zz <- cbind(zz,z*cos(z)*seq(3,1,length=length(z)))

dist.grafik <- as.matrix(dist(zz))
pca.grafik <- princomp(zz)

isomap.resul <-  function (d, eps=NA, k=NA) {
  a <- isomap.incidence.matrix(d, eps, k)
  b <- dist.grafik
  res <- a * b
  return(res)
}

a <- graph_from_adjacency_matrix(isomap.resul(dist.grafik, k=3), 
                                 mode = c("undirected"), weight = TRUE)
sp <- shortest_paths(a, 5, to = 66, mode = c("out", "all", "in"),
                     weights = NULL, output = c("vpath", "epath", "both"),
                     predecessors = FALSE, inbound.edges = FALSE)

path <- sp$vpath[[1]] 

x11(8, 8)
par(mar=c(0,0,0,0))
plot.graph(isomap.incidence.matrix(dist.grafik, k=3), pca.grafik$scores[,1], pca.grafik$scores[,2],
           xaxt = "n", yaxt = "n", xlab = "", ylab = "", cex = 1.3)
legend("topright", inset = 0.02, legend = "k = 3", col = "grey", lty = 1, cex = 1.3) 
segments(x0 = -8.57, y0 = -1.11, x1 = -10.83, y1 = -5.6, col = "black", lwd = 2, lty = "dashed")
text(x = -8.2, y = -1.4, labels = "A", font = 2, cex = 1.2)
text(x = -11, y = -5.1, labels = "B", font = 2, cex = 1.2)

for(i in 2:length(path)){
  aa <- pca.grafik$scores[path[i-1], 1]
  bb <- pca.grafik$scores[path[i-1], 2]
  cc <- pca.grafik$scores[path[i], 1]
  dd <- pca.grafik$scores[path[i], 2]
  segments(aa, bb, cc , dd, lwd = 2)
}

要运行这个脚本，你显然需要包igraph。

对我来说，根据测地距离，这似乎是最短的路径。

希望对你有帮助。