R:查找二维点云两点间的最短测地线路径

I created the following graph with help of two functions written by Vincent Zoonekynd (you can find them here) (find my code at the end of the post).

我利用Vincent Zoonekynd编写的两个函数(您可以在这里找到它们)创建了下面的图表(请在本文末尾找到我的代码)。

In order to be able to explain what that neighbourhood graph and that parameter "k" is, which the Isometric Feature Mapping uses. "k" specifies how many points each point is directly connected to. Their distance is just the euclidian distance to each other. The distance between any point and its (k + 1)-nearest point (or any point farther away) is called "geodesic", and is the smallest sum of all the lengths of the edges needed to get there. This is sometimes much longer than the euclidian distance. This is the case for points A and B in my figure.

为了能够解释邻域图和参数k是什么，等距特征映射使用了它。“k”指定每个点直接连接到多少个点。它们之间的距离就是欧几里得距离。任何点与最近点(k + 1)之间的距离称为“测地线”，是到达那里所需的所有边的长度的最小和。这有时比欧几里得距离长得多。这是我图中A点和B点的情况。

Now I want to add a black line showing the geodesic distance from point A to point B. I know about the command segments(), which will probably be the best for adding the line, and I know that one algorithm to find the shortest path (geodesic distance) is Dijkstra's Algorithm and that it is implemented in the package igraph. However, I'm neither able to have igraph interpret my graph nor to find out the points (vertices) that need to be passed (and their coordinates) on my own.

现在我想添加一个黑线显示测地距离从a点到b点我知道命令段(),这可能是最好的添加,我知道找到最短路径的一个算法(测地距离)Dijkstra算法的算法,它在包igraph实现。然而，我既不能让igraph解释我的图，也不能找出需要自己传递的点(顶点)(以及它们的坐标)。

By the way, if k = 18, i.e. if every point is directly connected to the 18 nearest points, the geodesic distance between A and B will just be the euclidian distance.

顺便说一下，如果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)

1 个解决方案

#1

The following code may help you, it use your data to create an igraph object with weight that are in your case, the euclidean distances between nodes. Then you find the weighted shortest path which is returned by sp$vpath[[1]]. In the following example, it is the shortest path between nodes number 5 and 66. I edited the code with the solution to plot from mattu

下面的代码可以帮助您，它使用您的数据创建一个具有权值的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)
}

To run this script, you obviously need the package igraph.

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

To me it seems the shortest path according to the geodesic distance.

在我看来，它似乎是根据测地线距离的最短路径。

Hope it helps.

希望它可以帮助。