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

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 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)
{
  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 * r2;
  //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);
}