### 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) xyr[pos0:pos1, 2] <- runif(config[i, 2], 0, size) 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), ylab="y", ylim=c(0, size)) 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) }`

### Circle packing in R (again)

Back in 2010 I posted some R code for circle packing . Now, just five years later, I've ported the code to Rcpp and created a little package which you can find at GitHub . The main function is circleLayout which takes a set of overlapping circles and tries to find a non-overlapping arrangement for them. Here's an example: And here's the code: # Create some random circles, positioned within the central portion # of a bounding square, with smaller circles being more common than # larger ones. ncircles <- 200 limits <- c(-50, 50) inset <- diff(limits) / 3 rmax <- 20 xyr <- data.frame( x = runif(ncircles, min(limits) + inset, max(limits) - inset), y = runif(ncircles, min(limits) + inset, max(limits) - inset), r = rbeta(ncircles, 1, 10) * rmax) # Next, we use the `circleLayout` function to try to find a non-overlapping # arrangement, allowing the circles to occupy any part of the bounding square. # The returned value is a list with elements for

### Fitting an ellipse to point data

Some time ago I wrote an R function to fit an ellipse to point data, using an algorithm developed by Radim Halíř and Jan Flusser 1 in Matlab, and posted it to the r-help list . The implementation was a bit hacky, returning odd results for some data. A couple of days ago, an email arrived from John Minter asking for a pointer to the original code. I replied with a link and mentioned that I'd be interested to know if John made any improvements to the code. About ten minutes later, John emailed again with a much improved version ! Not only is it more reliable, but also more efficient. So with many thanks to John, here is the improved code: fit.ellipse <- function (x, y = NULL) { # from: # http://r.789695.n4.nabble.com/Fitting-a-half-ellipse-curve-tp2719037p2720560.html # # Least squares fitting of an ellipse to point data # using the algorithm described in: # Radim Halir & Jan Flusser. 1998. # Numerically stable direct least squares fitting of ellipses

### Graph-based circle packing

The previous two posts showed examples of a simple circle packing algorithm using the packcircles package (available from CRAN and GitHub ). The algorithm involved iterative pair-repulsion to jiggle the circles until (hopefully) a non-overlapping arrangement emerged. In this post we'll look an alternative approach. An algorithm to find an arrangement of circles satisfying a prior specification of circle sizes and tangencies was described by Collins and Stephenson in their 2003 paper in Computation Geometry Theory and Applications. A version of their algorithm was implemented in Python by David Eppstein as part of his PADS library (see CirclePack.py ). I've now ported David's version to R/Rcpp and included it in the packcircles package. In the figure below, the graph on the left represents the desired pattern of circle tangencies: e.g. circle 7 should touch all of, and only, circles 1, 2, 6 and 8. Circles 5, 7, 8 and 9 are internal , while the remaining circles are exter