Find the perimeter of a subset of a near-regular g

问题:

Let us consider a set of near-regular grids in 2-D. These grids are adjacent (neighbouring grids have one or more same vertices) to the neighbouring grids. Here are the sample of 10 grids with the coordinates of the vertices (longitude,latitude) are as follows

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

The scatter-plot of the above sample set of points (vertices):

The grids are constructed as in the following picture.

My question is how to get the perimeter (red contour passing through all boundary points)??

Note that: The red encircled points (1,3,7,9,10,13,17,18,19,16,15,12,11,6,5,2) in figure 1 are the boundary points.

Note: It is observed the sides of the grids are less than 6000 metres and length of diagonals of all grids are more than 6000 metres.

I am using distHaversine from the geosphere package function in R to find the distance between two points.

回答1:

In outline: all pairs of points closer than 6000m form a graph in the form of grid squares. Construct that graph, and then find all subgraphs isomorphic to a square (a loop of size 4). The external edges will appear less often than internal edges since they are only part of one square (internal edges are shared by multiple squares). Hence find all the internal edges and drop them, then traverse the resulting graph which should be a simple loop.

Code:

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)
}

Usage:

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

Results:

Warnings:

This is untested on grid fragments with holes or where grid cells are only touching at a single corner points. Maybe that can't happen. Also I'm not sure what the efficiency of the algorithms are, but it seems pretty quick.

回答2:

One way of solving this problem is to compute the alpha hull around your points.

The alphahull package can do this. The package has excellent documentation, complete with animation, at http://yihui.name/en/2010/04/alphahull-an-r-package-for-alpha-convex-hull/

It is well worth taking a look at this documentation, in particular to understand the meaning of the alpha parameter.

First, replicate your data

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")

Now compute and plot the alpha hull. You have to provide the alpha value. I had to experiment to find a value small enough to capture the ouline.

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

回答3:

Since the point set is a grid, or at least very close to a grid, seems like a good approach would be to make use of the Delaunay triangulation.

1. Find the Delaunay triangulation.
2. Find the two valid edge angles. If the grid was axis aligned, then these would be vertical and horizontal lines.
3. Remove all edges that exceed error tolerances to be a valid edge on the grid based on the angle of the edge. An error tolerance of a few degrees is probably enough.
4. Do a walk around remaining edges to find the perimeter. The left-most or right-most point is a good place to start the walk.

Performance should be O(N log N).