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)

	[[alternative HTML version deleted]]

