R - Finding least cost path through raster image (maze)?

Question

How can I find a non-linear path through raster image data? e.g., least cost algorithm? Starting and ending points are known and given as:

Start point = (0,0)
End point = (12,-5)

For example, extract the approximate path of a winding river through a (greyscale) raster image.

# fake up some noisy, but reproducible, "winding river" data
set.seed(123)
df <- data.frame(x=seq(0,12,by=.01), 
                 y=sapply(seq(0,12,by=.01), FUN = function(i) 10*sin(i)+rnorm(1)))

# convert to "pixels" of raster data
# assumption: image color is greyscale, only need one numeric value, v
img <- data.frame(table(round(df$y,0), round(df$x,1)))
names(img) <- c("y","x","v")
img$y <- as.numeric(as.character(img$y))
img$x <- as.numeric(as.character(img$x))


## take a look at the fake "winding river" raster image...
library(ggplot2)
ggplot(img) +
  geom_raster(aes(x=x,y=y,fill=v))

Answer 1

As I was writing up my example, I stumbled upon an answer using the 'gdistance' r package... hopefully others will find this useful.

library(gdistance)
library(sp)
library(ggplot2)

# convert to something rasterFromXYZ() understands
spdf <- SpatialPixelsDataFrame(points = img[c("x","y")], data = img["v"])

# use rasterFromXYZ to make a RasterLayer 
r <- rasterFromXYZ(spdf)

# make a transition layer, specifying a sensible function and the number of connection directions
tl <- transition(r, function(x) min(x), 8)
## mean(x), min(x), and max(x) produced similar results for me

# extract the shortest path as something we can plot
sPath <- shortestPath(tl, c(0,0), c(12,-5), output = "SpatialLines")

# conversion for ggplot
sldf <- fortify(SpatialLinesDataFrame(sPath, data = data.frame(ID = 1)))

# plot the original raster, truth (white), and the shortest path solution (green)   
ggplot(img) +
  geom_raster(aes(x=x,y=y,fill=v)) +
  stat_function(data=img, aes(x=x), fun = function(x) 10*sin(x), geom="line", color="white") +
  geom_path(data=sldf, aes(x=long,y=lat), color="green") 

I wanted to make sure that I wasn't just giving myself too easy of a problem... so I made a noisier version of the image.

img2 <- img
img2$v <- ifelse(img2$v==0, runif(sum(img2$v==0),3,8), img2$v)

spdf2 <- SpatialPixelsDataFrame(points = img2[c("x","y")], data = img2["v"])
r2 <- rasterFromXYZ(spdf2)

# for this noisier image, I needed a different transition function. 
# The one from the vignette worked well enough for this example.
tl2 <- transition(r2, function(x) 1/mean(x), 8)
sPath2 <- shortestPath(tl2, c(0,0), c(12,-5), output = "SpatialLines")
sldf2 <- fortify(SpatialLinesDataFrame(sPath2, data = data.frame(ID = 1)))

ggplot(img2) +
  geom_raster(aes(x=x,y=y,fill=v)) +
  stat_function(data=img2, aes(x=x), fun = function(x) 10*sin(x), geom="line", color="white") +
  geom_path(data=sldf2, aes(x=long,y=lat), color="green") 

UPDATE: using real raster data...
I wanted to see if the same workflow would work on an actual real-world raster image and not just fake data, so...

library(jpeg)
# grab some river image...
url <- "https://c8.alamy.com/comp/AMDPJ6/fiji-big-island-winding-river-aerial-AMDPJ6.jpg"
download.file(url, "river.jpg", mode = "wb")
jpg <- readJPEG("./river.jpg")
img3 <- melt(jpg, varnames = c("y","x","rgb"))
img3$rgb <- as.character(factor(img3$rgb, levels = c(1,2,3), labels=c("r","g","b")))
img3 <- dcast(img3, x + y ~ rgb)

# convert rgb to greyscale 
img3$v <- img3$r*.21 + img3$g*.72 + img3$b*.07

For rgb to greyscale, see: https://stackoverflow.com/a/27491947/2371031

# define some start/end point coordinates
pts_df <- data.frame(x = c(920, 500), 
                     y = c(880, 50))

# set a reference "grey" value as the mean of the start and end point "v"s
ref_val <- mean(c(subset(img3, x==pts_df[1,1] & y==pts_df[1,2])$v,
                  subset(img3, x==pts_df[2,1] & y==pts_df[2,2])$v))

spdf3 <- SpatialPixelsDataFrame(points = img3[c("x","y")], data = img3["v"])
r3 <- rasterFromXYZ(spdf3)

# transition layer defines "conductance" between two points
# x is the two point values, "v" = c(v1, v2) 
# 0 = no conductance, >>1 = good conductance, so
# make a transition function that encourages only small changes in v compared to the reference value. 
tl3 <- transition(r3, function(x) (1/max(abs((x/ref_val)-1))^2)-1, 8) 

sPath3 <- shortestPath(tl3, as.numeric(pts_df[1,]), as.numeric(pts_df[2,]), output = "SpatialLines")
sldf3 <- fortify(SpatialLinesDataFrame(sPath3, data = data.frame(ID = 1)))

# plot greyscale with points and path
ggplot(img3) +
  geom_raster(aes(x,y, fill=v)) +
  scale_fill_continuous(high="white", low="black") + 
  scale_y_reverse() +
  geom_point(data=pts_df, aes(x,y), color="red") + 
  geom_path(data=sldf3, aes(x=long,y=lat), color="green")

I played around with different transition functions before finding one that worked. This one is probably more complex than it needs to be, but it works. You can increase the power term (from 2 to 3,4,5,6...) and it continues to work. It did not find a correct solution with the power term removed.

---

Alternative solution using igraph package.

Found an alternative set of answers using 'igraph' r package. I think it is important to note that one of the big differences here is that 'igraph' supports n-dimensional graphs whereas 'gdistance' only supports 2D graphs. So, for example, extending this answer into 3D is relatively easy.

library(igraph)

# make a 2D lattice graph, with same dimensions as "img"
l <- make_lattice(dimvector = c(length(unique(img$y)), 
                                length(unique(img$x))), directed=F, circular=F)
summary(l)
# > IGRAPH ba0963d U--- 3267 6386 -- Lattice graph
# > + attr: name (g/c), dimvector (g/n), nei (g/n), mutual (g/l), circular (g/l)

# set vertex attributes
V(l)$x = img$x
V(l)$y = img$y
V(l)$v = img$v

# "color" is a known attribute that will be used by plot.igraph()
V(l)$color = grey.colors(length(unique(img$v)))[img$v+1]

# compute edge weights as a function of attributes of the two connected vertices
el <- get.edgelist(l)

# "weight" is a known edge attribute, and is used in shortest_path()
# I was confused about weights... lower weights are better, Inf weights will be avoided.
# also note from help: "if all weights are positive, then Dijkstra's algorithm is used."
E(l)$weight <- 1/(pmax(V(l)[el[, 1]]$v, V(l)[el[, 2]]$v))
E(l)$color = grey.colors(length(unique(E(l)$weight)))[E(l)$weight+1]

Edge weights calculation courtesy of: https://stackoverflow.com/a/27446127/2371031 (thanks!)

# find the start/end vertices
start = V(l)[V(l)$x == 0 & V(l)$y == 0]
end = V(l)[V(l)$x == 12 & V(l)$y == -5] 

# get the shortest path, returning "both" (vertices and edges)...
result <- shortest_paths(graph = l, from = start, to = end,  output = "both")

# color the edges that were part of the shortest path green
V(l)$color = ifelse(V(l) %in% result$vpath[[1]], "green", V(l)$color)
E(l)$color = ifelse(E(l) %in% result$epath[[1]], "green", E(l)$color)

# color the start and end vertices red
V(l)$color = ifelse(V(l) %in% c(start,end), "red", V(l)$color)

plot(l, vertex.shape = "square", vertex.size=2, vertex.frame.color=NA, vertex.label=NA, curved=F)

Second (noisier) example requires a different formula to compute edge weights.

img2 <- img
img2$v <- ifelse(img2$v==0, runif(sum(img2$v==0),3,8), img2$v)

l <- make_lattice(dimvector = c(length(unique(img2$y)),
                                length(unique(img2$x))), directed=F, circular=F)

# set vertex attributes
V(l)$x = img2$x
V(l)$y = img2$y
V(l)$v = img2$v
V(l)$color = grey.colors(length(unique(img2$v)))[factor(img2$v)]

# compute edge weights 
el <- get.edgelist(l)

# proper edge weight calculation is the key to a good solution...
E(l)$weight <- (pmin(V(l)[el[, 1]]$v, V(l)[el[, 2]]$v))
E(l)$color = grey.colors(length(unique(E(l)$weight)))[factor(E(l)$weight)]

start = V(l)[V(l)$x == 0 & V(l)$y == 0]
end = V(l)[V(l)$x == 12 & V(l)$y == -5]

# get the shortest path, returning "both" (vertices and edges)...
result <- shortest_paths(graph = l, from = start, to = end,  output = "both")

# color the edges that were part of the shortest path green
V(l)$color = ifelse(V(l) %in% result$vpath[[1]], "green", V(l)$color)
E(l)$color = ifelse(E(l) %in% result$epath[[1]], "green", E(l)$color)

# color the start and end vertices red
V(l)$color = ifelse(V(l) %in% c(start,end), "red", V(l)$color)

plot(l, vertex.shape = "square", vertex.size=2, vertex.frame.color=NA, vertex.label=NA, curved=F)