Generating a Hexagonal Lattice Using Spatstat

Question

I'm analyzing the growth pattern of certain particles, and want to compare the point pattern with that of a perfectly hexagonal lattice with the same intensity (same number of points per unit of area). I've written a function that does this, but it has an inherent error and I'm not sure from where it originates. Essentially, after the function has run its course, it produces a perfectly hexagonal point pattern that does NOT have the correct number of particles - it is usually off by about 1-4%. If you copy and paste the following code into R, you'll see this - for this particular example, the error is 11.25% - the original point pattern has 71 particles and the perfectly hexagonal point pattern that was generated has 80 particles. This seems very odd, since the code, as you will see, is designed to place the particles a specific distance from one another, thereby creating a window of the same size as the original with the same number of particles.

The following is the code for the function I wrote to generate the hexagonal lattice.

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

Now, I'm relatively new to R, so it may very well be a syntax error somewhere in the code that has caused this error. Or, it may be that I have some fault in my train of thought with this process. However, I don't think that is likely, since my results have been quite close to what I've been trying for (only 1-4% error most of the time is fairly good).

In summary, what I would like help with is how to take a point pattern and create another point pattern of the same window size with the same number of particles, but a perfectly hexagonal point pattern.

Please don't hesitate to ask me to elucidate anything if you feel that something is unclear.

Thanks!

Answer 1

Forgive me if I am wrong, but I believe what you are trying to do is impossible (in the general case), given the constraints that you have shown in the example. Put simply, can you think of a way to draw 71 points in a hexagon pattern on a page with the same height and width as your window? I don't think such a pattern exists.

To explain further, consider the last line in your code:

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

Now, since your window in the same size as the original data, to get the same intensity, you would need exactly the same number of points (71). In your hexagon arrangement of points, you have x rows, each containing y points. But there are no integers x and y that multiply to 71.

That being said, if you "stretched" the window width a little, then half of your rows would contain one more point. That is a slightly looser constraint, but there is no guarantee that in the general case there would be a solution.

So, to get the exact same point intensity, you would need to be able to change the relative window size. You would need to stretch it out to add some white space and get a lower point intensity. That still might not work in the general case... but it might, I haven't worked it out. It might be easiest to start with just a plain grid, and then expand your code to hexagons.

---

Looking at your code, I noticed that you are using while loops when you could have been using the seq function. For example, if you wanted to generate all the x points from 0 to MaxXValue increasing by sqrt(3)*delta, simply do:

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

instead of that big while. There might be some errors here, but I think you could reduce your whole code to:

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

Finally, I noticed in your comments that you mention that the area of a hexagon is ( sqrt(3) / 2 ) * delta^2, but isn't it (3*sqrt(3)/2) * delta^2? `

---

I was interested in Josh O'Brien's comment, and decided to implement a rotation function to get the exact number of points desired. Here is the code:

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

Answer 2

Just for completeness I will add that spatstat now has a function hexgrid to generate hexagonal grid, and there is rotate.ppp which could be applied afterwards to rotate the pattern.