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

  • 3Dでは純虚部分を空間3次元に、4Dでは実+虚4成分を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)
    • 回転は、実は2つのことなるreflectionの連結。\eta_2 ( \eta_1 p \eta_1) \eta_2 = \eta_2 \eta_1 p (\bar{\eta2 \eta1})だから
      • 二つの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 斜めに傾ける。二次元で言えば\begin{pmatrix}1 & k \\ 0 &1 \end{pmatrix}的な
      • 二種類あって・・・
      • Axial-shear
        • 正規直交基底3ベクトルを使って変換する。斜めに伸ばしつつ、回転的成分が入る。斜め関係は、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)
# 二つのベクトルの先を結んだベクトルの長さも保たれている
# 結局3角形の三辺の長さが保たれているから、三角形が合同な三角形に移っていることになり、角も保存
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という面に関する回転だという
  • 球面座標
    • 3次元単位球面座標を純虚四元数a,bで表すとき、abを含む平面は球の大円を定める。この平面に垂直な単位ベクトル(aから反時計回りにbに向かう)mと、回転角、\psiとには\lambda=b a^{-1} = \cos{\psi} + m \sin{\psi}として計算される\lambdaが取れる
    • この\lambdaは大円上の弧に対応する四元数であるが、これは、「ある大円」に対して、その相対的位置に関係なく、弧の長さで決まる数である。
    • また、球面上の3点がa,b,c(a->b->c)で表されているとき、弧a->cの長さはa->b + b->cと「迂回」すれば長くなるが、知りたいのはa->cを通る大円上の距離であるとすると、これを\lambdaを使って\lambda{a,c} = \lambda{a,b}\lambda{b,c}として計算できる
      • その結果、球面上の三角形a->b->c->aとぐるりと戻ると、その距離は\lambda=1となるが、実部が1なので\cos{\psi}\psiは0)となる。

# 球面座標

# オイラー角で単位球面上の2点を決める
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)
  • 射影幾何・同次座標
    • 3次元用の同次座標を四元数の4つの係数に用いる
    • 同次座標はX,Y,Z,WX/W,Y/W,Z/W,1と見るが、このWを四元数の実部、X,Y,Zを虚部とする
    • 無限遠点は実部が0の場合。それは、空間上の点、というよりは、作用方向としての向きを表しているとみてもよい
    • 2点を結ぶ直線上の点をパラメタ扱いするときに、通常なら、パラメタの増加の具合をいじらないといけないけれど、2点の実部を変えることによって、2点間の動きに傾斜を発生させることができる
    • 4元数の線形な変換は、同次座標の線形変換をもたらし、それは3次元的には、複雑な変換に対応してくれる(行列で3次元のときに4x4行列を使うことで、並行移動(第4列)やその他の共形変換が表せることに対応する)
  • 連立四元数線型方程式
    • 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は四元数
# 4元数のベクトル表示から行列表示を作る
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
}

# 4元数行列の2:4 x 2:4 部分の転置
my.trans.H3x3 <- function(H){
	ret <- H
	ret[2:4,2:4] <- t(H[2:4,2:4])
	ret
}
# 4つの四元数から行列を作る
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行列を4パートに分けることができる
      • 1x1はスケール
      • 1x3は射影
      • 3x1は平行移動
      • 3x3は回転
    • となることが知られているのだが、この4x4行列変換に対応する、4つの四元数があって、それは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)