Circle packing with R

To visualize the results of a simulation model of woodland trees within R, I needed an algorithm that could arrange a large number of circles within a rectangle such that no two circles overlapped by more than a specified amount. A colleague had approached this problem earlier by sorting the circles in order of descending size, then randomly dropping each one into the rectangle repeatedly until it landed in a position with acceptable overlap.

I suspected a faster and more robust algorithm could be constructed using some kind of "jiggling the circles" approach. Luckily for me, I discovered that Sean McCullough had written a really nice example of circles packing into a cluster using the Processing language. Sean's program is based on an iterative pair-repulsion algorithm in which overlapping circles move away from each other. Based on this, and modifying the algorithm a little, I came up with an R function to produce constrained random layouts of a given set of circles. Here are two layouts for the same input set...

The algorithm is wonderfully simple. Each circle in the input list is compared to those following it. If two circles overlap by more than the allowed amount, they are moved apart until the overlap is acceptable. The distance moved by each circle is proportional to the radius of the other to approximate inertia (very loosely); thus when a small circle is overlapped by a large circle, the small circle moves furthest. This process is repeated iteratively until no more movement takes place (acceptable layout) or a maximum number of iterations is reached (layout failure). To avoid edge effects, the bounding rectangle is treated as a toroid. For my purposes, I only require a circle's centre to be within the plotted rectangle.

You can see in the plots that the algorithm tends to produce clusters of small circles around big ones. For the woodland simulation model this is a nice property (saplings clustering around parent trees) but for other applications the algorithm could be further modified to lessen or avoid this effect.

The code for the basic function is below. The set of input circles are described by the config matrix argument. Although this function produces good results, it takes a long time to run when the number of circles is large. However, a second version of the function, with most of the calculations moved into C code, runs much faster (code not posted here but available on request).

pack.circles <- function(config, size=c(100, 100), max.iter=1000, overlap=0 ) {
  #
  # Simple circle packing algorithm based on inverse size weighted repulsion
  #
  # config   - matrix with two cols: radius, N
  # size     - width and height of bounding rectangle
  # max.iter - maximum number of iterations to try
  # overlap  - allowable overlap expressed as proportion of joint radii

  # ============================================================================
  #  Global constants
  # ============================================================================
  # round-off tolerance
  TOL <- 0.0001

  # convert overlap to proportion of radius
  if (overlap < 0 | overlap >= 1) {
    stop("overlap should be in the range [0, 1)")
  }
  PRADIUS <- 1 - overlap

  NCIRCLES <- sum(config[,2])

  # ============================================================================
  #  Helper function - Draw a circle
  # ============================================================================
  draw.circle <- function(x, y, r, col) { 
    lines( cos(seq(0, 2*pi, pi/180)) * r + x, sin(seq(0, 2*pi, pi/180)) * r + y , col=col )
  }


  # ============================================================================
  #  Helper function - Move two circles apart. The proportion of the required
  #  distance moved by each circle is proportional to the size of the other 
  #  circle. For example, If a c1 with radius r1 overlaps c2 with radius r2,
  #  and the movement distance required to separate them is ds, then c1 will
  #  move ds * r2 / (r1 + r2) while c2 will move ds * r1 / (r1 + r2). Thus,
  #  when a big circle overlaps a little one, the little one moves a lot while
  #  the big one moves a little.
  # ============================================================================
  repel <- function(xyr, c0, c1) {
    dx <- xyr[c1, 1] - xyr[c0, 1]
    dy <- xyr[c1, 2] - xyr[c0, 2]
    d <- sqrt(dx*dx + dy*dy)
    r <- xyr[c1, 3] + xyr[c0, 3]
    w0 <- xyr[c1, 3] / r
    w1 <- xyr[c0, 3] / r

    if (d < r - TOL) {
      p <- (r - d) / d
      xyr[c1, 1] <<- toroid(xyr[c1, 1] + p*dx*w1, 1)
      xyr[c1, 2] <<- toroid(xyr[c1, 2] + p*dy*w1, 2)
      xyr[c0, 1] <<- toroid(xyr[c0, 1] - p*dx*w0, 1)
      xyr[c0, 2] <<- toroid(xyr[c0, 2] - p*dy*w0, 2)

      return(TRUE)
    }

    return(FALSE)
  }


  # ============================================================================
  #  Helper function - Adjust a coordinate such that if it is distance d beyond
  #  an edge (ie. outside the area) it is moved to be distance d inside the 
  #  opposite edge. This has the effect of treating the area as a toroid.
  # ============================================================================
  toroid <- function(coord, axis) {
    tcoord <- coord

    if (coord < 0) {
      tcoord <- coord + size[axis]
    } else if (coord >= size[axis]) {
      tcoord <- coord - size[axis]
    }

    tcoord
  }


  # ============================================================================
  #  Main program
  # ============================================================================

  # ------------------------------------------
  # create a random initial layout
  # ------------------------------------------
  xyr <- matrix(0, NCIRCLES, 3)

  pos0 <- 1
  for (i in 1:nrow(config)) {
    pos1 <- pos0 + config[i,2] - 1
    xyr[pos0:pos1, 1] <- runif(config[i, 2], 0, size[1])
    xyr[pos0:pos1, 2] <- runif(config[i, 2], 0, size[2])
    xyr[pos0:pos1, 3] <- config[i, 1] * PRADIUS
    pos0 <- pos1 + 1
  }

  # ------------------------------------------
  # iteratively adjust the layout
  # ------------------------------------------
  for (iter in 1:max.iter) {
    moved <- FALSE
    for (i in 1:(NCIRCLES-1)) {
      for (j in (i+1):NCIRCLES) {
        if (repel(xyr, i, j)) {
          moved <- TRUE
        }
      }
    }
    if (!moved) break
  }

  cat(paste(iter, "iterations\n"));
  
  # ------------------------------------------
  # draw the layout
  # ------------------------------------------
  plot(0, type="n", xlab="x", xlim=c(0,size[1]), ylab="y", ylim=c(0, size[2]))

  xyr[, 3] <- xyr[, 3] / PRADIUS
  for (i in 1:nrow(xyr)) {
    draw.circle(xyr[i, 1], xyr[i, 2], xyr[i, 3], "gray")
  }

  # ------------------------------------------
  # return the layout
  # ------------------------------------------
  colnames(xyr) <- c("x", "y", "radius")
  invisible(xyr) 
}