Rでペンローズタイルを描く

この記事は R Advent Calendar 2014 の25日目の記事です。


ペンローズ・タイルって皆様ご存知でしょうか。
ペンローズ・タイル - Wikipedia


簡単にまとめると以下のような特徴を持つ不思議なタイルとその平面充填方法です。

  • 非周期的であり、同じ充填方法が現れない。
  • kiteとdirtと呼ばれる二つの四角形から構成される。
  • それぞれの四角形は同じ形のロビンソンの三角形という二等辺三角形を2つ組み合わせたもので出来ている。
  • 三角形の長辺との長さは黄金比になっている。
  • 再帰的に分割することができ、フラクタル図形の一種と言える。


このペンローズタイルと同じ性質を持つ図形はイスラム建築の模様や、電動シェーバーの網刃、準結晶等様々なところに見られます。
とっても不思議な図形ですね。まとまった解説記事がないのが残念です。


ということで、このペンローズタイルをRで描いてみました。せっかくなのでクリスマスツリー状にしています。

PHI<<-(1+sqrt(5))/2   
L<<-2

decompose_kite<-function(t1,k,rev=FALSE){
  if(rev){
    x1<-c((t1$coord$x[2]/PHI+t1$coord$x[1])/PHI,t1$coord$x[3],t1$coord$x[1])
    y1<-c((t1$coord$y[2]/PHI+t1$coord$y[1])/PHI,t1$coord$y[3],t1$coord$y[1])
    x2<-c((t1$coord$x[2]/PHI+t1$coord$x[1])/PHI,t1$coord$x[3],(t1$coord$x[2]+t1$coord$x[3]/PHI)/PHI)
    y2<-c((t1$coord$y[2]/PHI+t1$coord$y[1])/PHI,t1$coord$y[3],(t1$coord$y[2]+t1$coord$y[3]/PHI)/PHI)
    x3<-c(t1$coord$x[2],(t1$coord$x[2]+t1$coord$x[3]/PHI)/PHI,(t1$coord$x[2]/PHI+t1$coord$x[1])/PHI)
    y3<-c(t1$coord$y[2],(t1$coord$y[2]+t1$coord$y[3]/PHI)/PHI,(t1$coord$y[2]/PHI+t1$coord$y[1])/PHI)
  }else{
    x1<-c(t1$coord$x[3],t1$coord$x[1],(t1$coord$x[2]/PHI+t1$coord$x[3])/PHI)
    y1<-c(t1$coord$y[3],t1$coord$y[1],(t1$coord$y[2]/PHI+t1$coord$y[3])/PHI)
    x2<-c((t1$coord$x[2]/PHI+t1$coord$x[3])/PHI,t1$coord$x[1],(t1$coord$x[2]+t1$coord$x[1]/PHI)/PHI)
    y2<-c((t1$coord$y[2]/PHI+t1$coord$y[3])/PHI,t1$coord$y[1],(t1$coord$y[2]+t1$coord$y[1]/PHI)/PHI)
    x3<-c(t1$coord$x[2],(t1$coord$x[2]+t1$coord$x[1]/PHI)/PHI,(t1$coord$x[2]/PHI+t1$coord$x[3])/PHI)
    y3<-c(t1$coord$y[2],(t1$coord$y[2]+t1$coord$y[1]/PHI)/PHI,(t1$coord$y[2]/PHI+t1$coord$y[3])/PHI)
  }  
  small_coord <-list(x=x1,y=y1)
  small_coord2<-list(x=x2,y=y2)
  big_coord   <-list(x=x3,y=y3)
  
  polygon(small_coord$x, small_coord$y, col = "green")
  polygon(small_coord2$x, small_coord2$y, col = "green")
  polygon(big_coord$x, big_coord$y, col = "red")
  polygon(-small_coord$x, small_coord$y, col = "green")
  polygon(-small_coord2$x, small_coord2$y, col = "green")
  polygon(-big_coord$x, big_coord$y, col = "red")
  
  if(k<L && rev==FALSE){
    decompose_dart(  list(coord=big_coord),    k=k+1, rev)
    decompose_kite(list(coord=small_coord),  k=k+1,!rev)
    decompose_kite(list(coord=small_coord2), k=k+1, rev)
  }
  if(k<L && rev==TRUE){
    decompose_dart(  list(coord=big_coord),     k=k+1,rev)
    decompose_kite(list(coord=small_coord),   k=k+1,rev)
    decompose_kite(list(coord=small_coord2),  k=k+1,rev)
  }
  return(list(small=list(coord=small_coord), big=list(coord=big_coord)))
}

plot_Tri_kite<-function(a,b,r,theta,rev=FALSE){
  x<-c(0,-1*cos(pi*3/5),-2*cos(pi*3/5))
  y<-c(0,sin(pi*3/5),0)
  x1<-(r*x*cos(theta)-r*y*sin(theta))+a
  y1<-(r*x*sin(theta)+r*y*cos(theta))+b  
  coord<-list(x=x1,y=y1)
  decompose_kite(list(coord=coord),1,rev)
  return(list(coord=coord))
}


decompose_dart<-function(t1,k,rev=FALSE){  
  if(rev){
    x1<-c((t1$coord$x[3]+t1$coord$x[1]/PHI)/PHI,t1$coord$x[1],t1$coord$x[2])
    y1<-c((t1$coord$y[3]+t1$coord$y[1]/PHI)/PHI,t1$coord$y[1],t1$coord$y[2])
    x2<-c(t1$coord$x[3],(t1$coord$x[3]+t1$coord$x[1]/PHI)/PHI,t1$coord$x[2])
    y2<-c(t1$coord$y[3],(t1$coord$y[3]+t1$coord$y[1]/PHI)/PHI,t1$coord$y[2])
  }else{
    x1<-c(t1$coord$x[2],t1$coord$x[3],(t1$coord$x[3]/PHI+t1$coord$x[1])/PHI)
    y1<-c(t1$coord$y[2],t1$coord$y[3],(t1$coord$y[3]/PHI+t1$coord$y[1])/PHI)
    x2<-c(t1$coord$x[2],(t1$coord$x[3]/PHI+t1$coord$x[1])/PHI,t1$coord$x[1])
    y2<-c(t1$coord$y[2],(t1$coord$y[3]/PHI+t1$coord$y[1])/PHI,t1$coord$y[1])
  }
  small_coord<-list(x=x1,y=y1)
  big_coord<-list(x=x2,y=y2)
  
  polygon(small_coord$x, small_coord$y, col = "green")
  polygon(big_coord$x, big_coord$y, col = "red")
  polygon(-small_coord$x, small_coord$y, col = "green")
  polygon(-big_coord$x, big_coord$y, col = "red")
  
  if(k<L){
    decompose_dart(list(coord=big_coord), k=k+1)
    decompose_kite(list(coord=small_coord), k=k+1)
  }
  return(list(small=list(coord=small_coord), big=list(coord=big_coord)))
}

plot(0,0, xlim=c(-10,10), ylim=c(-10,10))
t1<-plot_Tri_dart(0,0,5,0)

plot_Tri_dart<-function(a,b,r,theta,rev=FALSE){
  x<-c(0,cos(pi*1/5),2*cos(pi*1/5))
  y<-c(0,sin(pi*1/5),0)
  x1<-(r*x*cos(theta)-r*y*sin(theta))+a
  y1<-(r*x*sin(theta)+r*y*cos(theta))+b 
  
  if(rev){coord<-list(x=x1[3:1],y=y1[3:1])}
  else{coord<-list(x=x1,y=y1)}

  decompose_dart(list(coord=coord),1)
  return(list(coord=coord))
}


plot(0,0, xlim=c(-5,5), ylim=c(-1,9))
t1<-plot_Tri_kite(0,0,1,pi/10)
t1<-plot_Tri_dart(t1$coord$x[2],t1$coord$y[2],1/PHI,-3*pi/10)
t2<-t1
t1<-plot_Tri_kite(t1$coord$x[2],t1$coord$y[2],1,-5*pi/10)

t2<-plot_Tri_kite(t2$coord$x[1],t2$coord$y[1]+1/PHI,1,-5*pi/10)
t3<-t2
t2<-plot_Tri_dart(t2$coord$x[2],t2$coord$y[2],1/PHI,-9*pi/10)
t2<-plot_Tri_kite(t2$coord$x[1],t2$coord$y[1],1,-7*pi/10)

t3<-plot_Tri_dart(t3$coord$x[1],t3$coord$y[1]+1,1/PHI,-5*pi/10)
t4<-t3
t3<-plot_Tri_kite(t3$coord$x[2],t3$coord$y[2],1,-7*pi/10)
t3<-plot_Tri_dart(t3$coord$x[1],t3$coord$y[1],1/PHI,-3*pi/10)
t3<-plot_Tri_dart(t3$coord$x[2]+cos(3*pi/10)/PHI,t3$coord$y[2]+sin(3*pi/10)/PHI,1/PHI,-9*pi/10)

t4<-plot_Tri_dart(t4$coord$x[1],t4$coord$y[1],1/PHI,-1*pi/10)
t5<-t4
t4<-plot_Tri_dart(t4$coord$x[3],t4$coord$y[3],1/PHI,5*pi/10)
t4<-plot_Tri_dart(t4$coord$x[1],t4$coord$y[1],1/PHI,9*pi/10)

t5<-plot_Tri_kite(t5$coord$x[1],t5$coord$y[1],1,pi/10)
t5<-plot_Tri_kite(t5$coord$x[3],t5$coord$y[3],1,3*pi/10)
t5<-plot_Tri_dart(t5$coord$x[2],t5$coord$y[2],1/PHI,-1*pi/10)
t6<-t5
t5<-plot_Tri_dart(t5$coord$x[3],t5$coord$y[3],1/PHI,5*pi/10)
t5<-plot_Tri_kite(t5$coord$x[1],t5$coord$y[1],1,pi/10)
t5<-plot_Tri_kite(t5$coord$x[3],t5$coord$y[3],1,3*pi/10)

t5<-t6
t5<-plot_Tri_kite(t5$coord$x[1],t5$coord$y[1],1,pi/10)
t5<-plot_Tri_kite(t5$coord$x[3],t5$coord$y[3],1,3*pi/10)
t5<-plot_Tri_dart(t5$coord$x[2],t5$coord$y[2],1/PHI,-1*pi/10)
t6<-t5
t5<-plot_Tri_dart(t5$coord$x[3],t5$coord$y[3],1/PHI,5*pi/10)
t5<-plot_Tri_kite(t5$coord$x[1],t5$coord$y[1],1,pi/10)
t5<-plot_Tri_kite(t5$coord$x[3],t5$coord$y[3],1,3*pi/10)

t6<-plot_Tri_kite(t6$coord$x[1],t6$coord$y[1],1,pi/10)
t6<-plot_Tri_kite(t6$coord$x[3],t6$coord$y[3],1,3*pi/10)
t6<-plot_Tri_dart(t6$coord$x[2],t6$coord$y[2],1/PHI,-1*pi/10)
t6<-plot_Tri_dart(t6$coord$x[1],t6$coord$y[1],1/PHI,3*pi/10)
t7<-plot_Tri_dart(t6$coord$x[3],t6$coord$y[3],1/PHI,-7*pi/10)
t6<-plot_Tri_kite(t6$coord$x[2],t6$coord$y[2],1,pi/10)


実行結果が以下の通り、



ところどころロビンソンの三角形がむき出しになっていますね。。。
失敗したようなので来年の宿題としてリベンジしたいところです。