Sphere inversionという３次元変形

• こちらに３次元に広がる曲面の変形の話がある。曲面の変形にあたって、Willmore energyを変えない変形をする話がある
• 「うまい感じ」の変形のために、制約を加える話があり、そのうちの一つの制約が"sphere inversion"をさせない、というものである
• ここでいうsphere inversionというのがなんなのかわからなかったので調べもの
• ２次元平面に(単位)円を置いて、その円周で反転するものがある、circleのinversionという。複素平面と考えるとメビウス変換とつながる
• これの３次元版
• ただし、原点を中心とした球面に関する反転ではなく
• 単位球面で反転した後、もう一つ別の単位球を置きそこに関して、もう一度球面反転をする
• このようにすると、原点中心の単位球面の反転を２回続けて行うことで、元に戻すこともできるし、２度目の単位球の中心を原点から外すことで、変形を２度目の反転の中心位置座標を用いてパラメタライズすることができる
• 原点中心での単位球での１度の反転はである
• この計算は四元数の虚部を３次元座標に対応付けることでとすることで簡単に実現できるので
• Sphere inversionはのように表せる。ただし、fは純虚四元数、cは純虚四元数で表された、２つ目の球面反転の中心

library(onion)
library(rgl)
my.plot.quaternion <- function(q,...){
	x <- i(q)
	y <- j(q)
	z <- k(q)
	plot3d(x,y,z,...)
}


my.inverse <- function(q,ctr){
	1/(1/q-ctr)
	#1/q
}
my.imH2R <-function(q){
	t(as.matrix(Im(q)))[,2:4]
}


q <- rnorm(1) *Hi + rnorm(1) *Hj + rnorm(1) * Hk
ctr <- rnorm(1) *Hi + rnorm(1) *Hj + rnorm(1) * Hk




bunnyQ <- bunny[,1]*Hi+bunny[,2]*Hj+bunny[,3]*Hk
bunnyQ <- bunnyQ 
plot3d(my.imH2R(1/bunnyQ))
points3d(my.imH2R(bunnyQ),col=2)


plot3d(my.imH2R(bunnyQ))
points3d(my.imH2R(bunnyQ.),col=2)

n <- 1000
q1 <- rnorm(n) * Hi + rnorm(n) * Hj + rnorm(n) *Hk
for(i in 1:n){
	q1[i] <- 1*q1[i]/Mod(q1[i])
}
ctr <- rnorm(1) * Hi + rnorm(1) * Hj + rnorm(1) *Hk

q1 <- q1
q2 <- my.inverse(q1,ctr)
plot3d(my.imH2R(q2))
points3d(my.imH2R(q2),col=2)

• さて。
• 今、原点から出発する３次元曲線を考える。この３次元曲線の座標をとする。このがsphere inversionの時刻tにおける中心とする
• このとき、時刻0では中心が原点なので無変換(動かない)だが、中心が移動するに従い、３次元オブジェクトが少しずつ変形していく
• ここでとする
• このとき、曲面の各所の時間変形〜各所の接平面に起きる回転を表すの時間微分は,ただしが各所の座標を表すマッピング関数になるという
• sphere inversionでは、球面は(その粗密に変化はしょうじるものの)、球面に移される(らしい)

libarary(devtools)
install_packages("ryamada22/spinSpherization")
library(spinSpherization)
n.psi <- 25
sp.mesh <- my.sphere.tri.mesh(n.psi)
library(rgl)
plot3d(sp.mesh$xyz)
segments3d(sp.mesh$xyz[c(t(sp.mesh$edge)),])

n.step <- 100
ctrs <- matrix(rnorm(3*n.step)*0.01,ncol=3)
ctrs[1,] <- 0
ctrs.cumsum <- apply(ctrs,2,cumsum)
plot3d(ctrs.cumsum,type="l")

my.inverse <- function(q,ctr){
	1/(1/q-ctr)
	#1/q
}

my.inverse.real <- function(x,p){
	q <- x[,1]*Hi+x[,2]*Hj+x[,3]*Hk
	ctr <- p[1]*Hi+p[2]*Hj+p[3]*Hk
	tmp <- my.inverse(q,ctr)
	return(cbind(i(tmp),j(tmp),k(tmp)))
}


meshs <- list()
meshs[[1]] <- sp.mesh$xyz * 3
meshs[[1]][,1] <- meshs[[1]][,1]+3
meshs[[1]][,2] <- meshs[[1]][,2]+10
meshs[[1]][,3] <- meshs[[1]][,3]+5
for(i in 2:n.step){
	meshs[[i]] <- my.inverse.real(meshs[[i-1]],ctrs.cumsum[i,])
	plot3d(meshs[[i]])
	segments3d(meshs[[i]][c(t(sp.mesh$edge)),])

}

meshs <- list()
library(onion)
meshs[[1]] <- bunny
#meshs[[1]][,1] <- meshs[[1]][,1]+3
for(i in 2:n.step){
	meshs[[i]] <- my.inverse.real(meshs[[i-1]],ctrs.cumsum[i,])
	plot3d(meshs[[i]])
	#segments3d(meshs[[i]][c(t(sp.mesh$edge)),])

}