找出近似规则网格点的子集的周长

Question

让我们考虑一组近似规则的二维网格。这些网格与相邻网格相邻（相邻网格具有一个或多个相同的顶点）。这是10个网格的示例，顶点坐标（经度，纬度）如下

A<-

        lon    lat
        [,1]     [,2]
  [1,] 85.30754 27.91250
  [2,] 85.32862 27.95735
  [3,] 85.34622 27.89880
  [4,] 85.36732 27.94364
  [5,] 85.34958 28.00202
  [6,] 85.38831 27.98830
  [7,] 85.38487 27.88508
  [8,] 85.40598 27.92991
  [9,] 85.42353 27.87134
 [10,] 85.44466 27.91616
 [11,] 85.42698 27.97456
 [12,] 85.46567 27.96081
 [13,] 85.48334 27.90239
 [14,] 85.50437 27.94703
 [15,] 85.48645 28.00502
 [16,] 85.52517 27.99123
 [17,] 85.52198 27.88862
 [18,] 85.54302 27.93325
 [19,] 85.56384 27.97745

上述样本点（顶点）集的散点图：

网格的构造如下图所示。

我的问题是如何得到周长（通过所有边界点的红色轮廓）??

注意：图1中红色圈出的点（1,3,7,9,10,13,17,18,19,16,15,12,11,6,5,2）是边界点。

注意：观察到格子的边长小于6000米，所有格子的对角线长度都大于6000米。

我正在使用 R 中 geosphere 包函数中的 distHaversine 来计算两点之间的距离。

Answer 1

解决此问题的一种方法是计算点周围的 alpha 外壳。

alphahull 包可以做到这一点。该软件包在 http://yihui.name/en/2010/04/alphahull-an-r-package-for-alpha-convex-hull/

处有出色的文档和完整的动画

阅读此文档非常值得，尤其是了解 alpha 参数的含义。

首先，复制您的数据

A <- read.table(header=TRUE, text="
lon    lat
[1,] 85.30754 27.91250
[2,] 85.32862 27.95735
[3,] 85.34622 27.89880
[4,] 85.36732 27.94364
[5,] 85.34958 28.00202
[6,] 85.38831 27.98830
[7,] 85.38487 27.88508
[8,] 85.40598 27.92991
[9,] 85.42353 27.87134
[10,] 85.44466 27.91616
[11,] 85.42698 27.97456
[12,] 85.46567 27.96081
[13,] 85.48334 27.90239
[14,] 85.50437 27.94703
[15,] 85.48645 28.00502
[16,] 85.52517 27.99123
[17,] 85.52198 27.88862
[18,] 85.54302 27.93325
[19,] 85.56384 27.97745")

现在计算并绘制 alpha 外壳。您必须提供 alpha 值。我不得不尝试找到一个足够小的值来捕获线。

library("alphahull")
hull <- with(A, ahull(lat, lon, alpha=0.033))
plot(hull)

Answer 2

简而言之：所有距离小于6000m的点对都以网格正方形的形式形成一个图形。构造该图，然后找到与正方形（大小为 4 的环）同构的所有子图。外部边缘出现的频率低于内部边缘，因为它们只是一个正方形的一部分（内部边缘由多个正方形共享）。因此找到所有内部边并丢弃它们，然后遍历结果图，这应该是一个简单的循环。

代码：

library(igraph); library(geosphere)

# main function
perimeterGrid <- function(pts, maxdist=6000, mindist=1){
    g = edgeP(makegrid(pts, maxdist=maxdist, mindist=mindist))
    loop = graph.dfs(minimum.spanning.tree(g),1)$order
    cbind(V(g)$x, V(g)$y)[loop,]
}

# haversine distance matrix
dmat <- function(pts){
    n=nrow(pts)
    do.call(rbind,lapply(1:n,function(i){distHaversine(pts[i,],pts)}))
}

# make the grid cells given a maxdist (and a mindist to stop self-self edges)    
makegrid <- function(pts, maxdist=6000,  mindist=1){
    dists = dmat(pts)
    g = graph.adjacency(dists<maxdist & dists>mindist,
        mode="undirected")
    V(g)$x=pts[,1]
    V(g)$y=pts[,2]
    g
}

# clever function that does the grid edge count etc
edgeP <- function(g){
    # find all the simple squares
    square=graph.ring(4)
    subs = graph.get.subisomorphisms.vf2(g,square)
    # expand all the edges
    subs = do.call(rbind, lapply(subs, function(s){
        rbind(s[1:2], s[2:3], s[3:4], s[c(4,1)])
    }))
    # make new graph of the edges of all the squares
    e = graph.edgelist(subs,directed=FALSE)
    # add the weight as the edge count
    E(e)$weight=count.multiple(e)

    # copy the coords from the source
    V(e)$x=V(g)$x
    V(e)$y=V(g)$y

    # remove multiple edges
    e=simplify(e)

    # internal edges now have weight 256.
    e = e - edges(which(E(e)$weight==256))
    # internal nodes how have degree 0
    e = e - vertices(degree(e)==0)
    return(e)
}

用法：

 plot(pts)
 polygon(perimeterGrid(pts),lwd=2)

结果：

警告：

这尚未在有孔的网格片段或网格单元仅在单个角点接触的网格片段上进行测试。也许那不可能发生。另外我不确定算法的效率如何，但看起来很快。

Answer 3

由于点集是网格，或者至少非常接近网格，因此使用 Delaunay 三角剖分似乎是一个不错的方法。

求 Delaunay 三角剖分。
找到两个有效的边角。如果网格是轴对齐的，那么这些将是垂直线和水平线。
根据边的角度将所有超出误差容限的边移除为网格上的有效边。几度的误差容限可能就足够了。
绕着剩余的边缘走一圈，找到周长。最左边或最右边的点是开始行走的好地方。

性能应该是 O(N log N)。

找出近似规则网格点的子集的周长

Find the perimeter of a subset of a near-regular grid of points

r

computational-geometry

r-grid