Mandelbrot image in R

mandelbrot.R

library(ggplot2)
library(grid)
library(doParallel)
library(itertools)

resolution = 1000
iterations = 100
threshold = 10

# image range
data = expand.grid(
  a = seq(-.65, -.5, length.out=resolution), 
  b = seq(.55, .70, length.out=resolution))

# initialize
data$ab = complex(real=data$a, imaginary=data$b)
data$d = iterations

# functions
distance = function(x, y) sqrt(x^2 + y^2)
mandelbrot = function(z, c) z^2 + c

# compute
compute = function (ab) {
  z = ab
  d = iterations # number of iterations until the point reaches a threashold distance
  for (i in 1:iterations) {
    z = mandelbrot(z, ab)
    real = Re(z)
    imaginary = Im(z)
    dist = distance(real, imaginary)
    d = ifelse(dist>threshold & d>i, i, d)
  }
  return(d)
}

# chunked parallel compute
computeChunk = function(x) {
  # Create a cluster with 1 fewer cores than are available. Adjust as necessary
  cl = makeCluster(detectCores() - 1) # 8-1 cores 
  registerDoParallel(cl)
  # foreach
  out = foreach(xc=isplitVector(x, chunks=getDoParWorkers()),
          .export=c('mandelbrot', 'distance', 'iterations', 'threshold', 'compute'), 
          .combine='c') %dopar% {
            compute(xc)
          }
  # finish
  stopCluster(cl)
  return(out)
}

# single threaded
#system.time(data$d <- compute(data$ab))
# multi-threaded
system.time(data$d <- computeChunk(data$ab))

# color gradient
colors = c('#000000', '#000000', '#cccc22', '#000000', '#ff0000', '#000000', '#112266', '#aaaa22', '#ffffff')
values = c(0,         0.1,       0.49,      .491,      .51,       .52,       .53,       .99,       1)

plot = ggplot(data) + 
  geom_raster(aes(x=a, y=b, fill=log(d))) +
  scale_fill_gradientn(colours=colors, values=values, guide='none') +
  scale_x_continuous(expand = c(0,0)) + 
  scale_y_continuous(expand = c(0,0))

# strip out everything but the plot
gt <- ggplot_gtable(ggplot_build(plot))
ge <- subset(gt$layout, name == "panel")
grid.draw(gt[ge$t:ge$b, ge$l:ge$r])

# save image
png('mandelbrot.png', width=resolution, height=resolution)
grid.draw(gt[ge$t:ge$b, ge$l:ge$r])
dev.off()

mandelbrot2k.png

mandelbrot4k.png