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 до оттенков серого см.: /questions/31903217/konvertirovat-jpg-v-csv-v-ottenkah-serogo-ispolzuya-r/31903233#31903233

# 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. Я думаю, что важно отметить, что одно из больших отличий заключается в том, что "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://stackru.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)