<div dir="ltr">Hi,<div><br></div><div>Finally, I inspired myself from a code by Barry Rownlingson intially dedicated to calculate the shortest path.</div><div style>But in my cas, it consisted in calculating the longest one.</div>
<div style><br></div><div style>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.</div><div style>
<br></div><div style><div>library(rgdal)</div><div>library(spatstat)</div><div>library(rgeos)</div><div>library(maptools)</div><div>library(igraph)</div><div><br></div><div>setwd("F:/METHODES/1305_SKELETON/")</div>
<div><br></div><div>pol <- readOGR("IN/polygone.shp", "polygone")</div><div><br></div><div># FUNCTION</div><div><br></div><div>findMedianLine <- function(pol, eps, tol) {</div><div><br></div><div>
l <- as(as(pol, "SpatialLines"), "psp")</div><div><br></div><div> pts <- pointsOnLines(l, eps=eps)</div><div><br></div><div> vor <- dirichlet(pts)</div><div> </div><div> </div><div> # LINES EXTRACTION</div>
<div><br></div><div> vor.l <- as(as(vor, "SpatialPolygons"), "SpatialLines")</div><div><br></div><div> vor.l <- gIntersection(vor.l, vor.l)@lines[[1]]@Lines</div><div><br></div><div> vor.split <- SpatialLines(lapply(1:length(vor.l), function(x) Lines(list(vor.l[[x]]), ID=as.character(x)))) </div>
<div> </div><div> <a href="http://vor.l.in">vor.l.in</a> <- vor.l.split[gContains(pol, vor.split, byid=T), ]</div><div> </div><div> </div><div> # GRAPH CONSTRUCTION</div><div> </div><div> edges = do.call(rbind,lapply(vor.l.in@lines,function(ls){as.vector(t(ls@Lines[[1]]@coords))}))</div>
<div> lengths = sqrt((edges[,1]-edges[,3])^2+(edges[,2]-edges[,4])^2)</div><div> </div><div> froms = paste(edges[,1],edges[,2])</div><div> tos = paste(edges[,3],edges[,4])</div><div> </div><div> g = graph.edgelist(cbind(froms,tos),directed=FALSE)</div>
<div> E(g)$weight=lengths</div><div> </div><div> </div><div> # LONGEST PATH</div><div> </div><div> sel <- get.diameter(g)</div><div> </div><div> vor.sel <- vor.l.in@lines[sel-1]</div><div> </div><div> spine <- SpatialLines(vor.sel) </div>
<div> </div><div> spine.spl <- gSimplify(gLineMerge(spine), tol=tol)</div><div><br></div><div> return(spine.spl)</div><div>}</div><div><br></div><div><br></div><div>myMedianLine <- findMedianLine(pol, eps=1000, tol=1000)</div>
</div></div><div class="gmail_extra"><br><br><div class="gmail_quote">2013/5/17 Mathieu Rajerison <span dir="ltr"><<a href="mailto:mathieu.rajerison@gmail.com" target="_blank">mathieu.rajerison@gmail.com</a>></span><br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">Hi,<div><br><div><br></div><div>I've written a code to skeletonize an area shape.</div><div><br></div>
<div>The example could be a watershed from which we'd like to compute the "guessed" stream, or find the median axis of a shape.</div>
<div><br></div><div>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.</div>
<div><br></div><div>Has anyone idea to deal with that?</div><div><br></div><div>I apoplogize in advance if it doesn't work in all cases (holes, rings, etc...).</div><div><br></div><div>
<div>library(rgdal)</div><div>library(spatstat)</div><div>library(rgeos)</div><div>library(maptools)</div><div><br></div><div>setwd("F:/METHODES/1305_SKELETON/")</div><div><br></div><div>pol <- readOGR("IN/polygone.shp", "polygone")</div>
<div><br></div><div># FUNCTION</div><div><br></div><div>findMedianLine <- function(pol, tol) {</div><div><br></div><div> l <- as(as(pol, "SpatialLines"), "psp")</div><div><br></div><div> pts <- pointsOnLines(l, eps=10000)</div>
<div><br></div><div> vor <- dirichlet(pts)</div><div><br></div><div> vor.l <- as(as(vor, "SpatialPolygons"), "SpatialLines")</div><div><br></div><div> vor.lines <- gIntersection(vor.l, vor.l)@lines[[1]]@Lines</div>
<div><br></div><div> vor.l.split <- SpatialLines(lapply(1:length(vor.lines), function(x) Lines(list(vor.lines[[x]]), ID=as.character(x))))</div><div><br></div><div> medianLine <- gLineMerge(vor.l.split[gContains(pol, vor.l.split, byid=T), ])</div>
<div><br></div><div> medianLine.spl <- gSimplify(medianLine, tol=tol); plot(medianLine.spl)</div><div><br></div><div> return(medianLine.spl)</div><div>}</div><div><br></div><div><br></div><div>myMedianLine <- findMedianLine(pol, tol=1000)</div>
<div><br></div><div>plot(pol); plot(myMedianLine, add=T)</div></div><div><br></div><div><br></div><div><br></div><div><br></div></div></div>
</blockquote></div><br></div>