ryamadaのコンピュータ・数学メモ

2014-06-07

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

座標の取り方をいじるとぐっと簡単になるのでそうする
斉次座標という
- １次元の座標 $x$ を、２つの数値を使って $(X_1,X_2)$ で表す
- $x=\frac{X_1}{X_2}$ という関係にする
- $X_2=0$ の点は $x=\infty$ 無限遠点になる
２次元の点も $(x_1,x_2)$ を $(X_1,X_2,X_3)$ とする。次元が上がっても同様
この斉次座標を使うと、射影直線による射影変換が
$\begin{pmatrix}X_1'\\ X_2'\end{pmatrix}=\begin{pmatrix}a & b\\ c & d \end{pmatrix} \begin{pmatrix}X_1\\ X_2\end{pmatrix}$ となって
$x' = \frac{a x + b}{c x +d}$ と書くより簡単だし、本質的

my.projection.transform2 <- function(x,M){
	X <- c(x,1)
	X. <- M %*% X
	x. <- X.[1]/X.[2]
	x.
}

M <- matrix(c(1,1,1,1),2,2)

x <- c(1)
n <- 50
xs <- rep(0,n)
xs[1] <- x
for(i in 2:n){
	xs[i] <- my.projection.transform2(xs[i-1],M)
}
plot(xs)

M <- matrix(c(1,0,1,1),2,2)
x <- c(1)
n <- 50
xs <- rep(0,n)
xs[1] <- x
for(i in 2:n){
	xs[i] <- my.projection.transform2(xs[i-1],M)
}
plot(xs)

M <- matrix(c(6,-2,2,1),2,2)
x <- c(0.001)
n <- 50
xs <- rep(0,n)
xs[1] <- x
for(i in 2:n){
	xs[i] <- my.projection.transform2(xs[i-1],M)
}
plot(xs)

M <- matrix(rnorm(4),2,2)
x <- c(0.001)
n <- 50
xs <- rep(0,n)
xs[1] <- x
for(i in 2:n){
	xs[i] <- my.projection.transform2(xs[i-1],M)
}
plot(xs,type="l")

行列の４成分、a,b,c,dによる $(a-d)^2+4bd$ の正負、ゼロが不動点数を決めるので、それに応じて色々な射影変換が得られる

n.iter.sq <- 4
n.iter <- n.iter.sq^2
par(mfrow=c(n.iter.sq,n.iter.sq))

for(ii in 1:n.iter){
M <- matrix(rnorm(4),2,2)
x <- c(0.001)
n <- 50
xs <- rep(0,n)
xs[1] <- x
for(i in 2:n){
	xs[i] <- my.projection.transform2(xs[i-1],M)
}
plot(xs,type="l",xlab="",ylab="",axes=FALSE)

}
par(mfrow=c(1,1))