フィボナッチ螺旋とボロノイ図

• で、に黄金角に近い値を用いて、n=1,2,...で打点する。
• その点分布に対してボロノイ図を描く

n <- 200 # 点の数
x <- (1+sqrt(5))/2 # 黄金比
delta <- 0.1 # 黄金比からちょっとずらす
phi <- 2*pi * (x + delta) # 黄金角
x <- sqrt(1:n) * cos((1:n) * phi) # 和田さんのx座標
y <- sqrt(1:n) * sin((1:n) * phi) # 和田さんのy座標

X <- cbind(x,y) # ２次元座標

plot(X,pch=20,cex=0.1) # プロットする。pchは点の形を指定する、cexは点の大きさを指定する
install.packages(c("deldir","sp")) # 二つのパッケージをインストールする
# http://stackoverflow.com/questions/9403660/how-to-create-thiessen-polygons-from-points-using-r-packages 
# ボロノイ図を描く関数をウェブから借りてくる
voronoipolygons <- function(x) {
  require(deldir)
  require(sp)
  if (.hasSlot(x, 'coords')) {
    crds <- x@coords  
  } else crds <- x
  z <- deldir(crds[,1], crds[,2])
  w <- tile.list(z)
  polys <- vector(mode='list', length=length(w))
  for (i in seq(along=polys)) {
    pcrds <- cbind(w[[i]]$x, w[[i]]$y)
    pcrds <- rbind(pcrds, pcrds[1,])
    polys[[i]] <- Polygons(list(Polygon(pcrds)), ID=as.character(i))
  }
  SP <- SpatialPolygons(polys)
  voronoi <- SpatialPolygonsDataFrame(SP, data=data.frame(x=crds[,1],
    y=crds[,2], row.names=sapply(slot(SP, 'polygons'), 
    function(x) slot(x, 'ID'))))
}
# ボロノイ多角形情報を取る
v <- voronoipolygons(X)
# ボロノイ多角形を描く
plot(v)
# 点を打つ
points(X,pch=20,cex=0.1,col="red") # 点の色はred