2013-09-10

フィボナッチ格子を球面に展開

単位正方形にフィボナッチ格子を作ってそれを球面に張ることもできる
それは球面にらせんを描き、そのらせん上に黄金比関連の角間隔で点を打つことに近似できる

# フィボナッチ数を返す関数
FN <- function(n){
	# 負の数のフィボナッチ数は、正のそれに符号をつけたもの
	if(n<0){
		return(both.fn.3(-n)*(-1)^(n+1))
	}
	# 0,1,2はフィボナッチ数発生の種なので指定する
	if(n==0){
		return(0)
	}else if (n==1 | n==2){
		return(1)
	}else{# 3以上は、計算する
		tmp <- rep(0,n)
		tmp[1] <- tmp[2] <- 1
		for(i in 3:n){
			tmp[i] <- tmp[i-1]+tmp[i-2]
		}
		return(tmp[length(tmp)])
	}
}
# 有限法
n <- 10

# 平面格子を0<=x<=1, 0<=y<=1に作る
# 演算子 %% は剰余(ここでは小数点以下を取り出す)
k <- 0:FN(n)
F <- cbind(k/FN(n),(k*FN(n-1)/FN(n))%%1)
# 平面上にプロットしてみる
plot(F)
# [0,1] x [0,1]の２次元格子点座標を３次元に変換
# sqrt(x-x^2)は 0<=x<=1なので問題なく納まる
# F[,1]は２次元格子点のx座標、F[,2]はそのy座標
S <- cbind(2*cos(2*pi*F[,2])*sqrt(F[,1]-F[,1]^2),2*sin(2*pi*F[,2])*sqrt(F[,1]-F[,1]^2),1-2*F[,1])

# ３次元プロットのためのパッケージrglをインストールして読み込み
install.packages("rgl")
library(rgl)
# ３次元プロット
plot3d(S)

# 球面らせん法
# f1は黄金比
fib.lattice.S2 <-function(N,f1=(sqrt(5)+1)/2,k=1){
	f2 <- f1-1 # 黄金比はx^2-x-1=0の１つの根。もう１つの根を取り出す
	# 点の数2N+1
	P <- 2*N+1
	# -N から Nまでの整数列
	i <- seq(from=-N, to=N,by=k)
	# zの座標は-1 から 1の均等割り
	theta <- asin(2*i/(2*N+1))
	z <- sin(theta)
	# そのようなz座標に対応する、(x,y)座標を取り出すにあたり
	# 黄金比から作ったf2に応じてxy平面での角座標を作る
	phi <- 2*pi*i*f2
	# その角を使って、z座標に対応した半径の座標とする
	x <- cos(theta)*cos(phi)
	y <- cos(theta)*sin(phi)
	return(cbind(x,y,z))
}
# FN(n)個の点をつくるために、引数をFN(n)/2とする
S2 <- fib.lattice.S2(FN(n)/2)
plot3d(S2)

# 両者の値は、対応づいている
plot(c(S),c(S2))

ちょっとメモ

# FN(n)個の点をつくるために、引数をFN(n)/2とする
n <- 15
#S2 <- fib.lattice.S2(FN(n)/2)
nn <- 5000
S2 <- fib.lattice.S2(nn)
plot3d(S2)

# 両者の値は、対応づいている
#plot(c(S),c(S2))

dist.S2 <- as.matrix(dist(S2))

ord.dist.S2 <- apply(dist.S2,1,order)

E <- t(ord.dist.S2[2:5,])

sorted.dist.S2 <- apply(dist.S2,1,sort)
matplot(t(sorted.dist.S2[1:6,]),type="l")
matplot(t(sorted.dist.S2[1:8,]),type="l")

adj.m <- dist.S2
thres <- max(sorted.dist.S2[4,])
thres <- max(sorted.dist.S2[5,])
thres <- mean(sorted.dist.S2[5,length(S2[,1])/2],sorted.dist.S2[6,length(S2[,1])/2])

thres <- mean(sorted.dist.S2[6,(length(S2[,1])/3):(length(S2[,1])*2/3)],sorted.dist.S2[6,length(S2[,1])/2])
thres <- max(sorted.dist.S2[5,(length(S2[,1])/3):(length(S2[,1])*2/3)])
thres <- min(sorted.dist.S2[6,])
thres <- mean(sorted.dist.S2[5:6,])
#thres <- sorted.dist.S2[5,length(S2[,1])/2]
#adj.m <- adj.m <= thres
adj.m <- adj.m < thres
diag(adj.m) <- 0
el <- which(adj.m==1,arr.ind=TRUE)

plot3d(S2)
for(i in 1:length(el[,1])){
	#segments3d(rbind(S2[el[i,1],],S2[el[i,2],]))
}
tabulate(apply(adj.m,1,sum))

g <- graph.adjacency(adj.m)
plot(g)

apply(adj.m,2,sum)

#adj.m <- matrix(0,length(S2[,1]),length(S2[,1]))
#el <- matrix(0,length(E),2)
#cnt <- 1
for(i in 1:length(E[,1])){
	for(j in 1:length(E[i,])){
		adj.m[i,E[i,j]] <- 1
		el[cnt,] <- c(i,E[i,j])
		cnt <- cnt+1
	}
}
library(igraph)
g <- graph.adjacency(adj.m)
plot(g)
range(apply(adj.m,1,sum))
range(adj.m-t(adj.m))

plot3d(S2)
for(i in 1:length(el[,1])){
	segments3d(rbind(S2[el[i,1],],S2[el[i,2],]))
}