Draw a polygon colored like this in R or Matlab

Question

http://www.texample.net/tikz/examples/lindenmayer-systems/

My sample code shown below, I don't know how to colored with hue color.

plot.koch <- function(k,col="blue"){ 
  plot(0,0,xlim=c(0,1), ylim=c(-sqrt(3)/6,sqrt(3)/2), asp = 1,type="n",xlab="", ylab="")
  plotkoch <- function(x1,y1,x2,y2,n){
    if (n > 1){
      plotkoch(x1,y1,(2*x1+x2)/3,(2*y1+y2)/3,n-1); 
      plotkoch((2*x1+x2)/3,(2*y1+y2)/3,(x1+x2)/2-(y1-y2)*sqrt(3)/6,(y1+y2)/2-(x2-x1) *sqrt(3)/6,n-1);
      plotkoch((x1+x2)/2-(y1-y2)*sqrt(3)/6,(y1+y2)/2-(x2-x1)*sqrt(3)/6,(2*x2+x1)/3,(2 *y2+y1)/3,n-1); 
      plotkoch((2*x2+x1)/3,(2*y2+y1)/3,x2,y2,n-1)
    }    
    else { 
      x=c(x1,(2*x1+x2)/3,(x1+x2)/2-(y1-y2)*sqrt(3)/6,(2*x2+x1)/3,x2); 
      y=c(y1,(2*y1+y2)/3,(y1+y2)/2-(x2-x1)*sqrt(3)/6,(2*y2+y1)/3,y2); 
      polygon(x,y,type="l",col=col) 
    }
  }
  plotkoch(0,0,1,0,k) 
  plotkoch(0.5,sqrt(3)/2,0,0,k) 
  plotkoch(1,0,0.5,sqrt(3)/2,k)
}
  plot.koch(3, col=3)  

Answer 1

Here's a method using spatial objects in R, with sp, rgeos and raster packages in the mix.

1. Slight modifications to the function to return the x,y coordinates to the user (and in the correct order):

   koch <- function(k) { 
     yy <- xx <- numeric(0)
     Koch <- function(x1, y1, x2, y2, n) {
       if (n > 1){
         Koch(x1, y1, (2*x1+x2)/3, (2*y1+y2)/3, n-1); 
         Koch((2*x1+x2)/3, (2*y1+y2)/3, (x1+x2)/2-(y1-y2)*sqrt(3)/6, 
              (y1+y2)/2-(x2-x1) *sqrt(3)/6, n-1);
         Koch((x1+x2)/2-(y1-y2)*sqrt(3)/6, (y1+y2)/2-(x2-x1)*sqrt(3)/6, 
              (2*x2+x1)/3, (2 *y2+y1)/3, n-1); 
         Koch((2*x2+x1)/3, (2*y2+y1)/3, x2, y2, n-1)
       }    
       else { 
         x <- c(x1, (2*x1+x2)/3, (x1+x2)/2-(y1-y2)*sqrt(3)/6, (2*x2+x1)/3, x2); 
         xx <<- c(xx, x)
         y <- c(y1, (2*y1+y2)/3, (y1+y2)/2-(x2-x1)*sqrt(3)/6, (2*y2+y1)/3, y2); 
         yy <<- c(yy, y)
        }
     }
     Koch(0, 0, 1, 0, k)
     Koch(1, 0, 0.5, sqrt(3)/2, k)
     Koch(0.5, sqrt(3)/2, 0, 0, k) 
     xy <- data.frame(x=xx, y=yy)
     rbind(unique(xy), xy[1, ])
   }
   

2. Create a colour ramp:

   colr <- colorRampPalette(hcl(h=seq(0, 360, len=100), c=100))
   

3. Use koch function to get vertices:

   xy <- koch(4)
   

4. Load spatial packages and create SpatialPolygons object from fractal and plot it once to set up the plot area.

   library(sp)
   library(rgeos)
   library(raster)
   
   poly <- SpatialPolygons(list(Polygons(list(Polygon(xy)), 1)))
   plot(poly)
   

5. Plot a series of segments with desired origin and large enough radius to cover the fractal polygon (here we use radius r <- 1).

   r <- 1
   mapply(function(theta, col) {
     segments(0.5, 0.3, 0.5 + r*cos(theta), 0.3 + r*sin(theta), lwd=3, col=col) 
   }, seq(0, 360, length=1000)*pi/180, colr(1000))       
   

6. Create a second polygon of the difference between the plot area and the fractal polygon, and plot this (with col='white') to mask out the unwanted gradient area.

   plot(gDifference(as(extent(par('usr')), 'SpatialPolygons'), poly), 
        col='white', border='white', add=TRUE)
   

7. Plot the polygon once more.

   plot(poly, add=TRUE)
   

Answer 2

Here's my attempt at solving your question. Currently it draws the color also outside of the snowflake. If you can figure out if points are inside or outside the snowflake, you should be able to just remove outside points in the df_fill. Here I'm first creating the data.frame used for plotting the polygon. Then I'm creating the data.frame for the background color. And finally I'm using ggplot2 to plot the data.

# creating relevant data
data.koch <- function(k){ 
  df <- data.frame(x = 0, 
                   y = 0, 
                   grp = 0)
  plotkoch <- function(x1, y1, x2, y2, n, data){
    if (n==1) {
      x=c(x1,(2*x1+x2)/3,(x1+x2)/2-(y1-y2)*sqrt(3)/6,(2*x2+x1)/3,x2) 
      y=c(y1,(2*y1+y2)/3,(y1+y2)/2-(x2-x1)*sqrt(3)/6,(2*y2+y1)/3,y2) 
      df <- rbind(data, data.frame(x, y, grp=max(data$grp)+1))
    }
    if (n > 1){
      df <- plotkoch(x1,y1,(2*x1+x2)/3,(2*y1+y2)/3,n-1, data = data)
      df <- plotkoch((2*x1+x2)/3,(2*y1+y2)/3,(x1+x2)/2-(y1-y2)*sqrt(3)/6,(y1+y2)/2-(x2-x1) *sqrt(3)/6,n-1, data=df)
      df <- plotkoch((x1+x2)/2-(y1-y2)*sqrt(3)/6,(y1+y2)/2-(x2-x1)*sqrt(3)/6,(2*x2+x1)/3,(2 *y2+y1)/3,n-1, data=df) 
      df <- plotkoch((2*x2+x1)/3,(2*y2+y1)/3,x2,y2,n-1, data=df)
    }    
    return(df)
  }
  df <- plotkoch(0,0,1,0,k, data = df) 
  df <- plotkoch(0.5,sqrt(3)/2,0,0,k, data = df) 
  df <- plotkoch(1,0,0.5,sqrt(3)/2,k, data = df)
  return(df)
}
# plotting functon
plot.koch <- function(k){
  stopifnot(require(ggplot2))
  if (is.data.frame(k)) df <- k
  else df <- data.koch(k)
  # filling data (CHANGE HERE TO GET ONLY INSIDE POINTS)
  l <- 500
  df_fill <- expand.grid(x=seq(0, 1, length=l), 
                         y=seq(-sqrt(3)/6, sqrt(3)/2, length=l))
  df_fill[, "z"] <- atan2(-df_fill[, "y"] + sqrt(3)/6, df_fill[, "x"] - 0.5) + pi/2
  df_fill[df_fill[, "z"] < 0, "z"] <- df_fill[df_fill[, "z"] < 0, "z"] + 2*pi
  # plotting
  ggplot(df, aes(x, y, group=grp)) + 
    geom_raster(data = df_fill, 
                aes(fill=z, group=NULL), 
                hjust = 0,
                vjust = 0,
                linetype='blank') + 
    geom_path(data=df, size=1) +
    scale_fill_gradientn(colours = rainbow(30), guide = 'none') + 
    scale_x_continuous(name = '', limits = c(0, 1), expand=c(0, 0)) + 
    scale_y_continuous(name = '', limits = c(-sqrt(3)/6,sqrt(3)/2), expand=c(0, 0)) +
    coord_fixed() +
    theme_bw() +
    theme(axis.line = element_blank(), 
          panel.grid = element_blank(), 
          axis.ticks = element_blank(), 
          axis.text = element_blank())
}
#
p <- plot.koch(4)
print(p)