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

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

R 射影幾何 スプライン曲線 円
  • 一昨日の記事で円を使ったスプライン曲線を引いてみた
  • まあまあいい感じだったのだが、その曲線を考えるにあたり、「曲率中心」がどのように動くか、ということが気になった
  • 曲率の中心が反転する部分というのは、一度直線化するわけだが、それは曲率中心が無限遠になること。無限遠をうまく扱うには射影幾何がよい、ということで、射影幾何版にするとどんなよいことがあるか(ないか)調べてみることにする

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

# 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)
}
# その逆
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)
}
t <- seq(from=0,to=1,length=10000)*50
x <- 0.5*sin(t)
y <- cos(1*t)
X <- cbind(x,y)
X. <- t(apply(X,1,n2sphere))
X.. <- t(apply(X.,1,sphere2n))

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)
points3d(X.,col=2)


plot(X)

# 垂直なベクトルを取り出す

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))
}
# 球面の３点を通る円
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))
}


X <- matrix(rnorm(3*3),ncol=3)
X <- X/sqrt(apply(X^2,1,sum))
my.circle.sphere(X)

mcs <- my.circle.sphere(X)
t <- seq(from=0,to=1,length=100)*2*pi
XYZ <- cbind(mcs$r*cos(t),mcs$r*sin(t),rep(mcs$z,length(t)))
XYZ. <- t(mcs$M %*% t(XYZ))

n.r <- 1000
R <- matrix(rnorm(n.r*3),ncol=3)
R <- R/sqrt(apply(R^2,1,sum))
XYZ. <- t((mcs$M) %*% t(XYZ))
plot3d(R)

points3d(XYZ.,col=3)
points3d(X,col=2,size=10)


t <- seq(from=0,to=1,length=100)
x <- exp(t)*cos(10*t)
y <- exp(t)*sin(10*t)
X <- cbind(x,y)
plot(X,type="b")

X. <- t(apply(X,1,n2sphere))
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
R <- rbind(R,c(0,0,-0.1))
plot3d(R)
lines3d(X.,col=2)
planes3d(0,0,1,0,alpha=0.3)

circles <- list()
X.. <- X.
X..[,3] <- X..[,3] - 0.5
X.. <- X.. * 2
for(i in 1:(length(t)-2)){
	circles[[i]] <- my.circle.sphere(X..[i:(i+2),])
}
ctr <- matrix(0,length(circles),3)
for(i in 1:length(ctr[,1])){
	ctr[i,] <- circles[[i]]$M[,3] * circles[[i]]$z
}
ctr. <- ctr/2
ctr.[,3] <- ctr.[,3] + 0.5
lines3d(ctr,col=3)


##
t <- seq(from=0,to=0.3,length=1000)
t <- t[-1]
x <- t
y <- exp(t)*sin(10*t)
X <- cbind(x,y)
plot(X,type="b")

X. <- t(apply(X,1,n2sphere))
X.p <- X.
X.p[,3] <- X.p[,3]-0.5

apply(X.p^2,1,sum)
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
R <- rbind(R,c(0,0,-0.1))
plot3d(R)
lines3d(X.,col=2)
planes3d(0,0,1,0,alpha=0.3)

circles <- list()
X.. <- X.
X..[,3] <- X..[,3] - 0.5
X.. <- X.. * 2
for(i in 1:(length(t)-2)){
	circles[[i]] <- my.circle.sphere(X..[i:(i+2),])
}
ctr <- matrix(0,length(circles),3)
for(i in 1:length(ctr[,1])){
	ctr[i,] <- circles[[i]]$M[,3]
}
ctr. <- ctr/2
ctr.[,3] <- ctr.[,3] + 0.5
points3d(ctr.,col=3,size=5)
lines3d(ctr.,col=4)