2 Geometric Applications ぱらぱらめくる『Quaternion Fourier Transforms for Signal and Image Processing』

3Dでは純虚部分を空間３次元に、4Dでは実＋虚４成分を4次元に対応させる。
3D 反転
- 原点から法線ベクトルを定め、その面での反転は、単位法線ベクトルnと3次元の点pを実成分0の四元数で表したうえで、 $npn$ という簡単な計算

n.pt <-5000
# 単位球面を作ってから、変形・移動
X <- matrix(rnorm(n.pt*3),ncol=3)
X <- X/sqrt(apply(X^2,1,sum))
X[,1] <- X[,1]+X[,2]+ X[,3]^4
X <- X + 1.4
library(rgl)

# 正規直交基底を作って、第3ベクトルを法線ベクトルとする
library(GPArotation)
V <- Random.Start(3)
v <- V[,3]
v <- v/sqrt(sum(v^2))

# 形の四元数作成
Xh <- quaternion(n.pt)
for(i in 1:n.pt){
	Im(Xh[i]) <- X[i,]
}
# 法線ベクトルの四元数作成
vh <- quaternion(1)
Im(vh) <- v

Xhv <- vh * Xh * vh

# 反転後の四元数座標を通常の3次元座標に戻す
Xv <- t(as.matrix(Im(Xhv)))
plot3d(rbind(X,rep(min(X),3),rep(max(X),3)))
points3d(Xv[,2:4],col=3)
lines3d(matrix(c(rep(0,3),v),byrow=TRUE,ncol=3),col=2)
# 反転平面をプロット
F <- cbind(rnorm(n.pt),rnorm(n.pt),rep(0,n.pt)) %*% t(V)
points3d(F,col=4)

- 回転は、単位純虚四元数 $q=e^{\mu \phi}$ と表されたとき、 $\mu$ を回転軸とした、角度 $2\phi$ の回転が $q p \bar{q}$ と計算される

# 回転軸と角度を指定して、ベクトルを作る
mu <- vh
phi <- pi/3
vr <- exp(mu * phi)

Xhv2 <- vr * Xh * Conj(vr)
Xv2 <- t(as.matrix(Im(Xhv2)))
plot3d(rbind(X,rep(min(c(X,Xv2)),3),rep(max(c(X,Xv2)),3)))
points3d(Xv2[,2:4],col=3)
points3d(F,col=4)
lines3d(matrix(c(rep(0,3),v),byrow=TRUE,ncol=3),col=2,lw=5)

- 回転は、実は２つのことなるreflectionの連結。だから
  - 二つのreflection面に垂直なベクトルを回転軸とする

n.pt <-5000
# 単位球面を作ってから、変形・移動
X <- matrix(rnorm(n.pt*3),ncol=3)
X <- X/sqrt(apply(X^2,1,sum))
X[,1] <- X[,1]+X[,2]+ X[,3]^4
X <- X + 1.4
library(rgl)

# 正規直交基底を作って、第3ベクトルを法線ベクトルとする
library(GPArotation)
V <- Random.Start(3)
v <- V[,3]
v <- v/sqrt(sum(v^2))
V2 <- Random.Start(3)
v2 <- V2[,3]

# 形の四元数作成
Xh <- quaternion(n.pt)
for(i in 1:n.pt){
	Im(Xh[i]) <- X[i,]
}
# 法線ベクトルの四元数作成
vh <- quaternion(1)
Im(vh) <- v
vh2 <- quaternion(2)
Im(vh2) <- v2
Xhv <- vh2 * (vh * Xh * vh) * vh2

# 反転後の四元数座標を通常の3次元座標に戻す
Xv <- t(as.matrix(Im(Xhv)))
plot3d(rbind(X,rep(min(X),3),rep(max(X),3)))
points3d(Xv[,2:4],col=3)
lines3d(matrix(c(rep(0,3),v),byrow=TRUE,ncol=3),col=2)
# 反転平面をプロット
F <- cbind(rnorm(n.pt),rnorm(n.pt),rep(0,n.pt)) %*% t(V)
F2 <- cbind(rnorm(n.pt),rnorm(n.pt),rep(0,n.pt)) %*% t(V2)
points3d(F,col=4)
points3d(F2,col=3)

- 回転を連続するときには、 $q_2(q_1 p \bar{q_1})\bar{q_2}=(q_2q_1) p \bar{q_2q_1}$ と言うように、回転の合成を四元数の積として扱える
- qの指数関数によって、回転角を連続的にも扱える
- Shearing 斜めに傾ける。二次元で言えば的な
  - 二種類あって・・・
  - Axial-shear
    - 正規直交基底３ベクトルを使って変換する。斜めに伸ばしつつ、回転的成分が入る。斜め関係は、vh1,vh2で、回転系をvh3で。

# 3D shear
vh1 <- vh2 <- vh3 <- quaternion(1)
Im(vh1) <- V[,1]
Im(vh2) <- V[,2]
Im(vh3) <- V[,3]
alpha <- 2
Xhv3 <- Xh * (1-alpha/2*vh3) - alpha/2*vh2*Xh*vh1
Xv3 <- t(as.matrix(Im(Xhv3)))
plot3d(rbind(X,rep(min(c(X,Xv3)),3),rep(max(c(X,Xv3)),3)))
points3d(Xv3[,2:4],col=3)
#points3d(F,col=4)
lines3d(matrix(c(rep(0,3),V[,1]),byrow=TRUE,ncol=3),col=2,lw=5)
lines3d(matrix(c(rep(0,3),V[,2]),byrow=TRUE,ncol=3),col=6,lw=5)
lines3d(matrix(c(rep(0,3),V[,3]),byrow=TRUE,ncol=3),col=5,lw=5)

- - Beam-shear
    - 拡縮係数が二つに増える
    - 計算式はRのソースを参照

alpha <- 2
beta <- 0.6
Xhv4 <- Xh * (1/2+alpha/2*vh3) - alpha/2*vh2*Xh*vh1 + (1/2+beta/2*vh2)*Xh - beta/2*vh1*Xh*vh3
Xv4 <- t(as.matrix(Im(Xhv4)))
plot3d(rbind(X,rep(min(c(X,Xv4)),3),rep(max(c(X,Xv4)),3)))
points3d(Xv4[,2:4],col=3)
#points3d(F,col=4)
lines3d(matrix(c(rep(0,3),V[,1]),byrow=TRUE,ncol=3),col=2,lw=5)
lines3d(matrix(c(rep(0,3),V[,2]),byrow=TRUE,ncol=3),col=6,lw=5)
lines3d(matrix(c(rep(0,3),V[,3]),byrow=TRUE,ncol=3),col=5,lw=5)

- 拡張
  - 拡張も2種類
  - 指定ベクトル方向に $\alpha$ 倍
  - $(1+\alpha)/2 p + (1-\alpha)/2 \mu p \mu$

alpha <- 1.2
Xhv5 <-(1+alpha)/2 * Xh + (1-alpha)/2 * vh1*Xh*vh1

Xv5 <- t(as.matrix(Im(Xhv5)))
plot3d(rbind(X,rep(min(c(X,Xv5)),3),rep(max(c(X,Xv5)),3)))
points3d(Xv5[,2:4],col=3)
lines3d(matrix(c(rep(0,3),V[,1]),byrow=TRUE,ncol=3),col=2,lw=5)

- - 別法
  - 指定ベクトルを法線とする面に関して、原点からの距離を( $\alpha$ )倍する
  - $(1+\alpha)/2 p - (1-\alpha)/2 \mu p \mu$

alpha <- 1
Xhv5 <-(2+alpha)/2 * Xh + (alpha)/2 * vh1*Xh*vh1

Xv5 <- t(as.matrix(Im(Xhv5)))
plot3d(rbind(X,rep(min(c(X,Xv5)),3),rep(max(c(X,Xv5)),3)))
points3d(Xv5[,2:4],col=3)
lines3d(matrix(c(rep(0,3),V[,1]),byrow=TRUE,ncol=3),col=2,lw=5)
lines3d(matrix(c(rep(0,3),V[,2]),byrow=TRUE,ncol=3),col=6,lw=5)
lines3d(matrix(c(rep(0,3),V[,3]),byrow=TRUE,ncol=3),col=5,lw=5)

4Dに移る
- 4D reflection
  - 対象を普通の四元数pとして、reflectionを表す法線ベクトルを単位四元数qとすれば
  - $-q \bar{p}q$ がreflection

# 適当に四元数を作って
p <- rquat(1)
q <- rquat(1)
q <- q/sqrt(Norm(q))
q
p
# reflection
p. <- -q*Conj(p)*q
p.-p
# reflectionする前から後へのベクトルは、reflection法線ベクトルの定数倍
Re(p.-p)/Re(q)
i(p.-p)/i(q)
j(p.-p)/j(q)
k(p.-p)/k(q)

- 4D rotation
  - 二つの単位四元数q,rで回転を定める。
  - $qpr$ が回転

# 適当に二つの4次元ベクトルを作り
p1 <- rquat(1)
p2 <- rquat(1)
# 単位4次元ベクトルを作り
q <- rquat(1)
q <- q/sqrt(Norm(q))
r <- rquat(1)
r <- r/sqrt(Norm(r))

# 回転する
p1.rot <- q * p1 * r
p2.rot <- q * p2 * r
# 回転後に二つのベクトルの長さが保たれ、
Norm(p1)
Norm(p1.rot)
Norm(p2)
Norm(p2.rot)
# 二つのベクトルの先を結んだベクトルの長さも保たれている
# 結局３角形の三辺の長さが保たれているから、三角形が合同な三角形に移っていることになり、角も保存
Norm(p1-p2)
Norm(p1.rot-p2.rot)

- 4D回転も、3D回転と同じで、二つの純虚四元数h1,h2によるreflectionの合成として作れるようなものを考えることがある
  - $h2(\bar{h1 \bar{p} h1})h2 = h2\bar{h1}p(\bar{h1}h2)$ のように。
  - こうすると、実部が等しく、虚部の符合が逆転した二つの単位四元数 $h2\bar{h1},\bar{h1}h2$ による回転となる。このとき、回転角 $2\phi$ は実部S(q)=S(r)と $\cos(\phi)=S(q)=S(r)$ の関係であるという（四元数の形に制約があることに注意)。
  - この回転は $S(h2\bar{h1}p)=S(\bar{h1}h2p=0$ という面に関する回転だという
球面座標
- ３次元単位球面座標を純虚四元数a,bで表すとき、abを含む平面は球の大円を定める。この平面に垂直な単位ベクトル(aから反時計回りにbに向かう)mと、回転角、 $\psi$ とには $\lambda=b a^{-1} = \cos{\psi} + m \sin{\psi}$ として計算される $\lambda$ が取れる
- この $\lambda$ は大円上の弧に対応する四元数であるが、これは、「ある大円」に対して、その相対的位置に関係なく、弧の長さで決まる数である。
- また、球面上の３点がa,b,c(a->b->c)で表されているとき、弧a->cの長さはa->b + b->cと「迂回」すれば長くなるが、知りたいのはa->cを通る大円上の距離であるとすると、これをを使ってとして計算できる
  - その結果、球面上の三角形a->b->c->aとぐるりと戻ると、その距離は $\lambda=1$ となるが、実部が1なので $\cos{\psi}$ の $\psi$ は0)となる。

# 球面座標

# オイラー角で単位球面上の２点を決める
thetas <- runif(2)*2*pi
phis <- runif(2)*2*pi

x1 <- cos(thetas[1])*cos(phis[1])
y1 <- sin(thetas[1])*cos(phis[1])
z1 <- sin(phis[1])
X1 <- c(x1,y1,z1)
x2 <- cos(thetas[2])*cos(phis[2])
y2 <- sin(thetas[2])*cos(phis[2])
z2 <- sin(phis[2])
X2 <- c(x2,y2,z2)
n.pt <- 1000
X <- matrix(rnorm(n.pt*3),ncol=3)
X <- X/sqrt(apply(X^2,1,sum))
library(rgl)
plot3d(X)
points3d(rbind(X1,X2), col=2,size=15)
lines3d(matrix(c(0,0,0,X1),byrow=TRUE,ncol=3),col=2)
lines3d(matrix(c(0,0,0,X2),byrow=TRUE,ncol=3),col=2)

v1 <- v2 <- quaternion(1)
Im(v1) <- X1
Im(v2) <- X2

# lambdaを計算する
v3 <- v2 * (v1^(-1))

# 回転角
psi <- acos(Re(v3))
psi
# ベクトルから計算した内積(cos(psi))
cos12 <- sum(X1*X2)
cos12
# lambdaの実部となっている
Re(v3)
# 虚部が表す法線ベクトルの単位ベクトルを作る
x3 <- i(v3)/sin(psi)
y3 <- j(v3)/sin(psi)
z3 <- k(v3)/sin(psi)
lines3d(matrix(c(0,0,0,x3,y3,z3),byrow=TRUE,ncol=3),col=4)

射影幾何・同次座標
- ３次元用の同次座標を四元数の４つの係数に用いる
- 同次座標は $X,Y,Z,W$ を $X/W,Y/W,Z/W,1$ と見るが、このWを四元数の実部、X,Y,Zを虚部とする
- 無限遠点は実部が0の場合。それは、空間上の点、というよりは、作用方向としての向きを表しているとみてもよい
- ２点を結ぶ直線上の点をパラメタ扱いするときに、通常なら、パラメタの増加の具合をいじらないといけないけれど、２点の実部を変えることによって、２点間の動きに傾斜を発生させることができる
- ４元数の線形な変換は、同次座標の線形変換をもたらし、それは３次元的には、複雑な変換に対応してくれる(行列で３次元のときに４ｘ４行列を使うことで、並行移動(第４列)やその他の共形変換が表せることに対応する)
連立四元数線型方程式
- $f(q) = \sum_{p=1}^P m_p q n_p$ ですべて(m,n,p)が四元数、というものだが
- $f(q) = Aq + Bq Hi + Cq Hj + Dq Hk$ と簡単に書ける特徴がある
- 今、四元数 $p=p0 + p1 i + p2 j + p3 k$ に対して、その「変換行列表示」として、 $pM=\begin{pmatrix}p0 & -p1 & -p2 & -p3\\ p1 & p0 & -p3 & p2 \\ p2 & p3 & p0 %& -p1 \\ p3 & -p2 & p1 & p0 \end{pmatrix}$ とすると
- $f(q) = (AM + BM HiM' + CM HjM' + DM HkM')q$ となると言う。ただし、HiM'は純四元数Hiに対する変換行列表示の2-4 x 2-4 部分の転置行列である
- 実際に確かめてみる

# Linear quaternion function
# f(q) = Aq + Bq Hi + Cq Hj + Dq Hk
# A,B,C,Dは四元数
# ４元数のベクトル表示から行列表示を作る
my.H.mat <- function(h){
	h <- c(Re(h),i(h),j(h),k(h))
	ret <- diag(rep(h[1],4))
	ret[1,2] <- ret[3,4] <- -h[2]
	ret[2,1] <- ret[4,3] <- h[2]
	ret[1,3] <- ret[4,2] <- -h[3]
	ret[2,4] <- ret[3,1] <- h[3]
	ret[1,4] <- ret[2,3] <- -h[4]
	ret[3,2] <- ret[4,1] <- h[4]
	ret
}

# ４元数行列の2:4 x 2:4 部分の転置
my.trans.H3x3 <- function(H){
	ret <- H
	ret[2:4,2:4] <- t(H[2:4,2:4])
	ret
}
# ４つの四元数から行列を作る
my.H.linear.mat <- function(A,B,C,D){
	my.H.mat(A) + my.H.mat(B) %*% my.trans.H3x3(my.H.mat(Hi)) + my.H.mat(C) %*% my.trans.H3x3(my.H.mat(Hj)) + my.H.mat(D) %*% my.trans.H3x3(my.H.mat(Hk))
}
A <- rquat(1)
B <- rquat(1)
C <- rquat(1)
D <- rquat(1)

q <- rquat(1)
q.v <- c(Re(q),i(q),j(q),k(q))
q. <- A*q + B*q*Hi + C*q*Hj + D*q*Hk

Q <- my.H.linear.mat(A,B,C,D)
Q %*% q.v
q.

- 実は、-４ｘ４行列を４パートに分けることができる
  - １ｘ１はスケール
  - １ｘ３は射影
  - ３ｘ１は平行移動
  - ３ｘ３は回転
- となることが知られているのだが、この4x4行列変換に対応する、４つの四元数があって、それは $f(q) = Aq + Bq Hi + Cq Hj + Dq Hk$ という四元数の線形変換式に対応する(という)
- 行列での変換は以下のように

n.pt <-5000
# 単位球面を作ってから、変形・移動
X <- matrix(rnorm(n.pt*3),ncol=3)
X <- X/sqrt(apply(X^2,1,sum))
X[,1] <- X[,1]+X[,2]+ X[,3]^4
X <- X + 1.4
library(rgl)

XX <- (cbind(rep(1,n.pt),X))

X.p <- XX[,2:4]/XX[,1]


s <- 1
pers <-c (0.5,0.4,0.1)
trans <- c(3,3,3)
library(GPArotation)
rot <- Random.Start(3)

M <- matrix(0,4,4)
M[1,1] <- s
M[1,2:4] <- pers
M[2:4,1] <- trans
M[2:4,2:4] <- rot

postXX <- t(M %*% t(XX))

post.X.p <- postXX[,2:4]/postXX[,1]

plot3d(rbind(X.p,rep(max(c(X.p,post.X.p)),3),rep(min(c(X.p,post.X.p)),3)))
points3d(post.X.p,col=2)