- こちらに記事がある
- 四元数 を実部スカラーと、虚部ベクトルに分ける
- このときだという
- Rのonionパッケージには、四元数のログを取る関数log()があるので、四元数qの
library(onion)
my.exp.q <- function(q,t){
qt <- q *t
a <- Re(qt)
v <- Im(qt)
v.len <- Mod(v)
if(v.len==0){
return(exp(a))
}else{
return(exp(a) * (cos(v.len) + v/v.len * sin(v.len)))
}
}
q <- rnorm(1) + rnorm(1) * Hi + rnorm(1) * Hj + rnorm(1) * Hk
q
logq <- log(q)
my.exp.q(logq,0)
my.exp.q(logq,1)
my.exp.q(logq,2)
q * q
my.exp.q(logq,-1)
q * my.exp.q(logq,-1)
> q
[1]
Re -1.0920124
i -1.5109492
j 0.6109826
k 0.1439245
> logq <- log(q)
> my.exp.q(logq,0)
[1] 1
> my.exp.q(logq,1)
[1]
Re -1.0920124
i -1.5109492
j 0.6109826
k 0.1439245
> my.exp.q(logq,2)
[1]
Re -1.4844904
i 3.2999504
j -1.3344011
k -0.3143346
> q * q
[1]
Re -1.4844904
i 3.2999504
j -1.3344011
k -0.3143346
> my.exp.q(logq,-1)
[1]
Re -0.28221221
i 0.39047938
j -0.15789816
k -0.03719486
> q * my.exp.q(logq,-1)
[1]
Re 1.000000e+00
i -1.734723e-16
j 4.163336e-17
k -3.469447e-17
>
library(onion)
library(rgl)
library(Ronlyryamada)
require(RFOC)
my.exp.q <- function(q,t){
qt <- q *t
a <- Re(qt)
v <- Im(qt)
v.len <- Mod(v)
v.len0 <- which(v.len==0)
ret <- exp(a) * (cos(v.len) + v/v.len * sin(v.len))
if(length(v.len0)>0){
ret[v.len0] <- exp(a[v.len0])
}
return(ret)
}
ts <- seq(from=-5,to=5,length=2000)
A <- rnorm(4) * 2
q <- A[1] + A[2] * Hi + A[3] * Hj + A[4] * Hk
X0 <- matrix(rnorm(3),nrow=1)
X0.q <- Hi * X0[1] + Hj * X0[2] + Hk * X0[3]
Xs <- matrix(0,length(ts),3)
for(i in 1:length(ts)){
tmp.q <- my.exp.q(q,ts[i])
Xs[i,] <- as.matrix(Conj(tmp.q) * X0.q * tmp.q)[2:4]
}
plot3d(Xs,type="l")
- 四元数の自由度は4
- もう少し自由度が高いこともあるだろう
- 3次元空間のある点を基点とした正規直行基底が、その起点を移動し、さらに回転し、全体として拡縮するとすると
- 起点の移動分で自由度が3、3次元空間での回転は自由度が2、拡縮係数で自由度が1の、併せて6
- こんな移動だと、平行移動のために、同次座標を使うことにして、4x4行列を作って、その指数関数にするという手もある
library(onion)
library(rgl)
library(Ronlyryamada)
require(RFOC)
my.exp.q <- function(q,t){
qt <- q *t
a <- Re(qt)
v <- Im(qt)
v.len <- Mod(v)
v.len0 <- which(v.len==0)
ret <- exp(a) * (cos(v.len) + v/v.len * sin(v.len))
if(length(v.len0)>0){
ret[v.len0] <- exp(a[v.len0])
}
return(ret)
}
ts <- seq(from=-5,to=5,length=2000)
A <- rnorm(4) * 2
q <- A[1] + A[2] * Hi + A[3] * Hj + A[4] * Hk
X0 <- matrix(rnorm(3),nrow=1)
X0.q <- Hi * X0[1] + Hj * X0[2] + Hk * X0[3]
Xs <- matrix(0,length(ts),3)
for(i in 1:length(ts)){
tmp.q <- my.exp.q(q,ts[i])
Xs[i,] <- as.matrix(Conj(tmp.q) * X0.q * tmp.q)[2:4]
}
plot3d(Xs,type="l")
my.rot.3d <- function(theta,phi){
M1 <- M2 <- diag(rep(1,3))
M1[1:2,1:2] <- matrix(c(cos(theta),-sin(theta),sin(theta),cos(theta)),byrow=TRUE,2,2)
M2[2:3,2:3] <- matrix(c(cos(phi),-sin(phi),sin(phi),cos(phi)),byrow=TRUE,2,2)
return(M2 %*% M1)
}
n <- my.rot.3d(runif(1)*2*pi,runif(1)*2*pi)
x <- n[,3]
X <- rnorm(3) * 3
k <- runif(1)*3
N <- k * my.rot.3d(runif(1)*2*pi,runif(1)*2*pi)
my.trans.matrix <- function(x,X,n,N){
P <- Q <- R <- diag(rep(1,4))
P[1:3,4] <- X
R[1:3,4] <- -x
Q[1:3,1:3] <- N %*% solve(n)
return(P %*% Q %*% R)
}
exp.m2 <- function(A,n){
eigen.out<-eigen(A)
V<-eigen.out[[2]]
U<-solve(V)
B<-diag(exp(log(eigen.out[[1]])*n))
X <- V%*%B%*%U
return(list(matrix = X, eigen.vs <- eigen.out[[1]]))
}
trans.out <- my.trans.matrix(x,X,n,N)
trans.out %*% c(x,1)
X
trans.out %*% rbind(n+x,rep(1,3))
X + N
xXnN <- rbind(x,X)
for(i in 1:3){
xXnN <- rbind(xXnN,x+n[,i],X+N[,i])
}
xXnN <- rbind(xXnN,rep(min(xXnN),3),rep(max(xXnN),3))
plot3d(xXnN)
for(i in 1:3){
segments3d(rbind(x,x+n[,i]),col=2)
segments3d(rbind(X,X+N[,i]),col=3)
}
ts <- seq(from=0,to=1,length=30)
for(i in 1:length(ts)){
tmp.m <- Re(exp.m2(trans.out,ts[i])[[1]])
tmp.ori <- matrix(0,4,4)
tmp.ori[1:3,1:4] <- cbind(x,n+x)
tmp.ori[4,] <- 1
tmp.post <- tmp.m %*% tmp.ori
tmp.post.2 <- tmp.post[1:3,]
tmp.X <- tmp.post.2[,1]
tmp.N.tip <- tmp.post.2[,2:4]
for(ii in 1:3){
segments3d(rbind(tmp.X,tmp.N.tip[,ii]))
}
}