多次元図形 - ryamadaの遺伝学・遺伝統計学メモ

多次元多様体を経時撮影するデータをシミュレーションしたい
まず、球とそれの連続変形体を作るとしよう
一つの作り方は、球をその法線方向に伸び縮みさせる関数を周期関数を用いて滑らかに作ることにする。3次元の例

# 極座標化する
my.sp.coords <- function(x){
	d <- length(x)
	r <- sqrt(sum(x^2))
	thetas <- rep(0,d-1)
	tmp <- x
	for(i in 1:(d-1)){
		thetas[i] <- asin(tmp[length(tmp)]/sqrt(sum(tmp^2)))
		tmp <- tmp[-length(tmp)]
	}
	return(c(thetas[(d-1):1],r))
}

my.shape <- function(X,as,bs,r){
	X.sp <- t(apply(X,1,my.sp.coords))
	d <- length(X[1,])
	R <- rep(1,length(X[,1]))
	for(i in 1:(d-1)){
		for(j in 1:length(as[[i]])){
			R <- R + as[[i]][j] * cos(X.sp[,i]*j+bs[[i]][j])
		}
	}
	(X.sp[,d] < R*r)
}

as <- list(c(0.5,0,0,0.2),c(1,0,0,0,0.1))
bs <- list(c(pi/2,0,0,pi/3),c(0,0,0,0,pi/3))

n.pt <- 1000000
d <- 3
X<- matrix(runif(n.pt*d,min=-3,max=3),ncol=d)
InOut <- my.shape(X,as,bs,1)
X. <- X[which(InOut),]
library(rgl)
plot3d(X.)

撮像しよう。
- 撮像は、ある範囲を測定範囲として、それを等間隔グリッドにしたものとする。物体からはシグナルがランダムに出ているが、それがグリッド小立方に当たれば、「あり」そうでなければ「なし」の２値とする（量的観察にしてもよい）

my.nd.image <- function(X,R,n){
	X. <- X
	for(i in 1:length(R)){
		X.[,i] <- (X[,i]-R[[i]][1])/(R[[i]][2]-R[[i]][1])
	}
	X.. <- (X. >= 0) & (X. <= 1)
	s <- apply(X..,1,prod)
	X... <- X.[which(s==1),]*n
	X... <- round(X...,0)
	d <- length(X[1,])
	tmp <- apply(t(X...-1) * (n)^(0:(d-1)),2,sum) + 1
	ret <- tabulate(tmp,nbins=n^d)
	array(ret,rep(n,d))
}
R <- list(c(-3,3),c(-3,3),c(-3,3))
n <- 2^6
x <- seq(from=-3,to=3,length=n+1)

out <- my.nd.image(X.,R,n)
par(ask=TRUE)
for(i in 1:n){
	image(x,x,sign(out[,i,]))
}

さて、動かそう
- 回転と並進だけなら、次元dに対して、(d+1)^2行列で動かせる
- 動かしたうえで、「撮像」すれば、「映画」が出来上がる

# 小さい角の回転行列を作る(rに小さい値をあてる）
my.Random.Start <- function(d,r=1,n=d){
	ret <- diag(rep(1,d))
	for(i in 1:n){
		js <- sample(1:d,2)
		theta <- runif(1)*2*pi*r
		tmp <- diag(rep(0,d))
		tmp[js,js] <- matrix(c(cos(theta),sin(theta),-sin(theta),cos(theta)),2,2)
		for(k in (1:d)[-js]){
			tmp[k,k] <- 1
		}
		ret <- tmp %*% ret
	}
	ret
}
# 回転成分はdxd、行列は(d+1)x(d+1)
M <- my.Random.Start(d,0.02)
M. <- diag(rep(1,d+1))
M.[1:d,1:d] <- M
# 並進成分
shift <- rnorm(d)*0.01
M.[1:d,d+1] <- shift
M.
# 時刻数を、画像のグリッドと同じ数にすれば、超立方体データができる
t <- seq(from=0,to=1,length=n)
Image.hx <- array(0,rep(n,d+1))
obj.hx <- list()
X.. <- cbind(X.,rep(1,length(X.[,1])))
for(i in 1:length(t)){
	if(i==1){
		tmp <- Re(t(M. %*% t(X..)))
	}else{
		tmp <- Re(t(M. %*% t(tmp)))
	}
	
	obj.hx[[i]] <- tmp[,1:d]
}

minmax.obj <- range(unlist(obj.hx))
R <- list(minmax.obj,minmax.obj,minmax.obj)
for(i in 1:length(t)){
	Image.hx[,,,i] <- my.nd.image(obj.hx[[i]],R,n)
}
for(i in 1:length(t)){
	tmp <- rbind(obj.hx[[i]],rep(minmax.obj[1],3),rep(minmax.obj[2],3))
	plot3d(tmp)
}