四元数の指数関数

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
> 
  • 四元数の指数関数による3次元空間軌道
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
#q <- 1+ 3* Hi
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]
	#Xs[i,] <- rotate(X0,my.exp.q(q,ts[i]))
}

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
#q <- 1+ 3* Hi
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]
	#Xs[i,] <- rotate(X0,my.exp.q(q,ts[i]))
}

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)
	# P=V,P^{-1}=U
	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]))
	}
}