これが私の努力です。計画は次のとおりです。
- 線を密にする
- ドロネー三角形分割を計算する
- 中点を取得し、ポリゴン内にあるそれらの点を取得します
- 距離加重最小全域木を構築する
- グラフの直径のパスを見つける
手始めに高密度化コード:
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」関数は、1 つのリングだけであるべき SpatialPolygons のものから 1 つのリングの座標を取得するだけです。
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)))}
