Circle packing with R

July 26, 2010
By

(This article was first published on Last Resort Software, and kindly contributed to R-bloggers)

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)
}

To leave a comment for the author, please follow the link and comment on his blog: Last Resort Software.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.