Вот мои усилия. План:
- уплотнение линий
- вычислить триангуляцию Делоне
- возьмите средние точки и те точки, которые находятся в многоугольнике
- построение взвешенного на расстоянии минимального остовного дерева
- найди свой график диаметра пути
Уплотняющий код для стартеров:
densify <- function(xy,n=5){
## densify a 2-col matrix
cbind(dens(xy[,1],n=n),dens(xy[,2],n=n))
}
dens <- function(x,n=5){
## densify a vector
out = rep(NA,1+(length(x)-1)*(n+1))
ss = seq(1,length(out),by=(n+1))
out[ss]=x
for(s in 1:(length(x)-1)){
out[(1+ss[s]):(ss[s+1]-1)]=seq(x[s],x[s+1],len=(n+2))[-c(1,n+2)]
}
out
}
А теперь основное блюдо:
simplecentre <- function(xyP,dense){
require(deldir)
require(splancs)
require(igraph)
require(rgeos)
### optionally add extra points
if(!missing(dense)){
xy = densify(xyP,dense)
} else {
xy = xyP
}
### compute triangulation
d=deldir(xy[,1],xy[,2])
### find midpoints of triangle sides
mids=cbind((d$delsgs[,'x1']+d$delsgs[,'x2'])/2,
(d$delsgs[,'y1']+d$delsgs[,'y2'])/2)
### get points that are inside the polygon
sr = SpatialPolygons(list(Polygons(list(Polygon(xyP)),ID=1)))
ins = over(SpatialPoints(mids),sr)
### select the points
pts = mids[!is.na(ins),]
dPoly = gDistance(as(sr,"SpatialLines"),SpatialPoints(pts),byid=TRUE)
pts = pts[dPoly > max(dPoly/1.5),]
### now build a minimum spanning tree weighted on the distance
G = graph.adjacency(as.matrix(dist(pts)),weighted=TRUE,mode="upper")
T = minimum.spanning.tree(G,weighted=TRUE)
### get a diameter
path = get.diameter(T)
if(length(path)!=vcount(T)){
stop("Path not linear - try increasing dens parameter")
}
### path should be the sequence of points in order
list(pts=pts[path+1,],tree=T)
}
Вместо буферизации в более ранней версии я вычисляю расстояние от каждой средней точки до линии многоугольника и беру только точки, которые а) находятся внутри и б) дальше от края, чем 1,5 от расстояния изнутри точка, наиболее удаленная от края.
Проблемы могут возникнуть, если многоугольник изгибается сам по себе, с длинными сегментами и без уплотнения. В этом случае граф является деревом и код сообщает об этом.
В качестве теста я оцифровал линию (s, объект SpatialLines), буферизовал ее (p), затем вычислил осевую линию и наложил их:
s = capture()
p = gBuffer(s,width=0.2)
plot(p,col="#cdeaff")
plot(s,add=TRUE,lwd=3,col="red")
scp = simplecentre(onering(p))
lines(scp$pts,col="white")
Функция 'onering' просто получает координаты одного кольца от объекта SpatialPolygons, который должен быть только одним кольцом:
onering=function(p){p@polygons[[1]]@Polygons[[1]]@coords}
Захват пространственных линий с помощью функции «захват»:
capture = function(){p=locator(type="l")
SpatialLines(list(Lines(list(Line(cbind(p$x,p$y))),ID=1)))}