Hi,
Finally, I inspired myself from a code by Barry Rownlingson intially
dedicated to calculate the shortest path.
But in my cas, it consisted in calculating the longest one.
Here is the final code with an image as attached file. In the image, you
see the initial shape, the voronoi lines in rainbow colors and the
simplified median Line in bold black.
library(rgdal)
library(spatstat)
library(rgeos)
library(maptools)
library(igraph)
setwd("F:/METHODES/1305_SKELETON/")
pol <- readOGR("IN/polygone.shp", "polygone")
# FUNCTION
findMedianLine <- function(pol, eps, tol) {
l <- as(as(pol, "SpatialLines"), "psp")
pts <- pointsOnLines(l, eps=eps)
vor <- dirichlet(pts)
# LINES EXTRACTION
vor.l <- as(as(vor, "SpatialPolygons"), "SpatialLines")
vor.l <- gIntersection(vor.l, vor.l)@lines[[1]]@Lines
vor.split <- SpatialLines(lapply(1:length(vor.l), function(x)
Lines(list(vor.l[[x]]), ID=as.character(x))))
vor.l.in <- vor.l.split[gContains(pol, vor.split, byid=T), ]
# GRAPH CONSTRUCTION
edges = do.call(rbind,lapply(vor.l.in at lines
,function(ls){as.vector(t(ls at Lines[[1]]@coords))}))
lengths = sqrt((edges[,1]-edges[,3])^2+(edges[,2]-edges[,4])^2)
froms = paste(edges[,1],edges[,2])
tos = paste(edges[,3],edges[,4])
g = graph.edgelist(cbind(froms,tos),directed=FALSE)
E(g)$weight=lengths
# LONGEST PATH
sel <- get.diameter(g)
vor.sel <- vor.l.in at lines[sel-1]
spine <- SpatialLines(vor.sel)
spine.spl <- gSimplify(gLineMerge(spine), tol=tol)
return(spine.spl)
}
myMedianLine <- findMedianLine(pol, eps=1000, tol=1000)
2013/5/17 Mathieu Rajerison <mathieu.rajerison at gmail.com>
Hi,
I've written a code to skeletonize an area shape.
The example could be a watershed from which we'd like to compute the
"guessed" stream, or find the median axis of a shape.
I can find my median line but there are some adjacent segments connected
to it. The length criteria is not appropriate to deal with these adjacent
segments as some tiny ones could be part of to the median line.
Has anyone idea to deal with that?
I apoplogize in advance if it doesn't work in all cases (holes, rings,
etc...).
library(rgdal)
library(spatstat)
library(rgeos)
library(maptools)
setwd("F:/METHODES/1305_SKELETON/")
pol <- readOGR("IN/polygone.shp", "polygone")
# FUNCTION
findMedianLine <- function(pol, tol) {
l <- as(as(pol, "SpatialLines"), "psp")
pts <- pointsOnLines(l, eps=10000)
vor <- dirichlet(pts)
vor.l <- as(as(vor, "SpatialPolygons"), "SpatialLines")
vor.lines <- gIntersection(vor.l, vor.l)@lines[[1]]@Lines
vor.l.split <- SpatialLines(lapply(1:length(vor.lines), function(x)
Lines(list(vor.lines[[x]]), ID=as.character(x))))
medianLine <- gLineMerge(vor.l.split[gContains(pol, vor.l.split,
byid=T), ])
medianLine.spl <- gSimplify(medianLine, tol=tol); plot(medianLine.spl)
return(medianLine.spl)
}
myMedianLine <- findMedianLine(pol, tol=1000)
plot(pol); plot(myMedianLine, add=T)