遺伝学・遺伝統計学関連の姉妹ブログ『ryamadaの遺伝学・遺伝統計学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

その２〜射影幾何的に円を使ってスプライン曲線

R 射影幾何 スプライン曲線 円
  • 昨日の記事は、まあ、ユーティリティ関数をいじるだけ、今日はユーティリティ関数にコメントを加えたり、少し、いろいろと書いておこう
  • ２次元平面の曲線を、球面に変換する
# 二次元平面上の3点を通る円
# Xは3x2行列
# 返り値は中心座標：ctr、半径：R、３点の角度：angles、３点の角度が順周りか逆回りか：angle.dir

my.circle <- function(X){
	x1 <- X[1,1]
	x2 <- X[2,1]
	x3 <- X[3,1]
	y1 <- X[1,2]
	y2 <- X[2,2]
	y3 <- X[3,2]
	rot <- 0
	if(y2-y1 ==0){
		rot <- 1
	x1 <- X[1,2]
	x2 <- X[2,2]
	x3 <- X[3,2]
	y1 <- X[1,1]
	y2 <- X[2,1]
	y3 <- X[3,1]
	}
	x <- 1/2 * 1/(y1*(x2-x3)+y2*(x3-x1)+y3*(x1-x2)) * ((y1-y2)*(y2-y3)*(y3-y1)-((y2-y3)*x1^2+(y3-y1)*x2^2+(y1-y2)*x3^2))
	y <- -(x2-x1)/(y2-y1)*(x-(x1+x2)/2) + (y1+y2)/2
	ctr <- c(x,y)
	if(rot==1){
		ctr <- c(y,x)
	}
	# 原点と３点とを結ぶ半径ベクトルの角を求める
	tmp <- (t(t(X)-ctr))
	angles <- Arg(tmp[,1] + 1i*tmp[,2])
	angles[which(angles<0)] <- angles[which(angles<0)] + 2*pi
	#print(angles)
	angle.dir <- 1
	if(angles[2] > angles[1]){
		if(angles[3] > angles[2]){
			
		}else{
			if(angles[3] < angles[1]){
				angles[3] <- angles[3] + 2*pi
			}else{
				angle.dir <- -1
				angles[2] <- angles[2] - 2*pi
				angles[3] <- angles[3] - 2*pi
			}
		}
	}else{
		if(angles[3] < angles[2]){
			angle.dir <- -1
		}else{
			if(angles[3] > angles[1]){
				angle.dir <- -1
				angles[3] <- angles[3] - 2*pi
			}else{
				angles[2] <- angles[2] + 2*pi
				angles[3] <- angles[3] + 2*pi
			}
		}
	}
	#print(angles)
	R <- sqrt(sum((X[1,]-ctr)^2))
	return(list(ctr=ctr,R=R,angles=angles,angle.dir=angle.dir))
}

# 射影幾何的な円によるスプライン曲線

# 2次元カーブを考えるとき、2+1次元空間にあって、
# (0,0,0.5)を中心とする、半径0.5の球を考える

# n次元座標からn+1次元球面座標への変換
n2sphere <- function(x){
	r <- sqrt(sum(x^2))
	z <- r^2/(r^2+1)
	k <- sqrt(z*(1-z))
	ret <- c(k/r*x,z)
	ret
}
# その逆
sphere2n <- function(x){
	ret <- rep(0,length(x)-1)
	if(x[length(x)]==1){
		ret <- rep(Inf,length(ret))
	}else{
		tmp <- x
		tmp[length(x)] <- -(tmp[length(x)]-1)
		ret <- tmp[1:(length(x)-1)]/tmp[length(x)]
	}
	return(ret)
}
# 2次元平面に曲線を描いて
t <- seq(from=0,to=1,length=10000)*50
x <- exp(0.1*t)
y <- exp(0.1*t)*cos(1*t)
X <- cbind(x,y)
# 平面曲線
plot(X)
# 球面上に変換
library(rgl)
n.r <- 1000
R <- matrix(rnorm(n.r*3),ncol=3)
R <- R/sqrt(apply(R^2,1,sum))/2
R[,3] <- R[,3]+0.5
plot3d(R)
X. <- t(apply(X,1,n2sphere))

points3d(X.,col=2)


# 球面座標に変換してからもとに戻してもきちんと戻ることの確認
X. <- t(apply(X,1,n2sphere))
X.. <- t(apply(X.,1,sphere2n))

plot(X..)
  • (0,0,0.5)を中心とする半径0.5の球上の３点
# 垂直なベクトルを取り出す

my.perpendic <- function(X){
	n <- length(X[1,])
	tmp.m <- X[,1:(n-1)]
	if(det(tmp.m)!=0){
		ret <- solve(tmp.m) %*% ((-1)*X[,n])
		ret <- c(ret,1)
	}else{
		tmp.m <- X[,2:n]
		#if(det(tmp.m)==0){
		#	return(c(0,0,0))
		#}
		ret <- solve(tmp.m) %*% ((-1)*X[,1])
		ret <- c(1,ret)
	}
	ret/sqrt(sum(ret^2))
}
# 曲線を展開する球((0,0,0.5)を中心とする半径0.5の球)の表面の３点を通る円の
# 半径rと、円を乗せる平面の中心からの距離zと、その平面をx1,x2平面とする正規直交基底M
my.circle.sphere <- function(X){
	v1 <- X[1,]-X[3,]
	v2 <- X[2,]-X[3,]
	v1 <- v1/sqrt(sum(v1^2))
	v2 <- v2/sqrt(sum(v2^2))
	tmp.X <- rbind(v1,v2)
	v3 <- my.perpendic(tmp.X)
	v2. <- my.perpendic(rbind(v1,v3))
	R <- sqrt(sum(X[1,]^2))
	M <- cbind(v1,v2.,v3)
	
	ip <- sum(v3 * X[1,]/R)
	
	return(list(M=M,r=sqrt(1-ip^2),z=ip))
}
# (0,0,0.5)を中心とする球面上の３点を作る
X <- matrix(rnorm(6),ncol=2)
X2 <- t(apply(X,1,n2sphere))
# 原点中心の半径１の円に変換
X2[,3] <- X2[,3] - 0.5
X2 <- X2 * 2
# その３点を通る円を規定する情報を取り出す
tmp.out <- my.circle.sphere(X2)
# それに基づいて円を３次元空間に描く
t <- seq(from=0,to=1,length=100)*2*pi
x <- tmp.out$r*cos(t)
y <- tmp.out$r*sin(t)
z <- rep(tmp.out$z,length(t))
XYZ <- cbind(x,y,z)
XYZ. <- t(tmp.out$M %*% t(XYZ))
plot3d(XYZ.)
points3d(X2,col=2,size=20)
  • 一般にベクトル空間に３点あれば、その３点を通る円は確定する。３点を通る円の「第１点」から「第２点」までに円弧を引いてみよう
# 一般に３点を通る円
my.circle.n <- function(X,n=100){
	v1.ori <- X[1,]-X[3,]
	v2.ori <- X[2,]-X[3,]
	
	v1 <- v1.ori/sqrt(sum(v1.ori^2))
	v2 <- v2.ori/sqrt(sum(v2.ori^2))
	tmp.X <- rbind(v1,v2)
	v3 <- my.perpendic(tmp.X)
	v2. <- my.perpendic(rbind(v1,v3))
	#R <- sqrt(sum(X[1,]^2))
	M <- rbind(v1,v2.,v3)
	coords <- M %*% cbind(v1.ori,v2.ori)
	Y <- rbind(t(coords[1:2,1:2]),c(0,0))
	out <- my.circle(Y)
	ctr <- c(out$ctr,0)
	t <- seq(from=out$angles[1],to=out$angles[2],length=100)
	#t <- seq(from=0,to=2*pi,length=100)
	x <- out$R*cos(t)+out$ctr[1]
	y <- out$R*sin(t)+out$ctr[2]

	tmp2 <-t(t(M) %*% rbind(x,y,rep(0,length(t))) + X[3,])
	tmp2
}

X <- diag(rep(1,3))
X <- matrix(rnorm(9),3,3)
out <- my.circle.n(X)
plot3d(out)
points3d(X,col=2,size=5)
  • 似たような関数
# 2ベクトルを与えv1、v2を含む平面にv1を角０、v2のある方向に角をとって描く円の座標を返す
my.rotation.2vec <- function(v1,v2,theta=seq(from=0,to=1,length=100)*2*pi){
	L1 <- sqrt(sum(v1^2))
	L2 <- sqrt(sum(v2^2))

	v1. <- v1/L1
	v2. <- v2/L2plot.sphere.curve <- function(X,X.line=TRUE,Ctr=TRUE,Ctr.curve=TRUE,Circle=TRUE,n.r=1000,n.t=50,z0=TRUE,alphar=0.3,alphaz0=0.3){
	n.pt <- length(X[,1])
	X. <- t(apply(X,1,n2sphere))
	R <- matrix(rnorm(n.r*3),ncol=3)
	R <- R/sqrt(apply(R^2,1,sum))/2
	R[,3] <- R[,3]+0.5
	R <- rbind(R,c(0,0,-0.1))

	circles <- list()
	X.. <- X.
	X..[,3] <- X..[,3] - 0.5
	X.. <- X.. * 2
	for(i in 1:(n.pt-2)){
		circles[[i]] <- my.circle.sphere(X..[i:(i+2),])
	}
	ctr <- matrix(0,length(circles),3)
	t <- seq(from=0,to=1,length=n.t)*2*pi
	cic <- list()
	for(i in 1:length(ctr[,1])){
		ctr[i,] <- circles[[i]]$M[,3]*circles[[i]]$z
		tmp <- cbind(circles[[i]]$r*cos(t),circles[[i]]$r*sin(t),rep(circles[[i]]$z,length(t)))
		tmp <- t(circles[[i]]$M %*% t(tmp))
		cic[[i]] <- tmp
		cic[[i]] <- cic[[i]]/2
		cic[[i]][,3] <- cic[[i]][,3] + 0.5
	}
	ctr. <- ctr/2
	ctr.[,3] <- ctr.[,3] + 0.5
	
	plot3d(R,alpha=alphar)
	if(z0){
		planes3d(0,0,1,0,alpha=0.3)
	}
	if(Ctr){
		points3d(ctr.,col=3,size=5)
	}
	if(Ctr.curve){
		for(i in 1:(length(ctr.[,1])-1)){
			tmp.X <- rbind(ctr.[i,],ctr.[i+1,],c(0,0,0.5))
			out <- my.circle.n(tmp.X)
			lines3d(out,col=4)
		}
	}
	if(Circle){
		for(i in 1:length(cic)){
			lines3d(cic[[i]],col=6)
		}
	}
	
	points3d(X.,col=2,size=5)
	if(X.line){
		lines3d(X.,col=2)
	}

	plot(X,pch=20,type="b")
	points(t(apply(ctr.,1,sphere2n)),col=2)
	for(i in 1:length(cic)){
		points(t(apply(cic[[i]],1,sphere2n)),col=3,type="l")
	}
	lines(t(apply(ctr.,1,sphere2n)),col=2)


}

	v3. <- v2.-sum(v1.*v2.)*v1.
	v3. <- v3./sqrt(sum(v3.^2))
	x <- y <- rep(0,length(theta))
	x <- cos(theta)
	y <- sin(theta)
	m <- cbind(v1.,v3.)
	tmp <- t(m %*% rbind(x,y))
	print(L1)
	tmp * L1
}
v1 <- c(1,0,0)
v2 <- c(0,0,1)
X <- my.rotation.2vec(v1,v2)
plot3d(X)
  • さて。
  • 曲線を反映する2次元平面の点列が得られたとき、それを3次元球面に移して、その上での曲線ととらえつつ、滑らかに移行するように、曲率中心を「合い隣り合う２点の球面上に対応する２点を通る２つの球面上の円を取り、その２円のうちの第一のものから第二のものへと、その２点を通るままに移行させることを考える。そのような円の中心はもちろん連続的に変化する。その様子を描く関数

plot.sphere.curve <- function(X,X.line=TRUE,Ctr=TRUE,Ctr.curve=TRUE,Circle=TRUE,n.r=1000,n.t=50,z0=TRUE,alphar=0.3,alphaz0=0.3){
	n.pt <- length(X[,1])
	X. <- t(apply(X,1,n2sphere))
	R <- matrix(rnorm(n.r*3),ncol=3)
	R <- R/sqrt(apply(R^2,1,sum))/2
	R[,3] <- R[,3]+0.5
	R <- rbind(R,c(0,0,-0.1))

	circles <- list()
	X.. <- X.
	X..[,3] <- X..[,3] - 0.5
	X.. <- X.. * 2
	for(i in 1:(n.pt-2)){
		circles[[i]] <- my.circle.sphere(X..[i:(i+2),])
	}
	ctr <- matrix(0,length(circles),3)
	t <- seq(from=0,to=1,length=n.t)*2*pi
	cic <- list()
	for(i in 1:length(ctr[,1])){
		ctr[i,] <- circles[[i]]$M[,3]*circles[[i]]$z
		tmp <- cbind(circles[[i]]$r*cos(t),circles[[i]]$r*sin(t),rep(circles[[i]]$z,length(t)))
		tmp <- t(circles[[i]]$M %*% t(tmp))
		cic[[i]] <- tmp
		cic[[i]] <- cic[[i]]/2
		cic[[i]][,3] <- cic[[i]][,3] + 0.5
	}
	ctr. <- ctr/2
	ctr.[,3] <- ctr.[,3] + 0.5
	
	plot3d(R,alpha=alphar)
	if(z0){
		planes3d(0,0,1,0,alpha=0.3)
	}
	if(Ctr){
		points3d(ctr.,col=3,size=5)
	}
	if(Ctr.curve){
		for(i in 1:(length(ctr.[,1])-1)){
			tmp.X <- rbind(ctr.[i,],ctr.[i+1,],c(0,0,0.5))
			out <- my.circle.n(tmp.X)
			lines3d(out,col=4)
		}
	}
	if(Circle){
		for(i in 1:length(cic)){
			lines3d(cic[[i]],col=6)
		}
	}
	
	points3d(X.,col=2,size=5)
	if(X.line){
		lines3d(X.,col=2)
	}

	plot(X,pch=20,type="b")
	points(t(apply(ctr.,1,sphere2n)),col=2)
	for(i in 1:length(cic)){
		points(t(apply(cic[[i]],1,sphere2n)),col=3,type="l")
	}
	lines(t(apply(ctr.,1,sphere2n)),col=2)


}
t <- seq(from=0,to=3,length=200)
t <- t[-1]
x <- exp(t)*cos(10*t)
y <- exp(t)*sin(10*t)
X <- cbind(x,y)
plot(X,type="b")
plot.sphere.curve(X,Circle=FALSE)