從 chaos game 模擬分形 – Barnsley fern

上次透過 affine transformation 以及 Iterative Function System 來模擬Sierpinski triangle。

這次則是模擬 Barnsley fern

我們隨機選取一點初始點 p0 = (x0, y0) 後，隨機地將 p0 代入以下四個 affine transformation 其中之一:

$$ A:\begin{pmatrix} x_1\\ y_1 \end{pmatrix} = \begin{bmatrix} 0.85 & 0.04\\ -0.04 & 0.85 \end{bmatrix} \begin{pmatrix} x_0\\ y_0 \end{pmatrix} + \begin{pmatrix} 0.08\\ 0.18 \end{pmatrix} $$

$$ B:\begin{pmatrix} x_1\\ y_1 \end{pmatrix} = \begin{bmatrix} 0.2 & 0.23\\ -0.26 & 0.22 \end{bmatrix} \begin{pmatrix} x_0\\ y_0 \end{pmatrix} + \begin{pmatrix} 0.4\\ 0.05 \end{pmatrix} $$

$$ C: \begin{pmatrix} x_1\\ y_1 \end{pmatrix} = \begin{bmatrix} 0 & 0\\ 0 & 0.16 \end{bmatrix} \begin{pmatrix} x_0\\ y_0 \end{pmatrix} + \begin{pmatrix} 0.5\\ 0 \end{pmatrix} $$

$$ D: \begin{pmatrix} x_1\\ y_1 \end{pmatrix} = \begin{bmatrix} -0.15 & 0.26 \\ 0.28 & 0.14 \end{bmatrix} \begin{pmatrix} x_0\\ y_0 \end{pmatrix} + \begin{pmatrix} 0.58\\ -0.09 \end{pmatrix} $$

得到 p1 後再次隨機代入 A, B, C, D 其中一個 affine transformation ，如此往復，將所有的點繪製在平面座標上就可以得到：

R語言程式碼如下：

library(magrittr)
library(dplyr)
library(ggplot2)

# 分形 fractal
# Barnsley's fern

## affine transformations
matA <- matrix(c(0.85, -0.04, 0.04, 0.85), 
               nrow = 2, ncol = 2)
vecA0 <- matrix(c(0.08, 0.18))
matB <- matrix(c(0.2, 0.23, -0.26, 0.22), 
               nrow = 2, ncol = 2)
vecB0 <- matrix(c(0.4, 0.05))
matC <- matrix(c(0, -0, 0, 0.16), 
               nrow = 2, ncol = 2)
vecC0 <- matrix(c(0.5, 0))
matD <- matrix(c(-0.15, 0.26, 0.28, 0.14), 
               nrow = 2, ncol = 2)
vecD0 <- matrix(c(0.58, -0.09))

p0 <- matrix(runif(2), nrow = 2)

affineA <- function(p0){return(matA%*%p0 + vecA0)}
affineB <- function(p0){return(matB%*%p0 + vecB0)}
affineC <- function(p0){return(matC%*%p0 + vecC0)}
affineD <- function(p0){return(matD%*%p0 + vecD0)}

diceRoll <- function(){
  return(sample(c("A", "B", "C", "D"), size = 1, prob = rep(1/4, 4)))
}

affineTrans <- function(diceOutCome, initP0){
  if(diceOutCome == "A"){
    affineA(initP0)
  } else if(diceOutCome == "B"){
    affineB(initP0)
  } else if(diceOutCome == "C"){
    affineC(initP0)
  } else if(diceOutCome == "D"){
    affineD(initP0)
  }
}

# keep track of path
run <- function(n_runs){
  p0 <- matrix(runif(2), nrow = 2)
  df <- t(p0)
  for(i in 1:n_runs){
    p0 <- t(tail(df, 1))
    p1 <- affineTrans(diceRoll(), p0)
    df <- rbind(df, t(p1))
    rownames(df) <- NULL

    if( (t(p1)%*%p1) >= 10000){break}
  }
  colnames(df) <- c("X", "Y")
  return(df)
}
df <- run(10000)
plot(df, pch = 16)

# using ggplot2
df %>% 
  as_tibble() %>% 
  ggplot(aes(X, Y)) +
  geom_point()