R: Find shortest geodesic path between 2 points of a 2-D point cloud

Question

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).

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.

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.

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.

---

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

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

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.

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

Hope it helps.