彩色繪製這樣的多邊形R或Matlab的

Question

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

如下圖所示我的示例代碼中，我不知道怎麼用彩色色調的顏色。

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

我會做這樣的：

1. 任何drawed像素獲取其位置x,y

2. 計算angle=atan2(y-y0,x-x0)

   其中x0,y0是科赫雪花中旬位置

3. 計算基於角度

   如果使用HSV然後hue=angle和計算目標顏色值的顏色（我假設RGB）。如果你想在可見光譜的顏色，你可以試試我的：

   • RGB values of visible spectrum

   只是轉換的角度範圍angle=<0,2*Pi> [rad]到波長l=<400,700> [nm]這樣：

   l = 400.0 + (700.0-400.0)*angle/(2.0*M_PI) 
   

4. 渲染像素

[注意事項]

不使用[R也不Matlab的所以你需要自己編寫它。該角度可能需要一些轉移到你的座標系匹配，例如：

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; 

如果你畫這個多邊形，那麼你需要計算不是每個像素的每個頂點的顏色，但那麼你可以得到一些問題，因爲這不是凸多邊形。那麼如何確保中點顏色？恐怕你將需要使用某種三角形，因爲簡單的三角形風扇將失敗...

唯一明顯的是填充整個空間的顏色，然後用黑色繪製輪廓，然後填充所有來自外部的非黑色像素均爲白色

Answer 2

以下是我在解決您的問題時的嘗試。目前它也在雪花外繪製顏色。如果您可以確定點位於雪花內部還是外部，您應該能夠在df_fill中刪除外部點。 我在這裏首先創建用於繪製多邊形的data.frame。然後，我爲背景色創建data.frame。最後我使用ggplot2來繪製數據。

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

Answer 3

下面是使用在R個空間對象的方法，其中在混合sp，rgeos和raster包。

1. 稍微修改函數返回的x，y座標給用戶（和以正確的順序）：

   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. 創建色階：

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

3. 使用koch函數獲取頂點：

   xy <- koch(4) 
   

4. 加載空間包並從分形創建對象並將其繪製一次以設置繪圖區域。

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

5. 疊加了一系列所需的來源和足夠大的半徑段覆蓋多邊形分形（這裏我們使用半徑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. 創建繪圖區和分形多邊形之間的差的第二多邊形，並繪製這樣（col='white'），以屏蔽掉不需要的梯度區域。

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

7. 再次繪製多邊形。

   plot(poly, add=TRUE) 
   

Answer 4

這是我與電網一攬子解決方案。

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