[R-sig-Geo] Skeletonization Process - Eliminating Adjacent Segments

[Mathieu Rajerison]

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)
>
>
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20130521/a814264c/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Rplot01.png
Type: image/png
Size: 10392 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20130521/a814264c/attachment.png>