ぱらぱらめくる『射影幾何学入門〜生物の形態と数学〜

射影幾何学入門―生物の形態と数学

• ユークリッド幾何は、長さが大事。平行線が交わらないことも大事
• 射影幾何は長さはどうでもよくなって、平行線は無限遠で交わり、すべての直線が交点を持つ
  • 有限幾何が、点と線(と面と…)と組合せ論的関係性の抽象化したものであるように
  • 射影幾何も「ユークリッド幾何」で慣れ親しんでいるもろもろの関係性を部分的に抽出したものになっている
• 射影幾何では円も楕円も双曲線も放物線も相互に移行し合える「同じもの〜二次曲線」となる。三角形は直角三角形も鋭角三角形も鈍角三角形も「三角形」として同じものになる
• 射影は透視絵と関係していて、あるd次元空間の点をある視点・光源から見て、d-1次元(超)平面に存在するものとする
• 射影変換は、d次元空間において、d-1次元(超)平面の点(の集合)を別のd-1次元(超)平面に射影写像して、今度は視点・光源をずらして、また、次のd-1次元(超)平面に移す。これを２回以上繰り返して、もとのd-1次元(超)平面に戻してやることができる。そうすると、もとのd-1次元(超)平面上の「点の集合(形に意味があったりする)」が、「別の点の集合」に移される。移されるけれど、「射影幾何的には同じもの〜円が楕円に変わっていたりするけれど、それが『違う』のはユークリッドな感覚の場合であって、射影幾何的には、『同じ』で、『射影幾何的「合同」変換」しただけです」と言う
• そんな射影変換で説明できる生物形態があるよ、という話
• まずは、射影写像をRでやろう
  • 準備。任意次元空間で、法線ベクトルを指定して、(1,0,0,...)というベクトルをその法線ベクトルに回転して移す正規直交基底行列・回転行列を作れるようにしておく(こちらで作った『遺産』)

# あるベクトルを１つの軸とする正規直交基底
# vは接ベクトル
base.v <- function(v){
# 次元d
	d <- length(v)
# 単位ベクトル化
	v <- v/sqrt(sum(v^2))
	ret <- matrix(0,d,d)
# 第１列は接ベクトルそのもの
# 第２列以降は次の通り
# 第i列については、1,2,...,(i-1)行の成分はvの1,2,...,(i-1)成分に
# 比例した値としておき
# 第i行成分を接ベクトルと直交するための調節成分とする
# このようにすると、求めたいn-1個のベクトルと接ベクトルとの内積の計算と
# 求めたいn-1個のベクトル同士の内積の計算とが、同じになるので、
# 求めたいベクトルの値を接ベクトルとの関係で求めると、
# 漏れなく、他のベクトルとの関係も調整されることを使う
	for(i in 2:(d)){
		ret[1:(i-1),i] <- v[1:(i-1)]
		ret[i,i] <- -sum(v[1:(i-1)]^2)/v[i]
	}
	ret[,1] <- v
# 最後にノルムをそろえて回転行列化する
	q <- sqrt(apply(ret^2,2,sum))
	t(t(ret)/q)
}
out <- base.v(runif(4))
out %*% t(out) # 単位ベクトルが得られる

• • 射影しよう

# 射影空間の次元 d
# d+1次元空間にd1次元超平面が1つ(p)ある
# xはd+1次元空間の点を束ねた(d+1)行、n.pt列の行列
# Lはd+1次元空間の点を
# Lはpの上にないd+1次元空間の光源(視点)を表す点座標
# pはp0を通り、vを法線とする超平面
my.projection <- function(x,L,p0,v){
	if(!is.matrix(x)){
		x <- matrix(x,ncol=1)
	}
# 第１列がvであるような回転行列
	base1 <- base.v(v)
# x,Lを平行移動し、さらに逆回転する
# その空間では「マップ先」は第１成分が0であるような超平面になっているので…
	x. <- x-p0
	L. <- L-p0
	x.. <- t(base1) %*% x.
	L.. <- t(base1) %*% L.
	k <- -L..[1,]/(x..[1,]-L..[1,])
	X <- matrix(0,length(L),length(x[1,]))
	X <- c(L..) + t(k * t(x..-c(L..)))
	X. <- base1 %*% X + p0
	X.
}

• • 使ってみる。まず２次元

v <- c(1,2)
p0 <- c(1,1)
L <- c(4,5)
t <- seq(from=0,to=1,length=100)
xs <- t*1
a <- 3
b <- 2
ys <- a*xs+b

XYs <- rbind(xs,ys)

p.XYs <- my.projection(XYs,L,p0,v)
xlim <- ylim <- range(c(c(p.XYs),c(XYs)))
plot(t(p.XYs),type="b",asp=TRUE,xlim=xlim,ylim=ylim)
points(t(XYs),col=2)

• • • ３次元

v <- runif(3)
p0 <- runif(3)
L <- runif(3)
library(MCMCpack)
r <- rdirichlet(1000,rep(0.4,3))

xvs <- matrix(rnorm(3*3),3,3)

Xs <- xvs %*% t(r)
p.Xs <- my.projection(Xs,L,p0,v)

library(rgl)
all.Xs <- cbind(Xs,p.Xs)
col <- rep(1:2,each=length(Xs[1,]))
plot3d(t(all.Xs),col=col)

plot3d(t(Xs))
plot3d(t(p.Xs))