这是我的努力。计划是:
- 使线条变密
- 计算 delaunay 三角剖分
- 取中点,取多边形中的那些点
- 构建距离加权最小生成树
- 找到它的图形直径路径
初学者的致密代码:
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)
}
我计算了从每个中点到多边形线的距离,而不是早期版本的缓冲,并且只取 a) 内部和 b) 距离边缘比内部点距离的 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)))}