Voronoi 镶嵌的共享边表

Question

我正在尝试基于由 spatstat 库的 dirichlet() 函数生成的 Voronoi 镶嵌（又名 Dirichlet 镶嵌或泰森多边形）创建多边形邻居表。例如，在下图中，右上角和右下角的瓷砖各有 2 个邻居，中间右侧的瓷砖有 4 个邻居，其余两个瓷砖各有 3 个邻居。我想捕获表中的邻居对，理想情况下捕获它们共享的边界线的长度：例如，'Tile1'、'Tile2'、'shared_edge_length'。

最初，我尝试使用intersect.tess()、intersect.own() 和polyclip 函数循环并比较镶嵌中的每一对多边形，但我猜这些不起作用，因为瓷砖是根据定义，尽管共享边缘，但在区域上不重叠。是否有一个简单的功能来实现这一点（替代方法是循环遍历$bdry 点）？ regeos 包似乎有 gTouches，但我找不到 spatstat 的类似内容。

这是我目前不起作用的方法：

library(spatstat)
points <- ppp(x=c(-77.308703, -77.256582, -77.290600,  -77.135668, -77.097144),
         y=c(39.288603, 39.147019, 39.372818, 39.401898, 39.689203),
         window=owin(xrange=c(-77.7,-77), yrange=c(39.1, 39.7)))
vt <- dirichlet(points) # Dirichlet tesselation
plot(vt)

tilesA <- tiles(vt)
n_tiles <- length(tilesA)
boundary_calcs <- data.frame('area1_id'=numeric(), 'area2_id'=numeric(), 'neighbor'=logical()) # Store boundary pairs
for (i in 1:n_tiles) {
  for (j in 1:n_tiles) {
    intersection <- intersect.owin(tilesA[[i]], tilesA[[j]], fatal=FALSE) # does not work
    if (!is.empty(intersection)) {
      boundary_calcs[nrow(boundary_calcs)+1, ] <- c(i, j, TRUE) # add to data table as new row
} } }

Answer 1

以下功能由包作者提供。它使用deldir() 函数的dirsgs 结构输出镶嵌中每条线的开始/结束坐标以及点索引的事实。这些可以转换为psp 线段模式，可以使用lengths.psp() 轻松提供每个线段的长度。下面的代码生成了一个表格，其中 7 条边中的每条边都有一行，可以在上面的图中看到。

library(spatstat)
library(deldir)
points <- ppp(x=c(-77.308703, -77.256582, -77.290600,  -77.135668, -77.097144),
         y=c(39.288603, 39.147019, 39.372818, 39.401898, 39.689203),
         window=owin(xrange=c(-77.7,-77), yrange=c(39.1, 39.7)))

sharededge <- function(X) {
   verifyclass(X, "ppp")
   Y <- X[as.rectangle(X)]
   dX <- deldir(Y)
   DS <- dX$dirsgs
   xyxy <- DS[,1:4]
   names(xyxy) <- c("x0","y0","x1","y1")
   sX <- as.psp(xyxy,window=dX$rw)
   marks(sX) <- 1:nobjects(sX)
   sX <- sX[as.owin(X)]
   tX <- tapply(lengths.psp(sX), marks(sX), sum)
   jj <- as.integer(names(tX))
   ans <- data.frame(ind1=DS[jj,5], 
                     ind2=DS[jj,6], 
                     leng=as.numeric(tX))
   return(ans)
}

shared_edge_lengths <- sharededge(points)

输出存储在shared_edge_lengths:

  ind1 ind2       leng
1    2    1 0.17387212
2    3    1 0.13444458
3    4    1 0.05791519
4    4    2 0.10039321
5    4    3 0.25842530
6    5    3 0.09818828
7    5    4 0.17162429

Answer 2

我刚刚问了 deldir 包的作者 Rolf Turner，他提供了下面的代码来查找每个图块的邻居数量。我还没有看到您希望“理想地捕获他们共享的边界线的长度”。

x <- c(-77.308703, -77.256582, -77.290600,  -77.135668, -77.097144)
y <- c(39.288603, 39.147019, 39.372818, 39.401898, 39.689203)
rw <- c(-77.7,-77,39.1,39.7)
require(deldir)
dxy <- deldir(x,y,rw=rw)

dxy$summary$n.tside