title: Generating a Hexagonal Lattice Using Spatstat
body: <p>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.  </p>
  
  <p>The following is the code for the function I wrote to generate the hexagonal lattice.</p>
  
  <pre><code>library(spatstat)
  
  data(swedishpines)
  
  swedishpines.df &lt;- as.data.frame(swedishpines)
  
  MaxXValue &lt;- max(swedishpines.df[1])
  MaxYValue &lt;- max(swedishpines.df[2])
  #The above two lines dictate the window size
  NumberOfParticles &lt;- nrow(swedishpines.df[1])
  #Number of entries = number of particles
  
  #calculate delta
  intensity &lt;- NumberOfParticles / (MaxXValue*MaxYValue)
  #Intensity ---&gt; 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 &lt;- sqrt(2/(intensity*(sqrt(3))))
  #This is derived from the equation for the area of a regular hexagon. 
  </code></pre>
  
  <p><img src="https://i.stack.imgur.com/2Zcpc.jpg" alt="Diagram of a unit cell of the hexagonal lattice"></p>
  
  <pre><code>#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 &lt;- c()
  ycoords &lt;- 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.
  </code></pre>
  
  <p><img src="https://i.stack.imgur.com/z0UBO.jpg" alt="First set of points"></p>
  
  <pre><code>y_shift &lt;- 0
  while (y_shift &lt; MaxYValue) {
  
      x &lt;- c(0)
      x_shift &lt;- 0 + delta
      count &lt;- 0
  
      while (x_shift &lt; MaxXValue) {
          x &lt;- append(x, x_shift)
          x_shift &lt;- x_shift + delta
          count &lt;- count + 1
      }
  
      y &lt;- c(y_shift)
  
      for (i in seq(0,(count-1))) {
          y &lt;- append(y, y_shift)
      }
  
      y_shift &lt;- y_shift + sqrt(3)*delta
      xcoords &lt;- append(xcoords,x)
      ycoords &lt;- 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. 
  </code></pre>
  
  <p><img src="https://i.stack.imgur.com/dOXyF.jpg" alt="Second set of points"></p>
  
  <pre><code>y_shift &lt;- 0 + (delta*(sqrt(3)))/2
  while (y_shift &lt; MaxYValue) {
  
      x &lt;- c(0 + (1/2)*delta)
      x_shift &lt;- 0 + (1/2)*delta + delta
      count &lt;- 0
  
      while (x_shift &lt; MaxXValue) {
          x &lt;- append(x, x_shift)
          x_shift &lt;- x_shift + delta
          count &lt;- count + 1
      }
  
      y &lt;- c(y_shift)
  
      for (i in seq(0,(count-1))) {
          y &lt;- append(y, y_shift)
      }
  
      y_shift &lt;- y_shift + sqrt(3)*delta
      xcoords &lt;- append(xcoords,x)
      ycoords &lt;- append(ycoords,y)
  
  }
  
  hexplot &lt;- ppp(xcoords, ycoords, window=owin(c(0,MaxXValue),c(0,MaxYValue)))
  </code></pre>
  
  <p><img src="https://i.stack.imgur.com/qd3RH.jpg" alt="Both sets combined"></p>
  
  <p>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).  </p>
  
  <p>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.</p>
  
  <p>Please don't hesitate to ask me to elucidate anything if you feel that something is unclear.</p>
  
  <p>Thanks!</p>