Draw a polygon colored like this in R or Matlab

Question

http://www.texample.net/tikz/examples/lindenmayer-systems/ enter image description here

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)

score 3 · Accepted Answer · answered Feb 10 '15 at 13:07

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

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

Create a colour ramp:

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

Use koch function to get vertices:
```
xy <- koch(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)

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

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)
```
Plot the polygon once more.
```
plot(poly, add=TRUE)
```

enter image description here

score 1 · Answer 2 · edited May 23 '17 at 12:10

I would do it like this:

for any drawed pixel obtain its position x,y
compute the angle=atan2(y-y0,x-x0)

where x0,y0 is the koch's snowflake mid position
compute the color based on angle

if you use HSV then hue=angle and compute the target color value (I assume RGB). If you want the visible spectra colors you can try mine:
- RGB values of visible spectrum
Just convert the angle range angle=<0,2*Pi> [rad] to wavelength l=<400,700> [nm] so:
```
l = 400.0 + (700.0-400.0)*angle/(2.0*M_PI)
```
render the pixel

[Notes]

not using R nor Matlab so you need to code it yourself. The angle may need some shifting to match your coordinate system for example:

const angle0=???; // some shift constant [rad]
angle+=angle0; // or angle=angle0-angle; if the direction is oposite
if (angle>=2.0*M_PI) angle-=2.0*M_PI;
if (angle<      0.0) angle+=2.0*M_PI;

If you drawing this as polygon then you need to compute color per vertex not per pixel but then you can get to problems because this is not convex polygon. So how to ensure the mid point color ??? I am afraid you will need to use some sort of triangulation because simple triangle fan will fail ...

The only thing that is obvious is to fill the color for whole space and then draw the outline with black color then flood fill all non black pixels from outside with white color

score 1 · Answer 3 · answered Feb 10 '15 at 13:05

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)

output from plot.koch(4)

score 1 · Answer 4 · edited Mar 26 '17 at 20:09

It's my solution with grid package.

##data
  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 - min(xx))/(max(xx) - min(xx)), y = (yy - 
            min(yy))/(max(yy) - min(yy)))
        rbind(unique(xy), xy[1, ])
    }
xy <- koch(5)
##Plot
library(grid)
grid.newpage()
pushViewport(dataViewport(xy$x, xy$y), plotViewport(c(1, 1, 1, 1)))
    for (i in 1:nrow(xy)) {
       grid.path(x = c(xy[i, 1], xy[i + 1, 1], mean(xy$x)), 
                 y = c(xy[i, 2], xy[i + 1, 2], mean(xy$y)), 
                 gp = gpar(col = rainbow(nrow(xy))[i], 
                           fill = rainbow(nrow(xy))[i]))
        }

Here is the image

Draw a polygon colored like this in R or Matlab

4 Answers4

Linked