[R] optical art

optical.art.R

optical.art <- function(P, N, colors = brewer.pal(P, "Greys")) {

## This function creates an Op art figure as can be found on
## http://r-de-jeu.blogspot.com/
##
## Inputs:
##    - P: (integer) number of starting points
##    - N: (integer) number of semi-lines starting from each point.
##                   For best results, N should be a multiple of P 
##    - colors:      vector of colors for filling polygons.
##                   For best results, you should use P colors

   require(sp)
   require(rgeos)
   require(RColorBrewer)

   # plot area
   plot(c(0, 1), c(0, 1), "n", axes=FALSE, xlab="", ylab="")

   create.polygon <- function(x, y) {
      pol <- Polygon(cbind(c(x, head(x, 1)), c(y, head(y, 1))))
      set <- Polygons(list(pol), "pol")
      spP <- SpatialPolygons(list(set))
   }

   plot.area <- create.polygon(x = c(0,1,1,0), y = c(0,0,1,1))

   polygons <- vector("list", length = P)

   for (p in 1:P) {

      # randomly select the coordinates of a point and the
      # angles of N semi-lines starting from that point
      center.x <- runif(1)
      center.y <- runif(1)
      angles   <- 2*pi/N*(0:(N-1) + 2*pi*runif(1))

      # the semi-lines and frame will define N polygons
      polygons[[p]] <- vector("list", length = N)
      for (i in 1:N) {

         i1 <- i
         i2 <- ifelse((i + 1) > N, 1, i + 1)
         a1 <- angles[i1]
         a2 <- angles[i2]

         polygons[[p]][[i]] <-
            create.polygon(x = center.x + c(0, 2*cos(a1), 2*cos(a2), 0),
                           y = center.y + c(0, 2*sin(a1), 2*sin(a2), 0))
      }
   }

   all.comb <- expand.grid(data.frame(replicate(P, 1:N)))
   for (k in 1:nrow(all.comb)) {

      # compute the intersection of the polygons
      poly <- plot.area
      for (i in 1:P) {
         poly <- gIntersection(poly, polygons[[i]][[all.comb[k,i]]])
         if (gIsEmpty(poly)) break
      }
      if (gIsEmpty(poly)) next

      # plot the polygon and fill it with the appropriate color
      coords <- poly@polygons[[1]]@Polygons[[1]]@coords
      polygon(coords[, 1], coords[, 2], border = NA,
              col = colors[1 + (sum(all.comb[k,]) %% length(colors))])
   }
}