mathr / blog / #

Here is an algorithm for generating circle packings:

1. start with an empty image
2. while the image has gaps bigger than a pixel
   1. pick a random unoccupied point
   2. draw the largest circle centered on that point that doesn't overlap the previous circles or the image boundary

This probably has been rediscovered many times, but I don't have a reference. To pack shapes other than circles, find the largest circle and put the shape inside it, oriented to touch the tangent with the image, or at random, as you wish.

An alternative algorithm that picks the sizes up front and fits them into the image has nicer fractal properties, but also has halting problems. See:

An Algorithm for Random Fractal Filling of Space

John Shier and Paul Bourke

Computational experiments with a simple algorithm show that it is possible to fill any spatial region with a random fractalization of any shape, with a continuous range of pre-specified fractal dimensions D. The algorithm is presented here in 1, 2, or 3 physical dimensions. The size power- law exponent c or the fractal dimension D can be specified ab initio over a substantial range. The method creates an infinite set of shapes whose areas (lengths, volumes) obey a power law and sum to the area (length, volume) to be filled. The algorithm begins by randomly placing the largest shape and continues using random search to place each smaller shape where it does not overlap or touch any previously placed shape. The resulting gasket is a single connected object.

I implemented both (using the GNU Scientific Library function for Hurwitz Zeta as my naive summation was woefully inaccurate) but the halting problems of the paper are very annoying, even though it produces nicer images. The problem with both is how to make it fast - the bottle-neck I found is the "pick a random unoccupied point" step.

Image pyramids, aka mipmaps, consist of a collection of downscaled-by-2 images reduced from a base layer. By performing high quality low pass filtering at each reduction level, the resulting collection can be used for realtime texturing of 3D objects of varying sizes without aliasing.

Histogram pyramids are similar, though instead of containing image data, each cell contains a count of the number of active cells in the corresponding base layer region. Thus the smallest 1x1 histogram layer contains the total number of active cells in the base layer. It's similar to a quad tree, but without all the cache-busting pointer chasing.

There are two operations of interest: the first is "deactivate a cell in the base layer", which can be done by decrementing all the cells in the path through the layers from the base layer to the smallest 1x1 layer. The complexity of this operation is O(layers) = O(log(max{width, height})). For a flat image without a histogram pyramid it would be O(1).

The second thing we want to do is "pick an active cell uniformly at random", which is where the histogram pyramid comes into its own. Now the technique picks a random subcell of each cell, starting from the 1x1 layer and proceeding towards the base layer. The random numbers are weighted according to the totals stored in the histogram pyramid, which ensures uniformity. Assuming the pseudo-random number generator is O(1) (which is very likely), this algorithm is again O(layers). Without a histogram pyramid, the best alternatives I came up with were O(number of active cells), which is O(width * height) at the start of the packing algorithm, thus the histogram pyramid is a huge improvement (100mins vs 10secs in one small test), or O(number of previously drawn circles), which is O(width * height) at the end of the packing algorithm, again very poor.

Some code:

/*
gcc -std=c99 -Wall -Wextra -pedantic -O3 -o x x.c -lm
./x W H > x.pgm
*/

#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

struct histogram
{
  int levels;
  int **counts;
};

struct histogram *h_new(int width, int height)
{
  int levels = 0;
  int size = 1 << levels;
  while (size < width || size < height)
  {
    levels += 1;
    size <<= 1;
  }
  struct histogram *h = malloc(sizeof(struct histogram));
  h->levels = levels;
  h->counts = malloc(sizeof(int *) * (levels + 1));
  for (int l = 0; l <= levels; ++l)
  {
    int d = 1 << l;
    int n = d * d;
    h->counts[l] = malloc(sizeof(int) * n);
  }
  for (int y = 0; y < size; ++y)
    for (int x = 0; x < size; ++x)
      h->counts[levels][(y << levels) + x] = y < height && x < width;
  for (int l = levels - 1; l >= 0; --l)
    for (int y = 0; y < 1 << l; ++y)
      for (int x = 0; x < 1 << l; ++x)
        h->counts[l][(y << l) + x]
          = h->counts[l+1][(((y<<1) + 0) << (l + 1)) + (x << 1) + 0]
          + h->counts[l+1][(((y<<1) + 0) << (l + 1)) + (x << 1) + 1]
          + h->counts[l+1][(((y<<1) + 1) << (l + 1)) + (x << 1) + 0]
          + h->counts[l+1][(((y<<1) + 1) << (l + 1)) + (x << 1) + 1];
  assert(h->counts[0][0] == width * height);
  return h;
}

void h_free(struct histogram *h)
{
  for (int l = 0; l <= h->levels; ++l)
    free(h->counts[l]);
  free(h->counts);
  free(h);
}

void h_decrement(struct histogram *h, int x, int y)
{
  for (int l = h->levels; l >= 0; --l)
  {
    int k = (y << l) + x;
    h->counts[l][k] -= 1;
    x >>= 1;
    y >>= 1;
  }
}

int h_empty(struct histogram *h)
{
  return h->counts[0][0] == 0;
}

int h_choose(struct histogram *h, int *x, int *y)
{
  if (h_empty(h)) return 0;
  *x = 0;
  *y = 0;
  for (int l = 1; l <= h->levels; ++l)
  {
    *x <<= 1;
    *y <<= 1;
    int xs[4] = { *x, *x + 1, *x, *x + 1 };
    int ys[4] = { *y, *y, *y + 1, *y + 1 };
    int ks[4] =
      { (ys[0] << l) + xs[0]
      , (ys[1] << l) + xs[1]
      , (ys[2] << l) + xs[2]
      , (ys[3] << l) + xs[3]
      };
    int ss[4] =
      { h->counts[l][ks[0]]
      , h->counts[l][ks[1]]
      , h->counts[l][ks[2]]
      , h->counts[l][ks[3]]
      };
    int ts[4] =
      { ss[0]
      , ss[0] + ss[1]
      , ss[0] + ss[1] + ss[2]
      , ss[0] + ss[1] + ss[2] + ss[3]
      };
    int p = rand() % ts[3];
    int i;
    for (i = 0; i < 4; ++i)
      if (p < ts[i]) break;
    *x = xs[i];
    *y = ys[i];
  }
  return 1;
}

struct delta
{
  int dx;
  int dy;
};

int cmp_delta(const void *a, const void *b)
{
  const struct delta *p = a;
  const struct delta *q = b;
  double x = p->dx * p->dx + p->dy * p->dy;
  double y = q->dx * q->dx + q->dy * q->dy;
  if (x < y) return -1;
  if (x > y) return  1;
  return 0;
}

struct delta *deltas;
int *image;
unsigned char *pgm;

struct histogram *initialize(int width, int height)
{
  deltas = malloc(sizeof(struct delta) * width * height);
  image = malloc(sizeof(int) * width * height);
  pgm = malloc(sizeof(unsigned char) * width * height);
  int d = 0;
  for (int dx = 0; dx < width; ++dx)
    for (int dy = 0; dy < height; ++dy)
    {
      deltas[d].dx = dx;
      deltas[d].dy = dy;
      ++d;
    }
  qsort(deltas, width * height, sizeof(struct delta), cmp_delta);
  for (int k = 0; k < width * height; ++k)
    image[k] = 0;
  return h_new(width, height);
}

void draw_circle(struct histogram *h, int width, int height,
  double cx, double cy, double r, int c)
{
  double r2 = r * r;
  int x0 = floor(cx - r);
  int x1 = ceil (cx + r);
  int y0 = floor(cy - r);
  int y1 = ceil (cy + r);
  for (int y = y0; y <= y1; ++y)
    if (0 <= y && y < height)
      for (int x = x0; x <= x1; ++x)
        if (0 <= x && x < width)
        {
          double dx = x - cx;
          double dy = y - cy;
          double d2 = dx * dx + dy * dy;
          if (d2 <= r2)
          {
            h_decrement(h, x, y);
            image[y * width + x] = c;
          }
        }
}

double find_radius(int width, int height, int cx, int cy)
{
  for (int d = 0; d < width * height; ++d)
  {
    int dx = deltas[d].dx;
    int dy = deltas[d].dy;
    int r2 = dx * dx + dy * dy;
    int xs[2] = { cx - dx, cx + dx };
    int ys[2] = { cy - dy, cy + dy };
    for (int j = 0; j < 2; ++j)
    {
      int y = ys[j];
      if (y < 0) return sqrt(r2);
      if (y >= height) return sqrt(r2);
      for (int i = 0; i < 2; ++i)
      {
        int x = xs[i];
        if (x < 0) return sqrt(r2);
        if (x >= width) return sqrt(r2);
        if (image[y * width + x]) return sqrt(r2);
      }
    }
  }
  return 0;
}

void packing(struct histogram *h, int width, int height, double p)
{
  for (int c = 1; h->counts[0][0] > 0; ++c)
  {
    fprintf(stderr, "%16d / %d\r", h->counts[0][0], width * height);
    int cx = 0;
    int cy = 0;
    h_choose(h, &cx, &cy);
    double r = find_radius(width, height, cx, cy);
    draw_circle(h, width, height, cx, cy, r, c);
  }
}

void edges(int width, int height)
{
  for (int y = 0; y < height; ++y)
  {
    for (int x = 0; x < width; ++x)
    {
      if (x == 0) pgm[y * width + x] = 0;
      else if (x == width-1) pgm[y * width + x] = 0;
      else if (y == 0) pgm[y * width + x] = 0;
      else if (y == height-1) pgm[y * width + x] = 0;
      else
      {
        int e = (image[y * width + x] != image[(y+1) * width + x+1])
             || (image[(y+1) * width + x] != image[y * width + x+1]);
        pgm[y * width + x] = e ? 0 : 255;
      }
    }
  }
}

void write_pgm(unsigned char *p, int width, int height)
{
  fprintf(stdout, "P5\n%d %d\n255\n", width, height);
  fwrite(p, width * height, 1, stdout);
  fflush(stdout);
}

int main(int argc, char **argv)
{
  if (argc < 3) return 0;
  srand(time(0));
  int width = atoi(argv[1]);
  int height = atoi(argv[2]);
  struct histogram *h = initialize(width, height);
  packing(h, width, height);
  edges(width, height);
  write_pgm(pgm, width, height);
  return 0;
}