遺伝学・遺伝統計学関連の姉妹ブログ『ryamadaの遺伝学・遺伝統計学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

Rで曲率指定的な曲面共形変換変形

R 疎行列 四元数 曲率 DEC 離散微分幾何
  • ここ数日(数週間？)C++で書かれていて、MATLABにはすぐ連結する(らしい)疎行列ライブラリSuiteSparseを用いて離散微分幾何的アプローチの曲面変形について、調べてきた
  • RにもSuiteSparse準拠の疎行列用パッケージMatrixがあるので、それに連結してみようといろいろやってきた
  • それなりに(遅いけれど)動く模様なので、まずは、図(黒い正球を曲率濃淡ランダム二次元関数に基づいて変形した赤い閉曲面)のようなものが描けますよ、というソース
  • ごちゃごちゃしているので、このソースをRmd化しておくことにする(けれど、それは、ちょっと後)

my.readOBJ <- 
function (con, ...) 
{
    lines <- readLines(con)
    instrs <- sub(" .*", "", lines)
    vertices <- read.table(textConnection(lines[instrs == "v"]), 
        col.names = c("instr", "x", "y", "z"), colClasses = c(instr = "character", 
            x = "numeric", y = "numeric", z = "numeric"))
    vertices <- with(vertices, rbind(x, y, z))
    normals <- c()
    normals <- read.table(textConnection(lines[instrs == "vn"]), 
        col.names = c("instr", "x", "y", "z"), colClasses = c(instr = "character", 
            x = "numeric", y = "numeric", z = "numeric"))
    normals <- with(normals, rbind(x, y, z))
    textures <- c()
    textures <- read.table(textConnection(lines[instrs == "vt"]), 
        col.names = c("instr", "x", "y"), colClasses = c(instr = "character", 
            x = "numeric", y = "numeric"))
    textures <- with(textures, rbind(x, y))

    tfaces <- grepl("^f\\W+\\w+\\W+\\w+\\W+\\w+$", lines)
    triangles <- read.table(textConnection(lines[instrs == "f"]), col.names = c("instr", 
        "v1", "v2", "v3"), colClasses = "character")
    
   # triangles <- read.table(textConnection(lines[tfaces]), col.names = c("instr", 
    #    "v1", "v2", "v3"), colClasses = "character")
    triangles <- with(triangles, rbind(v1, v2, v3))
    triangles1 <- triangles2 <- triangles3 <- matrix(0,3,length(triangles[1,]))
    for(i in 1:length(triangles[1,])){
			for(j in 1:3){
				tmp <- unlist(strsplit(triangles[j,i],"/"))
				triangles1[j,i] <- as.numeric(tmp[1])
				triangles2[j,i] <- as.numeric(tmp[2])
				triangles3[j,i] <- as.numeric(tmp[3])
			}
		}
		ret <- list(vertices=vertices,normals=normals,textures=textures,triangles1=triangles1,triangles2=triangles2,triangles3=triangles3)
		return(ret)
}

con <- "sphere_original.obj"
obj.data <- my.readOBJ(con)
m <- tmesh3d(obj.data$vertices,obj.data$triangles1,homogeneous=FALSE)
shade3d(m)
library(onion)
vertices <- obj.data$vertices[1,] * Hi + obj.data$vertices[2,] * Hj + obj.data$vertices[3,] * Hk
uv <- obj.data$textures
faces.v <- obj.data$triangles1
faces.uv <- obj.data$triangles2

my.rho <- function(x,y,vs,Vs,Vs2){
	X <- x
	Y <- y
	if(Y < 0.5){
		Y <- y*2
	}else{
		Y <- 2*(y-0.5) 
	}
	X <- (X-0.5)*2
	Y <- (Y-0.5)*2
	if(X^2+Y^2>1){
		return(0)
	}
	if(y < 0.5){
		Y <- -Y
		Z <- sqrt(1-(X^2+Y^2))
	}else{
		Z <- -sqrt(1-(X^2+Y^2))
	}
	xyz <- c(X,Y,Z)
	tmp <- (t(vs) %*% xyz)

	for(i in 1:length(tmp)){
		if(tmp[i] > 1) tmp[i] <- 1
		if(tmp[i] < -1) tmp[i] <- -1
	}
	tmp <- acos(tmp)
	ret <- sum(Vs*exp(-Vs2*tmp^2))
	ret
}


x <- runif(30000)
y <- runif(30000)
np <- 4
vs <- matrix(rnorm(np*3),ncol=3)
vs <- vs/sqrt(apply(vs^2,1,sum))
vs <- t(vs)
Vs <- runif(np)*5
Vs2 <- runif(np)*2
rho <- rep(0,length(x))
for(i in 1:length(rho)){
	rho[i] <- my.rho(x[i],y[i],vs,Vs,Vs2)
}
rho[which(rho==0)] <- min(rho[rho!=0])-0.01
rho.st <- (max(rho)-rho)/(max(rho)-min(rho))
plot(x,y,col=gray(rho.st),pch=20,asp=TRUE)


rho <- rep(0,length(faces.uv[1,]))
for(i in 1:length(rho)){
	threes <- faces.uv[,i]
	tmp <- rep(0,3)
	for(j in 1:3){
		tmp[j] <- my.rho(obj.data$textures[1,threes[j]],obj.data$textures[2,threes[j]],vs,Vs,Vs2)
	}
	rho[i] <- mean(tmp)
}
rho <- ((rho-min(rho))/(max(rho)-min(rho)) - 0.5)*2*10
# 頂点にrhoを与え、そこから三角形のrhoを出してもよい
rho.v <- rep(0,length(vertices))
xyz.ori <- t(as.matrix(vertices)[2:4,])
np <- 4
NP <- matrix(rnorm(np*3),ncol=3)
NP <- t(t(NP)/sqrt(apply(NP^2,1,sum)))
peak <- rnorm(np)*50
sds <- runif(np)*20
rr <- sample(1:8,np)
for(i in 1:np){
	tmp <- (xyz.ori %*% NP[i,])
	tmp <- tmp/max(abs(tmp))
	tmp <- acos(tmp)
	tmp2 <- peak[i] * exp(-(tmp^2)^rr[i]*sds[i])
	rho.v <- rho.v + tmp2
}
tmp.rho <- matrix(rho.v[faces.v],ncol=3)
rho <- apply(tmp.rho,1,mean)
rho <- ((rho-min(rho))/(max(rho)-min(rho)) - 0.5)*2*10
rho.st <- (max(rho.v)-rho.v)/(max(rho.v)-min(rho.v))

plot3d(xyz.ori,col=rgb(rho.st,1-rho.st,1))

library(Matrix)
make.E <- function(vertices,faces.v,rho){
	# 三角形の面積
	edge1 <- vertices[faces.v[2,]]-vertices[faces.v[1,]]
	edge2 <- vertices[faces.v[3,]]-vertices[faces.v[1,]]
	tmp <- edge1 * edge2
	A <- abs(i(tmp)+j(tmp)+k(tmp))/2
	# 三角形ごとの計算用係数
	coef.a <- -1/(4*A)
	coef.b <- rho/6
	coef.c <- A*rho^2/9
	
	# Rでは四元数を要素とする行列がないので、re,i,j,kごとに正方行列を作ることにする
	E.re <- E.i <- E.j <- E.k <- sparseVector(c(0),i=c(1),length=length(vertices)^2)
	
	e.q <- list()
	e.q[[1]] <- vertices[faces.v[2,]]-vertices[faces.v[3,]]
	e.q[[2]] <- vertices[faces.v[3,]]-vertices[faces.v[1,]]
	e.q[[3]] <- vertices[faces.v[1,]]-vertices[faces.v[2,]]
	for(i in 1:3){
		for(j in 1:3){
			tmp <- coef.a * e.q[[i]] * e.q[[j]]+ coef.b * (e.q[[j]] -e.q[[i]] ) + coef.c
			addr <- faces.v[i,] + (length(vertices)*(faces.v[j,]-1))
			
			tmp.v <- t(as.matrix(tmp))
			tmp.out <- my.vector.access(tmp.v,addr)
			E.re <- E.re + sparseVector(tmp.out[[2]][,1],tmp.out[[1]],length(vertices)^2)
			E.i <- E.re + sparseVector(tmp.out[[2]][,2],tmp.out[[1]],length(vertices)^2)
			E.j <- E.re + sparseVector(tmp.out[[2]][,3],tmp.out[[1]],length(vertices)^2)
			E.k <- E.re + sparseVector(tmp.out[[2]][,4],tmp.out[[1]],length(vertices)^2)
			ord <- order(addr)
			rle.out <- rle(addr[ord])
			num.row <- length(rle.out$values)
			num.col <- max(rle.out$lengths)
			tmp1 <- rep(1:num.row,rle.out$lengths)
			tmp2 <- c()
			for(k in 1:num.row){
				tmp2 <- c(tmp2,1:rle.out$lengths[k])
			}
			addr2 <- tmp1 + num.row*(tmp2-1)
			tmp.M <- matrix(0,num.row,num.col)
			tmp.M[addr2] <- Re(tmp[ord])
			tmp.sum <- apply(tmp.M,1,sum)
			E.re <- E.re + sparseVector(tmp.sum,rle.out$values,length(vertices)^2)
			tmp.M[addr2] <- i(tmp)
			tmp.sum <- apply(tmp.M,1,sum)
			E.i <- E.i + sparseVector(tmp.sum,rle.out$values,length(vertices)^2)
			tmp.M[addr2] <- j(tmp)
			tmp.sum <- apply(tmp.M,1,sum)
			E.j <- E.j + sparseVector(tmp.sum,rle.out$values,length(vertices)^2)
			tmp.M[addr2] <- k(tmp)
			tmp.sum <- apply(tmp.M,1,sum)
			E.k <- E.k + sparseVector(tmp.sum,rle.out$values,length(vertices)^2)
			#E.re[addr] <- E.re[addr] + sparseVector(Re(tmp),1:length(tmp),length(tmp))
			#E.i[addr] <- E.i[addr] + sparseVector(i(tmp),1:length(tmp),length(tmp))
			#E.j[addr] <- E.j[addr] + sparseVector(j(tmp),1:length(tmp),length(tmp))
			#E.k[addr] <- E.k[addr] + sparseVector(k(tmp),1:length(tmp),length(tmp))
		}
	}

	return(list(E.re=Matrix(E.re,length(vertices),length(vertices)),E.i=Matrix(E.i,length(vertices),length(vertices)),E.j=Matrix(E.j,length(vertices),length(vertices)),E.k=Matrix(E.k,length(vertices),length(vertices))))
}


my.vector.access <- function(v,a,func=sum,zero=0){
	if(is.vector(v)){
		v <- matrix(v,ncol=1)
	}
	ord <- order(a)
	rle.out <- rle(a[ord])
	num.row <- length(rle.out[[1]])
	num.col <- max(rle.out[[1]])
	tmp1 <- rep(1:num.row,rle.out[[1]])
	tmp2 <- c()
	for(i in 1:num.row){
		tmp2 <- c(tmp2,1:rle.out[[1]][i])
	}
	addr <- tmp1 + num.row*(tmp2-1)
	ret.v <- matrix(0,num.row,ncol(v))
	for(i in 1:ncol(v)){
		if(zero==0){
			tmp.V <- sparseVector(v[ord,i],i=addr,length=num.row*num.col)
			M <- Matrix(tmp.V,num.row,num.col)
		}else{
			M <- matrix(zero,num.row,num.col)
			M[addr] <- v[ord,i]

		}
		ret.v[,i] <- apply(M,1,func)

	}
	return(list(A = rle.out[[2]],V = ret.v))
}

v <- 1:10
v <- cbind(1:10,1:10)
a <- banchi <- c(1,2,3,1,2,4,50,1,3,2)
oo <- my.vector.access(v,a)
oo2 <- my.vector.access(v,a,func=prod,zero=1)


make.E <- function(vertices,faces.v,rho){
	# 三角形の面積
	edge1 <- vertices[faces.v[2,]]-vertices[faces.v[1,]]
	edge2 <- vertices[faces.v[3,]]-vertices[faces.v[1,]]
	tmp <- edge1 * edge2
	A <- abs(i(tmp)+j(tmp)+k(tmp))/2
	# 三角形ごとの計算用係数
	coef.a <- -1/(4*A)
	coef.b <- rho/6
	coef.c <- A*rho^2/9
	
	# Rでは四元数を要素とする行列がないので、re,i,j,kごとに正方行列を作ることにする
	E.re <- E.i <- E.j <- E.k <- sparseVector(c(0),i=c(1),length=length(vertices)^2)
	
	e.q <- list()
	e.q[[1]] <- vertices[faces.v[2,]]-vertices[faces.v[3,]]
	e.q[[2]] <- vertices[faces.v[3,]]-vertices[faces.v[1,]]
	e.q[[3]] <- vertices[faces.v[1,]]-vertices[faces.v[2,]]
	for(i in 1:3){
		for(j in 1:3){
			tmp <- coef.a * e.q[[i]] * e.q[[j]]+ coef.b * (e.q[[j]] -e.q[[i]] ) + coef.c
			addr <- faces.v[i,] + (length(vertices)*(faces.v[j,]-1))
			
			tmp.v <- t(as.matrix(tmp))
			tmp.out <- my.vector.access(tmp.v,addr)
			E.re <- E.re + sparseVector(tmp.out[[2]][,1],tmp.out[[1]],length(vertices)^2)
			E.i <- E.re + sparseVector(tmp.out[[2]][,2],tmp.out[[1]],length(vertices)^2)
			E.j <- E.re + sparseVector(tmp.out[[2]][,3],tmp.out[[1]],length(vertices)^2)
			E.k <- E.re + sparseVector(tmp.out[[2]][,4],tmp.out[[1]],length(vertices)^2)
		}
	}

	return(list(E.re=Matrix(E.re,length(vertices),length(vertices)),E.i=Matrix(E.i,length(vertices),length(vertices)),E.j=Matrix(E.j,length(vertices),length(vertices)),E.k=Matrix(E.k,length(vertices),length(vertices))))
}
# 返り値を疎ベクトルのリストにしたもの)
make.E.v <- function(vertices,faces.v,rho){
	# 三角形の面積
	edge1 <- vertices[faces.v[2,]]-vertices[faces.v[1,]]
	edge2 <- vertices[faces.v[3,]]-vertices[faces.v[1,]]
	tmp <- edge1 * edge2
	A <- abs(i(tmp)+j(tmp)+k(tmp))/2
	# 三角形ごとの計算用係数
	coef.a <- -1/(4*A)
	coef.b <- rho/6
	coef.c <- A*rho^2/9
	
	# Rでは四元数を要素とする行列がないので、re,i,j,kごとに正方行列を作ることにする
	E.re <- E.i <- E.j <- E.k <- sparseVector(c(0),i=c(1),length=length(vertices)^2)
	
	e.q <- list()
	e.q[[1]] <- vertices[faces.v[2,]]-vertices[faces.v[3,]]
	e.q[[2]] <- vertices[faces.v[3,]]-vertices[faces.v[1,]]
	e.q[[3]] <- vertices[faces.v[1,]]-vertices[faces.v[2,]]
	for(i in 1:3){
		for(j in 1:3){
			tmp <- coef.a * e.q[[i]] * e.q[[j]]+ coef.b * (e.q[[j]] -e.q[[i]] ) + coef.c
			addr <- faces.v[i,] + (length(vertices)*(faces.v[j,]-1))
			
			tmp.v <- t(as.matrix(tmp))
			tmp.out <- my.vector.access(tmp.v,addr)
			E.re <- E.re + sparseVector(tmp.out[[2]][,1],tmp.out[[1]],length(vertices)^2)
			E.i <- E.re + sparseVector(tmp.out[[2]][,2],tmp.out[[1]],length(vertices)^2)
			E.j <- E.re + sparseVector(tmp.out[[2]][,3],tmp.out[[1]],length(vertices)^2)
			E.k <- E.re + sparseVector(tmp.out[[2]][,4],tmp.out[[1]],length(vertices)^2)
		}
	}

	return(list(E.re=E.re,E.i=E.i,E.j=E.j,E.k=E.k))
}

# 四元数対応実行列作成

my.qMtorM <- function(Es){
	n <- sqrt(length(Es[[1]]))
	N <- (n*4)^2
	init.id <- c(1:4,(1:4)+n*4,(1:4)+n*4*2,(1:4)+n*4*3)
	spacing.id <- c(outer((0:(n-1)*4),n*4*4*(0:(n-1)),"+"))
	ret <- sparseVector(c(0),i=c(1),N)
	a <- c(1,2,3,4,2,1,4,3,3,4,1,2,4,3,2,1)
	b <- c(1,1,1,1,-1,1,1,-1,-1,-1,1,1,-1,1,-1,1)
	for(j in 1:length(a)){
		tmp.v <- sparseVector(b[j] * Es[[a[j]]]@x,i=init.id[j]+spacing.id[Es[[a[j]]]@i],length=N)
		ret <- ret + tmp.v
	}
	Matrix(ret,n*4,n*4)
}

n <- 100
p <- 10
s <- sample(1:n,p)
h <- rquat(p)
tmp.Es.v <- list(sparseVector(Re(h),s,n),sparseVector(i(h),s,n),sparseVector(j(h),s,n),sparseVector(k(h),s,n))
N <- 10
n <- N^2
p <- 80
s <- sample(1:n,p)
h <- rquat(p)
tmp.Es.v <- list(sparseVector(Re(h),s,n),sparseVector(i(h),s,n),sparseVector(j(h),s,n),sparseVector(k(h),s,n))

tmp.Es.real <- my.qMtorM(tmp.Es.v)

image(tmp.Es.real)

my.qVtorV <- function(v){
	c(as.matrix(v))
}
tmp.v <- rquat(10)
tmp.v
my.qVtorV(tmp.v)

# 四元数行列(ベクトルリスト形式)Mと
# 四元数ベクトルbとの
# Mx = b の解

my.linear.solver <- function(Es,b){
	b.re <- my.qVtorV(b)
	E.re <- my.qMtorM(Es)
	solve(E.re,b.re)
}
N <- 10
n <- N^2
p <- 80
s <- sample(1:n,p)
h <- rquat(p)
tmp.Es.v <- list(sparseVector(Re(h),s,n),sparseVector(i(h),s,n),sparseVector(j(h),s,n),sparseVector(k(h),s,n))
tmp.v <- rquat(10)

my.linear.solver(tmp.Es.v,tmp.v)

b <- rep(1+1*Hi+1*Hj+1*Hk,nrow(Es[[1]]))
tmp.out <- my.linear.solver(Es.v,b)

plot(tmp.out)

# Aは実正方行列
my.inv.pow <- function(A,n.iter=3,b=rep(1,ncol(A)),log=FALSE){
	x <- b
	if(log){
		x.log <- matrix(0,n.iter+1,ncol(A))
		x.log[1,] <- x
	}
	#x <- x/sqrt(sum(x^2))
	A. <- solve(A)
	for(i in 1:n.iter){
		x <- A. %*% x
		x <- x/sqrt(sum(x^2))
		if(log){
			x.log[i+1,] <- x
		}
		
	}
	if(log){
		return(list(x=x,x.log=x.log))
	}else{
		return(list(x=x,x.log=matrix(0,0,ncol(A))))
	}
}

my.inv.pow.2 <- function(A,n.iter=3,b=rep(1,ncol(A)),log=FALSE){
	x <- b
	if(log){
		x.log <- matrix(0,n.iter+1,ncol(A))
		x.log[1,] <- x
	}
	#x <- x/sqrt(sum(x^2))
	#A. <- solve(A)
	for(i in 1:n.iter){
		x2 <- solve(A,x)
		x <- x2/sqrt(sum(x2^2))
		if(log){
			x.log[i+1,] <- x
		}
		
	}
	if(log){
		return(list(x=x,x.log=x.log))
	}else{
		return(list(x=x,x.log=matrix(0,0,ncol(A))))
	}
}




tmp.Es.v <- list()
for(i in 1:4){
	tmp.Es.v[[i]] <- sparseVector(rnorm(16),i=1:16,length=16)
}

N <- 20
X <- matrix(rnorm(N^2),N,N)
b <- rnorm(N)
tmp.outs <- my.inv.pow(X,n.iter=100,b=b,log=TRUE)
tmp.outs.2 <- my.inv.pow(X,n.iter=100,b=b,log=TRUE)
matplot(tmp.outs[[2]],type="l")

#Es <- make.E(vertices,faces.v,rho)
Es.v <- make.E.v(vertices,faces.v,rho)
E.re <- my.qMtorM(Es.v)
lambda <- my.inv.pow.2(E.re,n.iter=3)[[1]]
plot(lambda)

my.make.L <- function(vertices,faces.v){
	n.v <- length(vertices)
	L <- sparseVector(c(0),i=c(1),length=n.v^2)
	for(i in 1:3){
		v.ord <- ((1:3)+i+1) %% 3 + 1
		k1 <- faces.v[v.ord[1],]
		k2 <- faces.v[v.ord[2],]
		k3 <- faces.v[v.ord[3],]
		# 頂点四元数
		v1 <- vertices[k1]
		v2 <- vertices[k2]
		v3 <- vertices[k3]
		
		# edge 四元数
		u1 <- v2-v1
		u2 <- v3-v1
		# edge 四元数(純虚四元数)の積は実部がドット積、虚部がクロス積ベクトル
		u12 <- u1 * u2
		cotAlpha <- (-Re(u12))/Mod(Im(u12))
		# このcotAlphaを行列Lの[k2,k2],[k3,k3],[k2,k3],[k3,k2]に加算する
		# 疎ベクトルで格納する
		addrk2k2 <- k2 + (k2-1)*n.v
		addrk3k3 <- k3 + (k3-1)*n.v
		addrk2k3 <- k2 + (k3-1)*n.v
		addrk3k2 <- k3 + (k2-1)*n.v
		
		addr <- c(addrk2k2,addrk3k3,addrk2k3,addrk3k2)
		
		val <- c(cotAlpha,cotAlpha,-cotAlpha,-cotAlpha)/2
		
		tmp.out <- my.vector.access(val,addr)
		L <- L + sparseVector(tmp.out[[2]][,1],tmp.out[[1]],n.v^2)
	}
	L
}
L <- my.make.L(vertices,faces.v)

my.make.omega <- function(vertices,faces.v,lambda){
	n.v <- length(vertices)
	omega <- rep(0*Hi,n.v)
	for(i in 1:3){
		v.ord <- ((1:3)+i+1) %% 3 + 1
		k1 <- faces.v[v.ord[1],]
		# 対向辺の向きは頂点IDの大小順にそろえる
		k23 <- rbind(faces.v[v.ord[2],],faces.v[v.ord[3],])
		k23 <- apply(k23,2,sort)
		k2 <- k23[2,]
		k3 <- k23[1,]
		# 頂点四元数
		v1 <- vertices[k1]
		v2 <- vertices[k2]
		v3 <- vertices[k3]
		
		edge <- v3-v2
		
		# lambdaの四元数化、とその共役四元数化
		lambda.mat <- matrix(lambda,nrow=4)
		lambda.q <- as.quaternion(lambda.mat)
		lambda.mat.2 <- lambda.mat
		lambda.mat.2[2:4,] <- -lambda.mat.2[2:4,]
		lambda.q. <- as.quaternion(lambda.mat.2)
		lambda1 <- lambda.q[k2]
		lambda1. <- lambda.q.[k2]
		lambda2 <- lambda.q[k3]
		lambda2. <- lambda.q.[k3]
		
		val <- 1/3 * lambda1. * edge * lambda1 + 1/6 * lambda1. * edge * lambda2 + 1/6 * lambda2. * edge * lambda1 + 1/3 * lambda2. * edge * lambda2

		# edge 四元数
		u1 <- v2-v1
		u2 <- v3-v1
		# edge 四元数(純虚四元数)の積は実部がドット積、虚部がクロス積ベクトル
		u12 <- u1 * u2
		cotAlpha <- (-Re(u12))/Mod(Im(u12))

		Val <- cotAlpha * val /2
		Val.m <- t(as.matrix(Val))
		tmp.out2 <- my.vector.access(-Val.m,k2)
		tmp.out3 <- my.vector.access(Val.m,k3)
		omega[tmp.out2[[1]]] <- omega[tmp.out2[[1]]] + as.quaternion(t(tmp.out2[[2]]))
		omega[tmp.out3[[1]]] <- omega[tmp.out3[[1]]] + as.quaternion(t(tmp.out3[[2]]))
	}
	omega.re <- as.matrix(omega)
	omega.re <- omega.re-apply(omega.re,1,mean)
	c(omega.re)
}

# lambdaを四元数化
lambda.q <- as.quaternion(matrix(lambda,nrow=4))
# 点ごとに四元数を与え、それを点の数ｘ４の長さのベクトルにしたもの
omega <- my.make.omega(vertices,faces.v,lambda)

L.q <- list()
L.q[[1]] <- L
for(i in 2:4){
	L.q[[i]] <- L*0
}
L.re <- my.qMtorM(L.q)
new.vertices <- solve(L.re,omega)
new.vertices <- as.quaternion(matrix(new.vertices,nrow=4))


newVertices <- matrix(new.vertices,nrow=4)[2:4,]

xyz.ori <- t(as.matrix(vertices)[2:4,])
xyz.new <- t(as.matrix(new.vertices)[2:4,])

mean.ori <- apply(xyz.ori,2,mean)
xyz.ori. <- t(t(xyz.ori)-mean.ori)

mean.new <- apply(xyz.new,2,mean)
xyz.new. <- t(t(xyz.new)-mean.new)


max.ori. <- max(abs(xyz.ori.))
xyz.ori.st <- xyz.ori./max.ori.
max.new. <- max(abs(xyz.new.))
xyz.new.st <- xyz.new./max.new.

plot3d(rbind(xyz.ori.st,xyz.new.st),col=rep(1:2,each=length(xyz.ori[,1])))


plot(c(xyz.ori),c(xyz.new))


library(rgl)
plot3d(t(as.matrix(vertices)[2:4,]))
open3d()
plot3d(t(as.matrix(new.vertices)[2:4,]))