R - 通过光栅图像（迷宫）找到成本最低的路径？

Question

如何通过光栅图像数据找到非线性路径？例如，最低成本算法？起点和终点已知并给出为：

起点 = (0,0)
终点 = (12,-5)

例如，通过（灰度）光栅图像提取蜿蜒河流的大致路径。

# 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

在我写我的例子时，我偶然发现了一个使用 'gdistance' r 包的答案......希望其他人会发现这很有用。

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

我想确保我不会让自己太容易遇到问题......所以我制作了一个更嘈杂的图像版本。

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

更新：使用真实的栅格数据...
我想看看相同的工作流程是否适用于实际的真实光栅图像，而不仅仅是假数据，所以......

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

RGB转灰度，见：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")

在找到一个可行的过渡函数之前，我尝试了不同的过渡函数。这可能比它需要的更复杂，但它确实有效。您可以增加幂项（从 2 到 3、4、5、6...），它会继续工作。它没有找到去除幂项的正确解。

---

使用igraph 包的替代解决方案。

使用 'igraph' r 包找到了一组替代答案。我认为重要的是要注意，这里最大的区别之一是“igraph”支持 n 维图，而“gdistance”仅支持 2D 图。因此，例如，将这个答案扩展到 3D 是相对容易的。

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]

边缘权重计算来自：https://stackoverflow.com/a/27446127/2371031（谢谢！）

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

第二个（噪音较大的）示例需要不同的公式来计算边权重。

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)