2015-04-08
Willmoreエネルギーと閉曲面の正球化

昨日の記事は、正球三角形メッシュを適当にまげて共形変換的に曲面を作る話だった。そのときの「曲げ」の情報は、各点の属性である『平均曲率mean curvature』の増減量としてrhoを指定する、というやり方だった
今日のは逆で、ある曲面三角形メッシュを共形変換して正球マップにしようというもの
ネタは同じくkeenan DECサイトのこちらのペイパー
ソースの元としてはSpinTransformation
曲面を正球にするっていうのは、「曲率をいたるところ同じ」にすること
曲面積を最小化することでそれを実現する方法と、平均曲率を揃える(その分散を最小化する)という方法に２大別できる
平均曲率の方を使えば、曲面上で各所の平均曲率の二乗を積分してやったものを最小化することになるが、これがWillmore energy(Willmore energy)
上掲ペイパーの方法は、点の座標を変えるのではなく、点の平均曲率を変える。ただし、変えながら、「変えた結果がきちんと曲面として成立する」ことの担保をしておく。平均曲率の二乗の和がWillmore energyなので、それを減らすにあたって、曲率の山から、曲率が低い方に「曲率を拡散させる」。このように減衰は「曲率が高いところほど大きく減らす」作戦。 $u^2$ で評価しているから $-2u$ の勢いで減らす、というような単純なもの( $u^{k+1} = u^k - 2\tau u^k$ のように、単位時間 $\tau$ 後の値を時間発展させる
ステップで表せば
- 平均曲率を算出(ただし、fを埋め込みとすると、平均曲率Hは、曲面の法線ベクトルNを使って $\bigtriangleup f = 2 HN$ (法線単位ベクトルとの内積としてスカラー量が決まる)となる)
- 曲率スカラー場における、gradientを算出
- 正規直交基底のグラム-シュミット算出
- 各所での曲率の平坦化時間変化量を割り付け算出
- 適当な時間幅で変化
- 変化後の座標を共形変換計算(曲面変形(昨日の記事)と同じ)
このステップを何度か繰り返す
この方法でよいのなら、曲率の拡散の部分は、「グラフ上の拡散」でも良さそう…。曲率勾配に応じて各辺の移動測度が決まって…みたいな。
ちょっと寄り道。適当な時間幅で変化、というときに、この時間幅(離散時間変化の幅)には、「うまく行く」範囲と「おかしくなる」範囲がある。うまく行く範囲は生物の形状変化として、「形をおおよそ保持した変化」に相当し、おかしくなる範囲は、「形をシャッフルする変化」に相当するかもしれない
ただし、fを埋め込みとすると、平均曲率 $\bigtriangleup f = 2 HN$ であり、
cotangent の元論文
正球に戻すだけなら
- 曲率スカラー場における、gradientを算出
- 適当な時間幅で変化
- 変化後の座標を共形変換計算(曲面変形(昨日の記事)と同じ)
このステップを何度か繰り返す
だけでよい
それをepub化
作者: ryamada
発売日: 2015/04/03
メディア: Kindle版
この商品を含むブログ (1件) を見る
---
title: "球面の変形〜Rで学ぶ曲面の四元数共形変換"
author: "ryamada"
date: "Saturday, April 04, 2015"
output: html_document
---
# 球面の変形〜Rで学ぶ曲面の四元数共形変換〜

# はじめに

この文書は球面とそれを滑らかに変形してできる曲面(球的閉曲面)をコンピュータ的に扱う方法に関するものである。

次の順序で構成する。

-１．Rについての知識があることを前提として、球的閉曲面の操作の概要を述べる

-２．その操作をする関数について解説する

-３．実例を扱う

-４．最後にその関数がうまく機能する背景知識について概説する

# Rで球的閉曲面メッシュを扱うこと

## 使用するパッケージ一覧

```{r}
library(onion) # 四元数
library(geometry) # 3次元ベクトル・クロス積関数
library(Matrix) # 疎行列
library(rgl) # 3Dプロット
library(knitr) # 文書作成・画像埋め込み
library(igraph) # グラフ理論
```

## 球面の三角形メッシュを作る
北極から南極まで、等緯度間隔で球を輪切りにする。
各輪切りには、等経度間隔の点を取る。点間距離が、緯度経度の全体で同じくらいになるように、配置数を調節する。

同一緯線上の点には周回するように辺を置き、
隣接緯線上の点には、橋渡しの辺が経線方向に近くなるように配置し、三角化する。

引数
n.psiは単位球の緯度の刻み数。北極-南極の両端を含めてn.psi段

出力
xyz 頂点３次元座標
edge 三角形の辺の両端頂点ID ２列の行列
triangles 三角形の３頂点ID ３列の行列
n.psi 緯度の数
thetas 緯度ごとに配した点の経度
n.triangles 隣接する緯度の間に作成された三角形の個数


以下に実装した関数を示し、それを用いた、球面三角形メッシュと、三角形タイル敷きつめの図を表示する。

```{r}
my.sphere.tri.mesh <- function(n.psi=30){
  thetas <- list()
	psis <- seq(from=-pi/2,to=pi/2,length=n.psi)
	d.psis <- psis[2]-psis[1]
	hs <- sin(psis)
	rs <- sqrt(1-hs^2)
	ls <- 2*pi*rs
	n.thetas <- floor(ls/d.psis)
	thetas[[1]] <- c(2*pi)
	for(i in 2:(n.psi-1)){
		thetas[[i]] <- seq(from=0,to=2*pi,length=n.thetas[i]+1)
		thetas[[i]] <- thetas[[i]][-(n.thetas[i]+1)]
	}
	thetas[[n.psi]] <- c(2*pi)
	sapply(thetas,length)

	bridge <- list()
	for(i in 1:(n.psi-1)){
		a <- c(thetas[[i]],2*pi)
		b <- c(thetas[[i+1]],2*pi)
		bridge[[i]] <- matrix(c(1,1),1,2)
		loop <- TRUE
		while(loop){
			n.r <- nrow(bridge[[i]])
			id.a <- bridge[[i]][n.r,1] + 1
			id.b <- bridge[[i]][n.r,2] + 1
			if(id.a > length(thetas[[i]]) & id.b > length(thetas[[i+1]])){
				if(id.a-1!=1 & id.b-1!=1){
					bridge[[i]] <- rbind(bridge[[i]],c(1,id.b-1))
				}
				
				loop <- FALSE
			}else{
				if(id.a > length(thetas[[i]])){
					tmp <- c(id.a-1,id.b)
				}else if(id.b > length(thetas[[i+1]])){
					tmp <- c(id.a,id.b-1)
				}else{
					if(a[id.a] < b[id.b]){
						tmp <- c(id.a,id.b-1)
					}else{
						tmp <- c(id.a-1,id.b)
					}
				}
				bridge[[i]] <- rbind(bridge[[i]],tmp)
			}
		}
	}
	xyz <- matrix(0,0,3)
	edge <- matrix(0,0,2)
	triangles <- matrix(0,0,3)
  n.triangles <- rep(0,n.psi-1)
	for(i in 1:n.psi){
		n.r <- nrow(xyz)
		if(i > 1){
			pre <- (n.r-length(thetas[[i-1]])+1):n.r
			post <- (n.r+1):(n.r+length(thetas[[i]]))
			edge <- rbind(edge,cbind(post,c(post[-1],post[1])))
			br <- bridge[[i-1]]
			new.edge <- cbind(pre[br[,1]],post[br[,2]])
			edge <- rbind(edge,new.edge)
			tmp.tri <- cbind(new.edge,rbind(new.edge[-1,],new.edge[1,]))
			tmp <- apply(tmp.tri,1,unique)
			triangles <- rbind(triangles,t(tmp))
      n.triangles[i-1] <- length(tmp[1,])
		}
		psi <- psis[i]
		theta <- thetas[[i]]
		xyz <- rbind(xyz,cbind(cos(psi) * cos(theta),cos(psi)*sin(theta),sin(psi)))
		
	}
	return(list(xyz=xyz,edge=edge,triangles=triangles,n.psi=n.psi,thetas=thetas,n.triangles=n.triangles))
}

```

```{r setup}
knit_hooks$set(rgl = hook_rgl)
```
```{r}
n.psi <- 15
sp.mesh <- my.sphere.tri.mesh(n.psi)
```

```{r,rgl=TRUE}
# 打点して
plot3d(sp.mesh$xyz)
# エッジを結ぶ
segments3d(sp.mesh$xyz[c(t(sp.mesh$edge)),])
```

三角形を配置する
```{r,rgl=TRUE}
# 打点して
plot3d(sp.mesh$xyz)
# 三角形メッシュオブジェクトを作り
mesh.tri <- tmesh3d(t(sp.mesh$xyz),t(sp.mesh$triangles),homogeneous=FALSE)
# 三角形を灰色で塗る
shade3d(mesh.tri,col="gray")
```

## 曲面メッシュ情報の四元数化

３次元座標を四元数の虚部の３要素に対応付けることで、座標の操作・ベクトル計算が容易になるので、四元数表現を取り入れる。
onionパッケージのHi,Hj,Hkはそれぞれ四元数虚部３軸の単位元に相当する。

四元数パッケージの特性から、頂点の３座標値を列ベクトルとしてvertices.matという行列に登録し、三角形の３頂点IDを列ベクトルとしてfaces.vに登録する。

```{r}
vertices.mat <- t(sp.mesh$xyz) # ３行の行列にする
# 座標を純虚四元数化
vertices <- vertices.mat[1,]*Hi + vertices.mat[2,]*Hj + vertices.mat[3,]*Hk
edges <- sp.mesh$edge
faces.v <- t(sp.mesh$triangles)
n.v <- length(vertices) # 頂点数
n.f <- length(faces.v[1,]) # 三角形数
```

## 共形変換処理のための諸要素、その関係

### 変形に関する情報

局所における平均曲率の変化量に相当する値rho(実数)を頂点・三角形に与え、平均曲率がその量で変化するように変形する。
実装にあたっては、頂点にrho値を定義しする場合には、
三角形が表す小領域のrho値は、三頂点のrho値の算術平均とする。
三角形自体にrho値を与えてもよい。

以下のRコマンドによる初期化は、各オブジェクトが実数であるのか、四元数であるのか、ベクトルであるのか、行列であるのか、そのサイズはいくつなのかをわかりやすくすることを目的としている。

たとえば、rho.vは頂点のrho値を持つベクトルで、その初期値は実数としての0であり、長さは、頂点数n.vである、と示してある。
同様に三角形のrho値、rho.f (f:face)は、三角形数n.f個の要素の実ベクトルである。

rho.vとメッシュ構成faces.vとからrho.fを計算する関数がrho.fromVtoTriである。
```{r}
rho.v <- rep(0,n.v) # 頂点における
rho.f <- rep(0,n.f) # 三角形における
# 関数 rho.fromVtoTri()は後掲
#rho.f <- rho.fromVtoTri(rho.v,faces.v)
```

### 共形変換に関する情報

曲面の各所で、指定されたrhoに近い値を満足しつつ、曲面全体で、単位球面からの共形変換であるように変形したい。

共形変換では、各点に局所回転を定め、その局所回転が滑らかに移行することである。
この点については、本書では詳述しない。
また、三次元空間座標の回転は、四元数とその共役とで挟んだ積であることから、各点に局所回転を定めることは、各点に四元数を定めることである。
四元数における三次元空間座標回転についても、詳述しない。
この局所回転に相当する四元数を算出するために、面の座標・三角メッシュ・三角形のrhoから作られる頂点数ｘ頂点数の四元数行列Eを作る。

実際には、Eは実、i,j,k要素ごとに頂点数ｘ頂点数の長さの疎ベクトルとして作成し、そこから(頂点数ｘ４)ｘ(頂点数ｘ４）の実数版行列E.reを作成する。

ここで言う、疎ベクトルはMatrixパッケージが提供する、疎ベクトル・疎行列オブジェクトのそれである。
疎ベクトル・疎行列オブジェクトでは、0の多いベクトル・行列の場合に、非0値とその番地を情報として把持することで、メモリを節約し、そのようなデータの持ち方をしながら、線形代数演算を効率に行えるような関数を実装することで、大型ながら疎な行列の演算が可能となっている。

以下の例では、初期疎ベクトルとして、値0を第i=1番地に持つ、長さが頂点数の二乗であるようなベクトルである。

```{r}
E <- list(rep(sparseVector(c(0),i=c(1),length=n.v^2),4))
E.re <- matrix(0,n.v*4,n.v*4)
# 関数 my.make.E.v()は後掲
# E <- my.make.E.v(vertices,faces.v,rho.f)
# 関数 my.qMtorM()は後掲
# E.re <- my.qMtorM(E)
```

Eの固有ベクトルの一つを近似的に求めるにあたり、逆冪乗法を用いる。
そのようにして求めた近似解がlambda.vである。(頂点数ｘ４）ｘ（頂点数ｘ４)の実行列の固有ベクトルであり、長さが(頂点数ｘ４)の実ベクトルである。
行列Eを頂点数ｘ頂点数の四元数行列であると考えれば、得られる固有ベクトルは頂点数の長さの四元数ベクトルである、ともみなせる。

逆冪乗法では、適当な初期ベクトルx0を与え、Ax1=x0を解き、次いでAx2=x1を解くという手順を繰り返す。

繰り返し数はそれほど多くなくても収束がよいことが知られているので、３回繰り返しをデフォルトとして採用している。

```{r}
lambda.v <- rep(0,n.v*4) # 四元数ベクトルを長さが４倍の実ベクトルとして表したもの
# lambda.v <- rep(0*Hi,n.v) # 四元数ベクトル
# lambda.v <- my.inv.pow.2(E.re)[[1]]
```

### 変形後座標の計算

Lは曲面の定義情報(頂点座標と三角形メッシュ情報)とから定まる曲面の広がりを表した頂点数ｘ頂点数の実行列。

変形後の頂点座標(純虚四元数)を最小二乗法で推定するにあたり、実行列の各要素を実部のみの四元数と見立てて、(頂点数ｘ４)ｘ(頂点数ｘ４)の実行列L.reに変換し、L x = omega をxについて解くことで、変形後頂点座標を得る。

omegaは各点に推定した回転相当の四元数のベクトルを用いて算出される、四元数ベクトル。

```{r}
L <- matrix(0,n.v,n.v)
L.re <- matrix(0*Hi,n.v*4,n.v*4)
# 関数 my.make.L(),my.make.quatL()は後掲
# L <- my.make.L(vertices,faces.v)
# L.q <- my.make.quatList(L)
# L.re <- my.qMtorM(L.q)
omega <- rep(0*Hi,n.v)
# 関数 my.make.omega()は後掲
# omega <- my.make.omega(vertices,faces.v,lambda.v)
new.vertices <- rep(0*Hi,n.v)
# new.vertices <- as.quaternion(matrix(solve(L.re,omega),nrow=4))
```

### vertices, faces.v,rho,L,E,lambda,omega,new.verticesの関係

入力は頂点座標vertices、頂点が作る三角形メッシュfaces.v、曲がり方情報rhoの３つ。

出力は、変形後頂点座標のnew.vertices。

中間オブジェクトは、固有ベクトルの元となる行列Eが３入力オブジェクトから作成され、それから一意な解(最小固有値に対応する固有ベクトル)の近似解lamdaが得られる。
もう一つの中間オブジェクトはラプラシアン行列で、vertices,faces.vの２つから作成される。
さらなる下流の中間オブジェクトはomegaで、これは、vertices,faces.v,lambdaから算出するが、lambda自体がvertices,faces.vから作成されているから、omegaはlambdaを介してvertices,faces.v,rhoから作成されるものである。

最終出力は、vertices,faces.vのみから作成されるLと、rhoも加わって作成されるomegaとの２つから算出される。

```{r}
lab <- c("v","f","rho","L","E","lmd","omg","n.v")
edge.list <- matrix(c("v","L","f","L","v","E","f","E","rho","E","E","lmb","v","omg","f","omg","lmb","omg","L","n.v","omg","n.v"),byrow=TRUE,ncol=2)
g.rel <- graph.edgelist(edge.list)
plot(g.rel)                 
```

# 関数

## 基礎となる諸関数

以下、関数ソースを記しながら、若干のコメントを加えておく。
```{r}
# 頂点のrho値から面のrho値を算出する
rho.fromVtoTri <- function(rho.v,faces.v){
  tmp.rho <- matrix(rho.v[faces.v],nrow=3)
  apply(tmp.rho,2,mean)
}
# Utility 関数
# 特に疎行列において、面単位で不定回数の値加算をするためのユーティリティ関数
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))
}
# E行列作成関数
# ただし、四元数行列なので、虚実４成分それぞれをリストの要素とし
# 個々のリストも行列ではなく、ベクトル化した形で返す
# 疎行列のリスト
# 三角形の面積を考慮して、面のrho値に応じて、rhoが作る曲面上の
# ベクトル場を作り
# その定常解の一つとして固有ベクトルを取る
my.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 <- Mod(Im(tmp))/2
	#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.i + sparseVector(tmp.out[[2]][,2],tmp.out[[1]],length(vertices)^2)
			E.j <- E.j + sparseVector(tmp.out[[2]][,3],tmp.out[[1]],length(vertices)^2)
			E.k <- E.k + 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)
}
# Ax = lambda b の固有ベクトルを逆冪乗法で近似する
# 初期ベクトル bは(1,1,...)がよいことが知られているので
# それをデフォルト値としてある
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))))
	}
}
# ラプラシアン行列の算出
# すべての点について、その点を含む三角形について、三角形の内角で定まる値を加算する
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.quatList <- function(L){
  L.q <- list()
  L.q[[1]] <- L
  for(i in 2:4){
    L.q[[i]] <- L*0
  }
  L.q
}
# omegaは回転を定めるlambdaを介して算出する
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)
}
# 局所平均曲率を計算する関数
# 入力は頂点の四元数座標ベクトルと３行で表された三角形頂点ID行列
# 返り値はnorm.v,norm.face=ret.face,dirの３つで
# norm.vは各頂点の法線ベクトルに相当する四元数の行列で、その虚部の絶対値が平均曲率の大きさ
# norm.faceは各三角形のそれで、虚部の絶対値が平均曲率の大きさ
# dirは各三角形の平均曲率が曲面の外向きか内向きか(今、扱っている曲面はすべて球面様の閉曲面なので、外向き・内向きはその意味での内外)を示すc(-1,0,1)
my.curvature.cot <- function(vertices,faces.v){
  n.v <- length(vertices)
  ret <- rep(0*Hk,n.v)
  # 三角形の面積
  edge1 <- vertices[faces.v[2,]]-vertices[faces.v[1,]]
  edge2 <- vertices[faces.v[3,]]-vertices[faces.v[1,]]
  tmp <- edge1 * edge2
  tmp2 <- Mod(Im(tmp))
  #tmp2 <- i(tmp)+j(tmp)+k(tmp)
  #inv.f <- which(tmp2<0)
  #faces.v[2:3,inv.f] <- rbind(faces.v[3,inv.f],faces.v[2,inv.f])
 	A.f <- (tmp2)/2
  # 頂点周囲の三角形面積の和
  val <- rep(A.f,each=3)
  addr <- c(faces.v)
  tmp.out <- my.vector.access(val,addr)
  A.v <- tmp.out[[2]][,1] # 頂点周囲面積和
  
  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
    
    # 対向辺
    u3 <- v3-v2
    #u3.len <- Mod(Im(u3))
		# edge 四元数(純虚四元数)の積は実部がドット積、虚部がクロス積ベクトル
		u12 <- u1 * u2
		cotAlpha <- ((-Re(u12))/Mod(Im(u12)))
    
    val <- c(cotAlpha * u3, cotAlpha * (-u3))
    addr <- c(k2,k3)
    tmp.out <- my.vector.access(t(as.matrix(val)),addr)
    tmp.val <- as.quaternion(t(tmp.out[[2]]))
    ret[tmp.out[[1]]] <- ret[tmp.out[[1]]] + tmp.val
  }
  ret.vec <- -ret/(4*A.v)
  ret.face.re <- rho.fromVtoTri(Re(ret.vec),faces.v)
  ret.face.i <- rho.fromVtoTri(i(ret.vec),faces.v)
  ret.face.j <- rho.fromVtoTri(j(ret.vec),faces.v)
  ret.face.k <- rho.fromVtoTri(k(ret.vec),faces.v)
  ret.face <- ret.face.re + Hi * ret.face.i + Hj * ret.face.j + Hk * ret.face.k
  
  
  dir <- sign(Re(tmp * ret.face))
  return(list(norm.v=ret.vec,norm.face=ret.face,dir=dir))
}

# ３列の実行列(３次元座標の行列)xyzと、３行の頂点ID行列faces.vとを用いて、三角形メッシュを塗り描く関数
# ただし、この関数は、元がsp.triは、元の三角形メッシュがmy.sphere.tri.mesh()によって作られていることを前提としている。そのような場合に、sp.triという引数によって、隣接緯線間ごとに色を塗り分けることで、元の正球からの変形の様子を色表現することができる
# また、各三角形の色はrho.fなる実ベクトルによってコントラストをつけることができる
# col1は、元の正球からの変化を縞模様化するための長さ２の整数ベクトルである。その縞模様を使いたくないとき(メッシュがmy.sphere.tri.mesh()由来ではないときを含む)は、この長さ２のベクトルの２つの値を同じにすればよい)
plot.sp.conformal <- function(xyz,faces.v,sp.tri,rho.f,col1=c(4,5)){
  plot3d(xyz,xlab="x",ylab="y",zlab="z")
  mesh.tri <- tmesh3d(t(xyz),faces.v,homogeneous=FALSE)
  
  # 縞のための値ベクトル
  col. <- rep(col1,length(sp.tri))[1:length(sp.tri)]
  col <- rep(col.,sp.tri*3)
  # rhoを反映した値ベクトル
  rho.f <- rep(rho.f,each=3)
  #rho.f2 <- sign(rho.f)
  rho.f <- (rho.f-min(rho.f))/(max(rho.f)-min(rho.f))
  
  col2 <- rgb(1-rho.f,1,col/6)
  #col3 <- gray((rho.f2+1)*0.5)
  shade3d(mesh.tri,col=col2)  
}
```

## 曲面の変形処理の関数

曲面が三角化されており、
初期頂点座標を四元数ベクトルverticesとして与え、三角形の頂点ID行列が行数３の行列faces.vとして与え、平均曲率の変化量をrho.vとして与える。ただし、face=FALSEの場合には、rho.vは頂点ごとにその変化量を与え、face=TRUEの場合には、三角形ごとにその変化量を与えるものとする。

```{r}
my.conformal.rho <- function(vertices,faces.v,rho.v,face=FALSE){
  xyz.ori <- t(as.matrix(vertices)[2:4,])
  n.v <- length(vertices) # 頂点数
  n.f <- length(faces.v[1,]) # 三角形数
  if(!face){
    rho.f <- rho.fromVtoTri(rho.v,faces.v)
  }else{
    rho.f <- rho.v
  }
  
  
  # 三角形の面積
  edge1 <- vertices[faces.v[2,]]-vertices[faces.v[1,]]
	edge2 <- vertices[faces.v[3,]]-vertices[faces.v[1,]]
	tmp <- edge1 * edge2
  A <- Mod(Im(tmp))/2
  
  # rho の面積重みつき総和は0でないと「閉じ」ない
  s <- sum(A*rho.f)
  rho.f <- rho.f -s*A/sum(A)
  
  
  E <- my.make.E.v(vertices,faces.v,rho.f)
  E.re <- my.qMtorM(E)
  
  lambda.v <- my.inv.pow.2(E.re)[[1]]
  
  L <- my.make.L(vertices,faces.v)
  L.q <- my.make.quatList(L)
  L.re <- my.qMtorM(L.q)
  
  omega <- my.make.omega(vertices,faces.v,lambda.v)
  
  new.vertices <- as.quaternion(matrix(solve(L.re,omega),nrow=4))
  xyz.new <- t(as.matrix(new.vertices)[2:4,])
  mean.new <- apply(xyz.new,2,mean)
  xyz.new. <- t(t(xyz.new)-mean.new)
  max.new. <- max(abs(xyz.new.))
  xyz.new.st <- xyz.new./max.new.
  
  new.q <- as.quaternion(t(cbind(rep(0,n.v),xyz.new.st)))
  #ret <- xyz.new.st[,1]*Hi + xyz.new.st[,2]*Hj + xyz.new.st[,3]*Hk
  
  ret <- list(xyz.new=xyz.new.st,xyz.ori=sp.mesh$xyz,xyz.new.q=new.q,xyz.ori.q=vertices,faces.v=faces.v,E=E,L=L,lambda.v=lambda.v,omega=omega,n.psi=n.psi,rho.fx=rho.fx,rho.v=rho.v,rho.f=rho.f,sp.mesh=sp.mesh)
  ret
}
```

また、緯度分割数n.psiを与えることで単位球面メッシュを作成し、その単位球面上の頂点にrho.fx()なる関数で産出されるrhoの変化量が与えられたとしたときの、my.conformal.rho()関数の結果を返す関数をmy.sphereCoformal()とする。

単位球面からの変形曲面の作成のための関数である。

```{r}
my.sphereConformal <- function(n.psi,rho.fx){
  sp.mesh <- my.sphere.tri.mesh(n.psi)
  vertices.mat <- t(sp.mesh$xyz) # ３行行列化
  vertices <- vertices.mat[1,]*Hi + vertices.mat[2,]*Hj + vertices.mat[3,]*Hk
  edges <- sp.mesh$edge
  faces.v <- t(sp.mesh$triangles)
  rho.v <- rho.fx(sp.mesh$xyz) # 頂点における
  out <- my.conformal.rho(vertices,faces.v,rho.v)
  ret <- list(xyz.new=out$xyz.new,xyz.ori=sp.mesh$xyz,xyz.new.q=out$xyz.new.q,xyz.ori.q=vertices,faces.v=faces.v,E=out$E,L=out$L,lambda.v=out$lambda.v,omega=out$omega,n.psi=n.psi,rho.fx=rho.fx,rho.v=rho.v,rho.f=out$rho.f,sp.mesh=sp.mesh)
  ret
}
```

# 実行例

## 例１

単位球状の三角形メッシュを作り、
平均曲率の変化量が球面座標の関数になっているものとして、適当な関数を決めて、共形変換した形を描いてみよう。

Z軸方向にサインカーブ様の平均曲率変化を与えてみる。
```{r}
n.psi <- 30 # 球面メッシュの緯度刻み数
rho.fx <- function(X){
  ret <- sin(X[,3]*pi*2 )*3
  return(ret)
}
out1 <- my.sphereConformal(n.psi,rho.fx)
```


緯度によって縞模様をつけた上で、rho値でコントラストをつけてrho値の勾配をつけてプロットする。

変形前。
```{r,rgl=TRUE}
plot.sp.conformal(out1$xyz.ori,out1$faces.v,out1$sp.mesh$n.triangles,out1$rho.f)
```

変形後。

青緑の濃い部分は、外に突き出していて、クリーム色の部分は、凹んでいる。


```{r,rgl=TRUE}
plot.sp.conformal(out1$xyz.new,out1$faces.v,out1$sp.mesh$n.triangles,out1$rho.f)
```

変形後の座標から、各三角形の平均曲率を算出し、それに応じて色を付け直してみる。

色の濃淡に若干の変化はあるも、突出部分で緑勝ちであり、凹みでクリーム色であることに違いはない。

```{r,rgl=TRUE}
# 平均曲率の計算
rho.cot1 <- my.curvature.cot(out1$xyz.new.q,out1$faces.v)
# 三角形の平均曲率を符号付きで出しなおす
rho.face1 <- rho.cot1[[3]] * Mod(Im(rho.cot1[[2]]))
plot.sp.conformal(out1$xyz.new,out1$faces.v,out1$sp.mesh$n.triangles,rho.face1)
```

変形のために指定したrho値と、変形後の座標から算出した平均曲率との関係を散布図で見てみる。

rhoと平均曲率とは大まかにそろっていて、特に、凸は凸、凹は凹の関係にある。

```{r}
plot(out1$rho.f,rho.face1,xlab="rho",ylab="mean curvature")
abline(h=0,col=2)
abline(v=0,col=2)
abline(0,1,col=3)
```

## 例２

盛り上がりとくぼみとを指定してみる。

```{r}
n.psi <- 30
rho.fx <- function(X){
  ret <- (2*(X[,1]-0.3)^2+1)*1 - (3*(X[,2]+0.1)^4+2)*1 
  return(ret)
}

out2 <- my.sphereConformal(n.psi,rho.fx)
```

変形前。
```{r,rgl=TRUE}
plot.sp.conformal(out2$xyz.ori,out2$faces.v,out2$sp.mesh$n.triangles,out2$rho.f)
```

変形後。

```{r,rgl=TRUE}
plot.sp.conformal(out2$xyz.new,out2$faces.v,out2$sp.mesh$n.triangles,out2$rho.f)
```

変形後の座標から、各三角形の平均曲率を算出し、それに応じて色を付け直してみる。

```{r,rgl=TRUE}
# 平均曲率の計算
rho.cot2 <- my.curvature.cot(out2$xyz.new.q,out2$faces.v)
# 三角形の平均曲率を符号付きで出しなおす
rho.face2 <- rho.cot2[[3]] * Mod(Im(rho.cot2[[2]]))
plot.sp.conformal(out2$xyz.new,out2$faces.v,out2$sp.mesh$n.triangles,rho.face2)
```

変形条件rhoと、変形後の平均曲率との関係。

初めの例よりは、指定値と変形後曲率との関係の一致が悪いが、努めて指定値を実現しようとした結果である様子とでもいえばよいだろう。


```{r}
plot(out2$rho.f,rho.face2,xlab="rho",ylab="mean curvature")
abline(h=0,col=2)
abline(v=0,col=2)
abline(0,1,col=3)
```

## 球に戻す変形

正球の各所に変形後の平均曲率のようなものを指定して変形することができた。
この指定値は、平均曲率の変化量であったから、球ではない閉曲面に、その平均曲率分だけ、平均曲率が減るようにrho値を与えれば、球に戻るだろう。

このような処理により、凹凸を減らした形を作ることができる。

このような球に戻る処理を段階的に行うために、関数を用意する。

ステップごとに、三角形の平均曲率を算出し、その大きさに応じて変化量rhoを与えることを繰り返す。

頂点座標の四元数ベクトルと、行数３の頂点ID行列とともに、平均曲率のk倍だけ減じるように指定する引数kと、繰り返し回数n.iterを与える。
```{r}
my.deform.k.serial <- function(vertices,faces.v,k,n.iter=10){
  v.list <- rho.cot.list <- list()
  v.list[[1]] <- vertices
  rho.cot.list[[1]] <- my.curvature.cot(v.list[[1]],faces.v)
  for(i in 1:n.iter){
    
    tmp.rho.f <- rho.cot.list[[i]][[3]] * Mod(Im(rho.cot.list[[i]][[2]]))
    tmp.out <- my.conformal.rho(v.list[[i]],faces.v,k*tmp.rho.f,face=TRUE)
    v.list[[i+1]] <- tmp.out$xyz.new.q
    rho.cot.list[[i+1]] <- my.curvature.cot(v.list[[i+1]],faces.v)
  }
  return(list(v.list=v.list,rho.cot.list=rho.cot.list,k=k))
}
```

```{r}
k <- -1.5
n.iter <- 10
out.series1 <- my.deform.k.serial(out1$xyz.new.q,out1$faces.v,k,n.iter=n.iter)
```
```{r,rgl=TRUE}
i <- 1
xyz <- t(as.matrix(out.series1$v.list[[i]])[2:4,])
rho.f <- out.series1$rho.cot.list[[i]][[3]]*Mod(Im(out.series1$rho.cot.list[[i]][[2]]))
plot.sp.conformal(xyz,out1$faces.v,out1$sp.mesh$n.triangles,rho.f)
i <- i+1
```
```{r,rgl=TRUE}
xyz <- t(as.matrix(out.series1$v.list[[i]])[2:4,])
rho.f <- out.series1$rho.cot.list[[i]][[3]]*Mod(Im(out.series1$rho.cot.list[[i]][[2]]))
plot.sp.conformal(xyz,out1$faces.v,out1$sp.mesh$n.triangles,rho.f)
i <- i+1
```
```{r,rgl=TRUE}
xyz <- t(as.matrix(out.series1$v.list[[i]])[2:4,])
rho.f <- out.series1$rho.cot.list[[i]][[3]]*Mod(Im(out.series1$rho.cot.list[[i]][[2]]))
plot.sp.conformal(xyz,out1$faces.v,out1$sp.mesh$n.triangles,rho.f)
i <- i+1
```
```{r,rgl=TRUE}
xyz <- t(as.matrix(out.series1$v.list[[i]])[2:4,])
rho.f <- out.series1$rho.cot.list[[i]][[3]]*Mod(Im(out.series1$rho.cot.list[[i]][[2]]))
plot.sp.conformal(xyz,out1$faces.v,out1$sp.mesh$n.triangles,rho.f)
i <- i+1
```
```{r,rgl=TRUE}
xyz <- t(as.matrix(out.series1$v.list[[i]])[2:4,])
rho.f <- out.series1$rho.cot.list[[i]][[3]]*Mod(Im(out.series1$rho.cot.list[[i]][[2]]))
plot.sp.conformal(xyz,out1$faces.v,out1$sp.mesh$n.triangles,rho.f)
i <- i+1
```
```{r,rgl=TRUE}
xyz <- t(as.matrix(out.series1$v.list[[i]])[2:4,])
rho.f <- out.series1$rho.cot.list[[i]][[3]]*Mod(Im(out.series1$rho.cot.list[[i]][[2]]))
plot.sp.conformal(xyz,out1$faces.v,out1$sp.mesh$n.triangles,rho.f)
i <- i+1
```
```{r,rgl=TRUE}
xyz <- t(as.matrix(out.series1$v.list[[i]])[2:4,])
rho.f <- out.series1$rho.cot.list[[i]][[3]]*Mod(Im(out.series1$rho.cot.list[[i]][[2]]))
plot.sp.conformal(xyz,out1$faces.v,out1$sp.mesh$n.triangles,rho.f)
i <- i+1
```
```{r,rgl=TRUE}
xyz <- t(as.matrix(out.series1$v.list[[i]])[2:4,])
rho.f <- out.series1$rho.cot.list[[i]][[3]]*Mod(Im(out.series1$rho.cot.list[[i]][[2]]))
plot.sp.conformal(xyz,out1$faces.v,out1$sp.mesh$n.triangles,rho.f)
i <- i+1
```
```{r,rgl=TRUE}
xyz <- t(as.matrix(out.series1$v.list[[i]])[2:4,])
rho.f <- out.series1$rho.cot.list[[i]][[3]]*Mod(Im(out.series1$rho.cot.list[[i]][[2]]))
plot.sp.conformal(xyz,out1$faces.v,out1$sp.mesh$n.triangles,rho.f)
i <- i+1
```
```{r,rgl=TRUE}
xyz <- t(as.matrix(out.series1$v.list[[i]])[2:4,])
rho.f <- out.series1$rho.cot.list[[i]][[3]]*Mod(Im(out.series1$rho.cot.list[[i]][[2]]))
plot.sp.conformal(xyz,out1$faces.v,out1$sp.mesh$n.triangles,rho.f)
i <- i+1
```


kの値を色々変えて変形してみる。ksでkの値を指定する。
```{r}
my.deform.k.multi <- function(vertices,faces.v,ks){
  v.list <- rho.cot.list <- list()
  v.list[[1]] <- vertices
  rho.cot.list[[1]] <- my.curvature.cot(v.list[[1]],faces.v)
  for(i in 1:length(ks)){
    
    tmp.rho.f <- rho.cot.list[[1]][[3]] * Mod(Im(rho.cot.list[[1]][[2]]))
    tmp.out <- my.conformal.rho(v.list[[1]],faces.v,ks[i]*tmp.rho.f,face=TRUE)
    v.list[[i+1]] <- tmp.out$xyz.new.q
    rho.cot.list[[i+1]] <- my.curvature.cot(v.list[[i+1]],faces.v)
  }
  return(list(v.list=v.list,rho.cot.list=rho.cot.list,k=k))
}
```

```{r}
ks <- seq(from=-3,to=3,length=7)
out.series <- my.deform.k.multi(out1$xyz.new.q,out1$faces.v,ks)
```

```{r,rgl=TRUE}
i <- 1
xyz <- t(as.matrix(out.series1$v.list[[i]])[2:4,])
rho.f <- out.series1$rho.cot.list[[i]][[3]]*Mod(Im(out.series1$rho.cot.list[[i]][[2]]))
plot.sp.conformal(xyz,out1$faces.v,out1$sp.mesh$n.triangles,rho.f)
i <- i+1
```
```{r,rgl=TRUE}
xyz <- t(as.matrix(out.series1$v.list[[i]])[2:4,])
rho.f <- out.series1$rho.cot.list[[i]][[3]]*Mod(Im(out.series1$rho.cot.list[[i]][[2]]))
plot.sp.conformal(xyz,out1$faces.v,out1$sp.mesh$n.triangles,rho.f)
i <- i+1
```
```{r,rgl=TRUE}
xyz <- t(as.matrix(out.series1$v.list[[i]])[2:4,])
rho.f <- out.series1$rho.cot.list[[i]][[3]]*Mod(Im(out.series1$rho.cot.list[[i]][[2]]))
plot.sp.conformal(xyz,out1$faces.v,out1$sp.mesh$n.triangles,rho.f)
i <- i+1
```
```{r,rgl=TRUE}
xyz <- t(as.matrix(out.series1$v.list[[i]])[2:4,])
rho.f <- out.series1$rho.cot.list[[i]][[3]]*Mod(Im(out.series1$rho.cot.list[[i]][[2]]))
plot.sp.conformal(xyz,out1$faces.v,out1$sp.mesh$n.triangles,rho.f)
i <- i+1
```
```{r,rgl=TRUE}
xyz <- t(as.matrix(out.series1$v.list[[i]])[2:4,])
rho.f <- out.series1$rho.cot.list[[i]][[3]]*Mod(Im(out.series1$rho.cot.list[[i]][[2]]))
plot.sp.conformal(xyz,out1$faces.v,out1$sp.mesh$n.triangles,rho.f)
i <- i+1
```
```{r,rgl=TRUE}
xyz <- t(as.matrix(out.series1$v.list[[i]])[2:4,])
rho.f <- out.series1$rho.cot.list[[i]][[3]]*Mod(Im(out.series1$rho.cot.list[[i]][[2]]))
plot.sp.conformal(xyz,out1$faces.v,out1$sp.mesh$n.triangles,rho.f)
i <- i+1
```
```{r,rgl=TRUE}
xyz <- t(as.matrix(out.series1$v.list[[i]])[2:4,])
rho.f <- out.series1$rho.cot.list[[i]][[3]]*Mod(Im(out.series1$rho.cot.list[[i]][[2]]))
plot.sp.conformal(xyz,out1$faces.v,out1$sp.mesh$n.triangles,rho.f)
i <- i+1
```

#

# 背景知識

## 四元数

### 四元数の基礎
四元数は１つの実部と３種類の虚部の計４成分からなる数である。
４成分をそれぞれ(1,i,j,k)と表し、Rのonionパッケージでは(1,Hi,Hj,Hk)がそれに相当する。
４成分同士は積で入れ替わる関係にあり、以下のような演算表で表される関係にある。

$$
\begin{pmatrix}1 & i & j & k // i & -1 & k & -j // j & -k & -1 & i // k & j -i & -1 \end{pmatrix}$$

```{r}
library(onion)
four <- c(1+0*Hi,0+1*Hi,0+1*Hj,0+1*Hk)
1 * four
Hi * four
Hj * four
Hk * four
```
### 四元数による３次元座標表示

四元数の３つの虚部(i,j,k)を３次元直交座標の(x,y,z)に見立てることで、四元数による３次元座標表現とする。

例えば(x,y,z)=(1,2,3)は、(i,j,k)=(1,2,3)であり、0+1i+2j+3kという四元数とみなす。

## 四元数表現された３次元ベクトルの積の幾何学的意味

二つの３次元ベクトルを２つの純虚四元数で表すとき、その積の実部はベクトルの内積に(-1)を掛けたものになり、虚部はベクトルのクロス積に一致する。

このクロス積は、２ベクトルが張る平面に垂直なベクトルである。
また、このクロス積の長さは２ベクトルが作る三角形の面積の２倍でもある。

２ベクトルのなす角のcotangentは内積を面積で割ったものである。

```{r}
imq <- c(Hi,Hj,Hk)
u <- runif(3)
v <- runif(3)
u.q <- sum(u*imq) # uの四元数表現
v.q <- sum(v*imq) # vの四元数表現
u.v.q.prod <- u.q * v.q # 四元数の積
Re(u.v.q.prod) # 四元数の積の実部
Im(u.v.q.prod) # 四元数の積の虚部
sum(u*v) # ベクトルの内積
library(geometry) # extprod3d()関数を持つパッケージ
extprod3d(u,v) # ３次元ベクトルのクロス積
# 三角形の面積
area.triangle <- Mod(Im(u.v.q.prod))/2
# cotangent
cot <- (-Re(u.v.q.prod))/Mod(Im(u.v.q.prod))
```

四元数を用いた三角形の面積の計算は関数my.make.E.v()やmy.conformal.rho()の中で使われている。

cotangentの計算は同じくmy.make.M()やmy.make.omega(),my.conformal.rho()の中で使われている。


### 四元数による３次元回転

３次元座標を純虚四元数pで表しているとき、ある四元数qとその共役四元数q.とで挟んだ積 q p q. は、ある回転を表す。その回転とはqが純虚単位四元数uを用いて$\cos{\frac{\theta}{2}} + \sin{\frac{\theta}{2}} u$と表せるときに、uが表す３次元ベクトルを軸とした$\theta$の回転に相当する。

```{r}
theta <- pi/6
u <- sum(runif(1) * imq)
u <- u/Mod(u) # 単位純虚四元数
q <- cos(theta/2) + sin(theta/2) * u
q. <- cos(theta/2) - sin(theta/2) * u
p <- sum(runif(1) * imq)
p.rot <- q * p * q.
# pとp.rotの長さはともに純虚四元数
Re(p)
Re(p.rot)
# 純虚四元数の長さは等しい
Mod(Im(p))
Mod(Im(p.rot))
# pとqの虚部、p.rotとqの虚部のなす角は等しい(内積が等しければ、長さが等しいので角も等しい)
Re(p*u)
Re(p.rot*u)
# 回転角を確認するにはpとp.rotとのqに垂直な成分を取り出し、そのなす角を調べる必要がある
p.n <- p - (-Re(p*q)*q)
p.rot.n <- p.rot - (-(Re(p.rot*q)*q))
acos(-Re(p.n*p.rot.n)/(Mod(p.n)*Mod(p.rot.n)))
theta
```

四元数を用いた回転は関数my.make.omega()内で使われている。


### 四元数の行列表現

四元数の行列表現とは、以下のような関数で作られる行列のことである。

```{r}
my.make.Matq <- function(q){
  v <- c(Re(q),i(q),j(q),k(q))
	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)
  matrix(v[a]*b,4,4)
}
q.m <- my.make.Matq(q)
p.m <- my.make.Matq(p)
q.m
```  

このようにすることで、四元数の和の行列表現は、四元数の行列表現の和に一致し、四元数の積の行列表現は、四元数の行列表現の積に一致する。

したがって、四元数の行列表現は、和と積について、四元数でやっても行列でやっても同じである、と言う意味で、「同じもの」と言える。

これにより、四元数の推定等において、実数行列演算を用いることが可能になる。

```{r}
my.make.Matq(q+p) - (q.m + p.m)
my.make.Matq(q*p) - q.m %*% p.m
my.make.Matq(p*q) - p.m %*% q.m
```

## その他の大事な点

E行列の作成と、その固有ベクトル算出が局所回転推定になることや、L行列の作成、omegaベクトルの算出などは、http://www.cs.columbia.edu/~keenan/index.html からのリンク先である
http://www.cs.columbia.edu/~keenan/Projects/DGPDEC/　
やhttp://www.cs.columbia.edu/~keenan/Projects/SpinTransformations/　を参考にすること。