r - 使用 Spatstat 生成六角晶格

Question

我正在分析某些粒子的生长模式，并想将点模式与强度相同(每单位面积的点数相同)的完美六角晶格的点模式进行比较。我已经编写了一个执行此操作的函数，但它有一个固有错误，我不确定它的来源。本质上，在函数运行完之后，它会生成一个完美的六边形点图案，但没有正确数量的粒子——通常会偏离大约 1-4%。如果将以下代码复制并粘贴到 R 中，您将看到这一点 - 对于这个特定示例，错误为 11.25% - 原始点图案有 71 个粒子，生成的完美六边形点图案有 80 个粒子。这看起来很奇怪，因为正如您将看到的，代码旨在将粒子彼此放置特定距离，从而创建一个与原始窗口大小相同且粒子数量相同的窗口。

下面是我写的生成六角点阵的函数代码。

library(spatstat)

data(swedishpines)

swedishpines.df <- as.data.frame(swedishpines)

MaxXValue <- max(swedishpines.df[1])
MaxYValue <- max(swedishpines.df[2])
#The above two lines dictate the window size
NumberOfParticles <- nrow(swedishpines.df[1])
#Number of entries = number of particles

#calculate delta
intensity <- NumberOfParticles / (MaxXValue*MaxYValue)
#Intensity ---> particles / unit^2
    #Area = ( sqrt(3) / 2 ) * delta^2
    #Now - in substituting intensity in for area, it is key to recognize 
    #that there are 3 particles associated with each hexagonal tile.
    #This is because each particle on the border of the tile is really 1/3 of a 
    #a particle due to it being associated with 3 different hexagonal tiles.  
    #Since intensity = 3 Particles / Area, 
delta <- sqrt(2/(intensity*(sqrt(3))))
#This is derived from the equation for the area of a regular hexagon. 

#xcoords and ycoords represent the x and y coordintes of all of the generated points.  The 'x' and 'y' are temporary holders for the x and y coordinates of a single horizontal line of points (they are appended to xcoords and ycoords at the end of each while loop).  

xcoords <- c()
ycoords <- c()

#The following large while loop calculates the coordinates of the first set of points that are vertically aligned with one another. (alternating horizontal lines of particles)  This is shown in the image below.

y_shift <- 0
while (y_shift < MaxYValue) {

    x <- c(0)
    x_shift <- 0 + delta
    count <- 0

    while (x_shift < MaxXValue) {
        x <- append(x, x_shift)
        x_shift <- x_shift + delta
        count <- count + 1
    }

    y <- c(y_shift)

    for (i in seq(0,(count-1))) {
        y <- append(y, y_shift)
    }

    y_shift <- y_shift + sqrt(3)*delta
    xcoords <- append(xcoords,x)
    ycoords <- append(ycoords,y)

}

#The following large while loop calculates the coordinates of the second set of points that are vertically aligned with one another. This is shown in the image below. 

y_shift <- 0 + (delta*(sqrt(3)))/2
while (y_shift < MaxYValue) {

    x <- c(0 + (1/2)*delta)
    x_shift <- 0 + (1/2)*delta + delta
    count <- 0

    while (x_shift < MaxXValue) {
        x <- append(x, x_shift)
        x_shift <- x_shift + delta
        count <- count + 1
    }

    y <- c(y_shift)

    for (i in seq(0,(count-1))) {
        y <- append(y, y_shift)
    }

    y_shift <- y_shift + sqrt(3)*delta
    xcoords <- append(xcoords,x)
    ycoords <- append(ycoords,y)

}

hexplot <- ppp(xcoords, ycoords, window=owin(c(0,MaxXValue),c(0,MaxYValue)))

现在，我对 R 还比较陌生，所以很可能是代码中某个地方的语法错误导致了这个错误。或者，可能是我在这个过程中的思路有问题。但是，我认为这不太可能，因为我的结果与我一直在尝试的结果非常接近(大多数时候只有 1-4% 的错误是相当不错的)。

总而言之，我需要帮助的是如何采用一个点图案并创建另一个具有相同窗口大小、具有相同数量粒子但完美的六边形点图案的点图案。

如果您觉得有什么不清楚的地方，请不要犹豫，让我解释清楚。

谢谢!

Answer 1

如果我错了，请原谅我，但我相信鉴于您在示例中显示的约束，您尝试做的事情是不可能的(在一般情况下)。简而言之，你能想出一种方法，在与你的窗口高度和宽度相同的页面上以六边形图案绘制 71 个点吗？我认为这种模式不存在。

为了进一步解释，请考虑代码中的最后一行:

hexplot <- ppp(xcoords, ycoords, window=owin(c(0,MaxXValue),c(0,MaxYValue)))

现在，由于您的窗口与原始数据大小相同，要获得相同的强度，您需要完全相同的点数 (71)。在你的六边形点排列中，你有 x 行，每行包含 y 个点。但是没有整数 x 和 y 可以乘以 71。

也就是说，如果您稍微“拉伸(stretch)”窗口宽度，那么一半的行将包含一个点。这是一个稍微宽松的约束，但不能保证在一般情况下会有解决方案。

因此，要获得完全相同的点强度，您需要能够更改相对窗口大小。您需要将其拉伸(stretch)以添加一些空白并获得较低的点强度。在一般情况下，这仍然可能行不通......但它可能，我还没有解决。从简单的网格开始，然后将代码扩展到六边形可能是最简单的。

---

查看您的代码，我注意到您正在使用 while 循环，而您本可以使用 seq 函数。例如，如果您想生成从 0 到 MaxXValue 的所有 x 点，增加 sqrt(3)*delta，只需执行以下操作:

x<-seq(0,MaxXValue,by=delta)

而不是那个大的 while。这里可能有一些错误，但我认为您可以将整个代码缩减为:

library(spatstat)
data(swedishpines)
swedishpines.df <- as.data.frame(swedishpines)
MaxXValue <- max(swedishpines.df[1])
MaxYValue <- max(swedishpines.df[2])
NumberOfParticles <- nrow(swedishpines.df)
intensity <- NumberOfParticles / (MaxXValue*MaxYValue)
delta <- sqrt(2/(intensity*(sqrt(3))))
x<-seq(0,MaxXValue,by=delta)
y<-seq(0,MaxYValue,by=sqrt(3)*delta)
first.coords=merge(x,y) # Find every combo of x and y
x=seq(delta/2,MaxXValue,by=delta)
y=delta*sqrt(3)/2 + (delta*sqrt(3)*seq(0,length(x)/2))
second.coords=merge(x,y)
coords=rbind(first.coords,second.coords)
ppp(coords$x, coords$y, window=owin(c(0,MaxXValue),c(0,MaxYValue)))

最后，我注意到您在评论中提到六边形的面积是 ( sqrt(3)/2 ) * delta^2，但不是吗 (3 *sqrt(3)/2) * delta^2?`

---

我对 Josh O'Brien 的评论很感兴趣，并决定实现旋转函数以获得所需的确切点数。这是代码:

# Include all above code
rotate=function(deg) {
    r=matrix(c(cos(deg),-sin(deg),sin(deg),cos(deg)),byrow=T,nrow=2)
    rotated.coords=data.frame(t(r %*% t(as.matrix(coords))))
    names(rotated.coords)=c('x','y')
    rotated.coords
}

rotate.optim=function(deg) {
    rotated.coords=rotate(deg)
    abs(NumberOfParticles-length(suppressWarnings(ppp(rotated.coords$x, rotated.coords$y, window=owin(c(0,MaxXValue),c(0,MaxYValue)))$x)))
}

o=optim(0,rotate.optim,lower=0,upper=2*pi,method='Brent')
rotated.coords=rotate(o$par)
rotated.coords.window=rotated.coords[rotated.coords$x >= 0 & rotated.coords$y >= 0 & rotated.coords$x <= MaxXValue & rotated.coords$y <= MaxYValue,]
final=ppp(rotated.coords.window$x,rotated.coords.window$y,window=owin(c(0,MaxXValue),c(0,MaxYValue)))
plot(final)