Poisson Disk Sampling

This article originally appeared in Dev.Mag Issue 21, released in March 2008.

One way to populate large worlds with objects is to simply place objects on a grid, or randomly. While fast and easy to implement, both these methods result in unsatisfying worlds: either too regular or too messy. In this article we look at an alternative algorithm that returns a random set of points with nice properties:

  • the points are tightly packed together; but
  • no closer to each other than a specified minimum distance.

Figure 1 shows an example of such a set, which is called Poisson-disk sample set. For comparison, two other sample sets are also shown: a uniform random sample, and a jittered grid sample.

Poisson disk sampling has many applications in games:

  • random object placement;
  • sampling for graphics applications;
  • procedural texture algorithms; and
  • mesh algorithms.

In this article we will look mostly at object placement and briefly at texture generation.

Figure 1

Figure 1

Implementation

There are several algorithms for producing a Poisson disk sample set. The one presented here is easy to implement, and runs reasonably fast. It is also easy adapted for specific applications (described in the next section).

The basic idea is to generate points around existing points, and to check whether they can be added so that they don’t disturb the minimum distance requirement. A grid is used to perform fast lookups of points. Two lists keep track of points that are being generated, and those that needs processing.

Here are the details:

  1. A grid is created such that every cell will contain at most one sampling point. If points are at least distance r from each other, then the cell size must be r/2?. The ill-rendered symbol ? is pi.
  2. The first point is randomly chosen, and put in the output list, processing list and grid.
  3. Until the processing list is empty, do the following:
    1. Choose a random point from the processing list.
    2. For this point, generate up to k points, randomly selected from the annulus surrounding the point. You can choose k – a value of 30 gives good results. In general, larger values give tighter packings, but make the algorithm run slower. For every generated point:
      1. Use the grid to check for points that are too close to this point. See below for more detail.
      2. If there is none, add the point to the output list, processing list, and grid.
    3. Remove the point from the processing list.
  4. Return the output list as the sample points.

Here is how all this look in pseudo code:

generate_poisson(width, height, min_dist, new_points_count)
{
  //Create the grid
  cellSize = min_dist/sqrt(2);

  grid = Grid2D(Point(
    (ceil(width/cell_size),         //grid width
    ceil(height/cell_size))));      //grid height

  //RandomQueue works like a queue, except that it
  //pops a random element from the queue instead of
  //the element at the head of the queue
  processList = RandomQueue();
  samplePoints = List();

  //generate the first point randomly
  //and updates 

  firstPoint = Point(rand(width), rand(height));

  //update containers
  processList.push(firstPoint);
  samplePoints.push(firstPoint);
  grid[imageToGrid(firstPoint, cellSize)] = firstPoint;

  //generate other points from points in queue.
  while (not processList.empty())
  {
    point = processList.pop();
    for (i = 0; i < new_points_count; i++)
    {
      newPoint = generateRandomPointAround(point, min_dist);
      //check that the point is in the image region
      //and no points exists in the point's neighbourhood
      if (inRectangle(newPoint) and
        not inNeighbourhood(grid, newPoint, min_dist,
          cellSize))
      {
        //update containers
        processList.push(newPoint);
        samplePoints.push(newPoint);
        grid[imageToGrid(newPoint, cellSize)] =  newPoint;
      }
    }
  }
  return samplePoints;
}

The grid coordinates of a point can be easily calculated:

imageToGrid(point, cellSize)
{
  gridX = (int)(point.x / cellSize);
  gridY = (int)(point.y / cellSize);
  return Point(gridX, gridY);
}

Figure 2 shows how a random point (red) is selected in the annulus around an existing point (blue). Two parameters determine the new point’s position: the angle (randomly chosen between 0 and 360 degrees), and the distance from the original point (randomly chosen between the minimum distance and twice the minimum distance). In pseudo code:

generateRandomPointAround(point, mindist)
{ //non-uniform, favours points closer to the inner ring, leads to denser packings
  r1 = Random.nextDouble(); //random point between 0 and 1
  r2 = Random.nextDouble();
  //random radius between mindist and 2 * mindist
  radius = mindist * (r1 + 1);
  //random angle
  angle = 2 * PI * r2;
  //the new point is generated around the point (x, y)
  newX = point.x + radius * cos(angle);
  newY = point.y + radius * sin(angle);
  return Point(newX, newY);
}
Figure 2 - Generating a new sample point.

Figure 2 - Generating a new sample point.

Before a newly generated point is admitted as a sample point, we have to check that no previously generated points are too close. Figure 3 shows a piece of the grid. The red dot is a potential new sample point. We have to check for existing points in the region contained by the red circles (they are the circles at the corners of the cell of the new point). The blue squares are cells that are partially or completely covered by a circle. We need only check these cells. However, to simplify the algorithm, we check all 25 cells.

Here is the pseudo code:

inNeighbourhood(grid, point, mindist, cellSize)
{
  gridPoint = imageToGrid(point, cellSize)
  //get the neighbourhood if the point in the grid
  cellsAroundPoint = squareAroundPoint(grid, gridPoint, 5)
  for every cell in cellsAroundPoint
    if (cell != null)
      if distance(cell, point) < mindist
        return true
  return false
}
Figure 3 - Checking the neighbourhood of a potential sample point.

Figure 3 - Checking the neighbourhood of a potential sample point.

Figure 4 - Spheres placed at points in a Poisson disk sample of 3D space

Figure 4 - Spheres placed at points in a Poisson disk sample of 3D space.

Implementation for 3D

The algorithm can easily be modified for 3D:

  • Change all points to 3D points.
  • Change the grid to a 3D grid. The neighbourhood of a point is now a cube of 125 cells around the cell of the point.
  • Change the code to generate a new point to the following:
generateRandomPointAround(point, minDist)
{ //non-uniform, leads to denser packing.
  r1 = Random.nextDouble(); //random point between 0 and 1
  r2 = Random.nextDouble();
  r3 = Random.nextDouble();
  //random radius between mindist and 2* mindist
  radius = mindist * (r1 + 1);
  //random angle
  angle1 = 2 * PI * r2;
  angle2 = 2 * PI * r3;
  //the new point is generated around the point (x, y, z)
  newX = point.x + radius * cos(angle1) * sin(angle2);
  newY = point.y + radius * sin(angle1) * sin(angle2);
  newZ = point.z + radius * cos(angle2);
  return Point(newX, newY, newZ);
}

Applications

Object Placement

Placing objects at the positions of a Poisson disk sample set is the simplest way to use this algorithm in games (Figure 5). Ideally, the algorithm can be built into your level editing tool with features that allows the artist to select the region to fill, and the models to fill them with.

Figure 5 - An example with shrubs placed at Poisson disk sample points.

Figure 5 - An example with shrubs placed at Poisson disk sample points.

Figure 5 – An example with shrubs placed at Poisson disk sample points.

One important variation of a Poisson sample set is one where the minimum distance between points is not constant, but varies across the image. In this variation, we feed the algorithm a greyscale image, which is used to modulate the minimum distance between points.

To make this work, you need to modify the algorithm as follows:

  • The grid must take a list of points.
  • The cell size must be computed from the maximumthe radius can be.
  • The neighbourhood function should iterate through all points in every cell of a point’s neighbourhood.
  • Where a new point is generated, check the grey scale value of the image at that point, and calculate a minimum distance from the grey scale value:
    min_dist = min_radius + grey * (max_radius - min_radius)

As an example, you can use Perlin noise to drive the minimum distance, giving you interesting clusters of objects (Figure 6). This method is especially useful for generating a field of plants.

Figure 6 - Poisson disk sample, where the minimum distance is driven by Perlin noise.

Figure 6 - Poisson disk sample, where the minimum distance is driven by Perlin noise.

When using different radii as explained above, you might run into some problems. Take the following precautions:

  • Ensure that your minimum radius is not too small. A very small radius might produce millions of points; a zero radius might prevent the algorithm from ever completing.
  • Built in a bail-out feature in the algorithm, forcing it to end after a certain number of points have been reduced. Make it a function parameter, so that it can be modified according to the purpose. If you make a tool on top of this algorithm, also expose it to the user in the tool.
  • Ensure that your maximum radius is not too big: if it is, no or very few points will be produced. This might seem obvious, but it can go wrong in a subtle way. Say, for example, you want to create a fall-off effect around a certain point. It would be natural to define your minimum distance array as follows:
    min_dist[i, j] = dist((i, j), (x0, y0))

    But because new points are generated at exactly this distance, many more points are excluded than expected, leading to a rather poor result. A better sample can be obtained by using the square root of the distance. (See Figure 7.)

Figure 7

Figure 7

In another situation, sections of large radii might be too close to other sections of large radii, so that no points are produced in sections of small radii (see Figure 8).

Figure 8 - The high radius is too great; some low radius regions are skipped.

Figure 8 - The high radius is too great; some low radius regions are skipped.

  • For best results, your radius should vary smoothly across the rectangle. The algorithm respects sharp edges only ruggedly – see the checker board examples in Figure 9.
  • Beware of the bleeding effect, as seen in Figure 9. You might want to run a dilation filter (use Photoshop or Gimp) on the radius grid before you do the sampling to compensate for this. Ideally the dilation should be a function of the minimum radius in a region, but in many cases you can used a fixed dilation.
Figure 9

Figure 9

Texture Generation

Poisson samples can also be used to generate certain textures. In Figure 10, the droplets on the bottle and glass have been created by combining three layers of Poisson disk sample points.

The modified algorithm creates three layers of points, each layer with a smaller minimum distance between points than the next. In every layer, the algorithm checks that there are no circles in the previous layers. To do this, the look-up grids for all generated layers are used in a final step to eliminate unwanted points.

The local minimum distance of a every sample point is stored, so that it can be used as a radius to draw a circle once all points have been found.

The raw texture is then further processed by clever artists through filters and shaders to produce the final result.

Figure 10.

Figure 10.

Poisson disk samples also form the basis of the procedural textures shown in Figure 11.

The first texture was produced by painting filled translucent circles for every Poisson sample point, for three separate samples. A different colour was used for every sample set, with small random variations.

The second texture was produced by painting circles for two sample sets; one with filled circles, the other with outlines only.

The third texture was created by painting daisies (randomly scaled and rotated) onto an existing grass texture.

Figure 11

Figure 11

Figure 11 - Procedural textures generated from Poisson disk samples.

Figure 11 - Procedural textures generated from Poisson disk samples.

Download

(Some of these links were added after the article was first published).

Resources

  • Pingback: Cellular Automata for Simulation in Games · code-spot

  • Roy

    Thanks for the great tutorial, this was very helpful indeed. I noticed one copy-paste error in “generateRandomPointAround”. You should use r3 as well ;-)

    • Herman Tulleken

      Thanks for spotting that one, Roy. It’s fixed now.

  • http://www.hlcao.com hailengc

    Nice article,i’ll study it in detail~~thanks a lot

  • http://www.vision.ime.usp.br/~creativision Jorge Leandro

    Hi, Tulleken

    Thanks for the tutorial. I’m still wondering what the grid is indeed: either a small matrix used to check the neighbourhood of the current point or a matrix with the same size of the great matrix to be sampled.

    Anyway, the java implementation worked out. Also, I managed to call it from Matlab.

    By the way, are you the author of the Java implementation?
    Regards from Brazil.

    • Herman Tulleken

      Hi Jorge,

      The matrix is just a partition of the sample space to perform fast lookups, so it is 25 cells from the big grid.

      We put the points in the grid as we generate them; when we want to examine the points in the set close to any given point, we only have to check a few points, instead of the entire set. What is nice about this scheme is that it is so simple – no weird tree structures as is often the case with spatial partitioning.

      The size covered by the cells is only a function of the minimum distance we want between points – in particular, it is not a function of the granularity of the points in the original space.

      The algorithm is not very concerned with the actual pixels of the space to sample – points are discretized either as they are generated, or as a final step before the list is returned, or even by the caller if you choose so.

      Yup, I am the author of the code… I hope your question is not a prelude to reporting a bug :P

      Thanks for your comment.

  • John Adams

    Thank you for the great explanation! When running the Python code, I was just wondering – in the uniform poisson based sampling, is there a suggested way to control the number of points generated? That is, if I want to say, generate 20 points in total for instance.

  • mammad

    Thank you for your exellent article. I am doing a project which requires me to understand the following concepts:
    tiled-poisson disk sampling
    variable density poisson disk sampling
    variable density tiled-poisson disk sampling
    Would you be able to kindly assisst me?

  • Pingback: Illyriad at Night | Illyriad - Beneath the Misted Land

  • Pingback: Tech Talk: The Voronoi Diagram - Arcade

  • Rhodri Cusack

    Thank you very much for the informative article. I’ve implemented your 2d and 3d algorithms in matlab, should that be useful to anyone! They’re linked from here http://cusacklabcomputing.blogspot.ca/2013/07/poisson-disc-2d-in-matlab.html

  • ymg

    Thanks, for an informative article. Since the grid normally as the same spacing as your min dist for points, Why do we need to check a 5×5 matrix. Seems to me that a 3×3 would be sufficient.

    • Herman Tulleken

      The image with the squares and circles should explain explain it. The circles show where we need to check (actually, the convex bounding curve of those four circles). They overlap slightly with 12 of the squares in the outer border of the 5×5 grid. So if we only checked the 3×3 grid, we would have the posibility of some points being too close. For example, if the red point was close to the top left corner of the center cell, and a blue point in the left column (of the 5×5 grid) close to the top left corner of the cell, it would be too close. If we only checked the 3×3 grid, we would miss it.

  • Pingback: Day 29 Game 29: 30 games in 30 days using Grids | Gamelogic

  • Pingback: Fun with 3D Time-of-Flight Cameras | Larrylisky's Wiki

  • Pingback: Fun with 3D Time-of-Flight Cameras | Larrylisky's Wiki

  • Pingback: Poisson Disk Sampling Example Code | code-spot

  • some guy

    boy, oh boy. i have just tried to implement this collection of poorly written and inaccurate code and the incompleteness and errors present in it gave me a 2 hours headache. why won’t you just plain copypate the original java code from the tokyo university instead? i ended up abandoning this mess on this blog in favor of the original one. i dare you, do try to implement the code from scratch, following your code. you will be ashamed.

    • http://www.gamelogic.co.za/ Herman Tulleken

      Since I wrote this 6 years ago I have definitely learned my lesson! Publishing pseudo-code is not as helpful (precisely because of the issues you mentioned: incompleteness and a tendency for errors to crop in). The Java, Ruby and Python code I wrote is all available from the link at the bottom, but it’s maybe a good idea to replace the pseudo-code in the text.

      (Of course, some have been able to overcome my messiness, so I am not too ashamed ;-))

      • Sei

        I just implemented this code from scratch, following this code, in Java.
        It took me about 2 hours including debugging my own stupid mistakes, and coding up the displaying of the distribution.

        Aside from minor difficulties parsing the pseudocode’s intent and translating that into Java (not too hard), my only problems were my own stupidity (mixing up variables, mostly).

        I’m not a very good programmer. I’m an amateur who messes around for a week or two every few months. The pseudocode was less than ideal, but it was perfectly fine, and lacked any discernable errors.

  • Pingback: Poisson disk points generator | Mobile Development

  • Sergey Kosarevsky
  • Scott Mitchell

    That’s great that a C++ implementation of Bridson’s method is provided by Sergey. I think that my more recent method is simple, and gets one truly to maximality. Sorry I don’t have permissions to share the code yet, but I hope it is simple enough to reimplement. http://www.cs.sandia.gov/~samitch/papers/eurographics_mps-final-with-appendix.pdf