---
title: "正八面体からの球面オイラー三角化とその木表現"
author: "ryamada"
date: "Sunday, June 28, 2015"
output: html_document
---
球面オイラー三角化の木表現について検討をしたいのだが、バリエーションのある球面オイラー三角化を作り出したり、それを視覚表現したりすることが面倒である。
この点を解決するために、正八面体から単純なアルゴリズムで球面オイラー三角化を作成するとともに、オイラー三角化グラフと１対１対応する木構造の取り出しをRで実行できるようにする。
球面を三角形でタイリングし、その三角形全体が、隣接三角形の色が同じにならないようにして、２色で塗り分けられるとき、球面オイラー三角化という。
オイラー三角化では、面数は偶数(２色のそれぞれの色の面の数が等しい)、すべての頂点は偶数次数を持つ、という特徴がある。
球面オイラー三角化グラフの１辺をルート辺と定め、その向きを決める。
そうすると、すべてのタイル三角形がサイクルになるようにすべてのエッジに向きを定めることができ、その定め方は、ルート辺とその向きの取り方によって一意に決まる。
具体的な手続きとしては、
(1) ルート辺とその向きを定める。
(2) ルート辺の始点からのグラフ距離を測る。
(3) オイラーの三角形の２色塗りをしたときに、ルート辺を含む時計回りのサイクルが作る三角形と同じ色になる三角形について、３頂点のグラフ距離値が大きい方の２頂点を結ぶエッジを取り出す。ルート辺も加える。
これが球面オイラー三角化を表す木構造である。
木構造の辺のうち、ルート辺以外は、すべてある色の三角形に１対１対応する。
"generating spherical eulerian triangulation" (http://www.sciencedirect.com/science/article/pii/S0012365X10000026) というタイトルの短い論文によれば
(1) シンプルな無向グラフの場合に限ると
(2) すべての球面に埋め込まれたオイラー三角化グラフは、正八面体から始めて、２つの処理を繰り返すことで作成できると言う
(3) ２つの処理をｐ処理とｑ処理と呼ぶことにする
(3-1) p処理
> 今、あるノードの次数が３以上であるとする(すべてのノードは偶数次数であるから実際には４以上)
> このとき、平面グラフとしてみたときに隣り合う３エッジについて次のような処理をする
> ３エッジの途中にノードを新たに作る。そのノードと残りの２エッジが連結しているノードとの間にエッジを引く
> 結果として、直線状に並んでいた３ノードが十字で仕切られたひし形に変わる
(3-2) q処理
> 今、あるノードの次数が５以上であるとする
> このとき、平面グラフとしてみたときに隣り合う５得次について次のような処理をする
> 時計回りに、第３、第４番目のエッジを除去する
> その除去によって空いたスペースに第２、第５エッジが隣接するノードを結ぶエッジを引く
> その新たに引いたエッジに関して着目ノードと対側に新たにノードをおき、この新ノードに、第２、３，４，５のエッジの端点との間にエッジを引く
> 結果として、直線状に並んでいた３ノードが縦に割られた形での菱形に変わる
まず、正八面体を適当に配置する。
エッジリストel。
3次元頂点座標x(ただし、座標はうまく機能していない)。
各頂点に連結するエッジのリストeofv、ただしこれは平面グラフであることを意識して、エッジの接続順が反時計回りとなるようにしたエッジの連結先頂点IDベクトルである。
三角形を表す3列の行列。３列頂点IDを納める。ただし、2種類を作成する。triは頂点IDを値でソートしたもの。tri.ordはサイクルの向きを考慮し頂点ID順にて作成したもの。
三角形の色を表すベクトルcols。
```{r}
library(igraph)
library(rgl)
# エッジリスト
el <- matrix(c(1,2,1,3,1,4,1,5,2,3,3,4,4,5,5,2,2,6,3,6,4,6,5,6),byrow=TRUE,ncol=2)
# 頂点座標
x <- matrix(c(0,0,1,1,0,0,0,1,0,-1,0,0,0,-1,0,0,0,-1),byrow=TRUE,ncol=3)
x <- x/sqrt(apply(x^2,1,sum))
# 隣接頂点、順序考慮のリスト
eofv <- list()
eofv[[1]] <- c(2,3,4,5)
eofv[[2]] <- c(1,5,6,3)
eofv[[3]] <- c(1,2,6,4)
eofv[[4]] <- c(1,3,6,5)
eofv[[5]] <- c(1,4,6,2)
eofv[[6]] <- c(5,4,3,2)
tris <- matrix(c(1,2,3,1,3,4,1,4,5,1,5,2,6,2,3,6,3,4,6,4,5,6,5,2),byrow=TRUE,ncol=3)
tris <- t(apply(tris,1,sort))
tris.ord <- matrix(c(1,2,3,1,4,3,1,4,5,1,2,5,6,2,3,6,4,3,6,4,5,6,2,5),byrow=TRUE,ncol=3)
cols <- c(FALSE,TRUE,FALSE,TRUE,TRUE,FALSE,TRUE,FALSE)
```
処理P。
隣接頂点リスト eofv、処理対象頂点ID v、三角形頂点ID行列 tris(ID値ソート)、三角形色ベクトル cols、三角形頂点ID行列(サイクル向き)を引数とし、
更新したeofv,x,tris,cols,tris.ordを返す。
```{r}
rule.p <- function(eofv,v,x,tris,cols,tris.ord){
# 処理対象頂点vのどの部分にひし形を挿入するかをランダムに決める
st.v <- sample(1:length(eofv[[v]]),1)
# 扱う隣接頂点のリスト
nb <- c(eofv[[v]],eofv[[v]])[st.v:(st.v+2)]
# 新規頂点のID
newv <- 1:2 + length(eofv)
# 隣接頂点ベクトルの更新
eofv[[nb[1]]] <- my.replace.nodes(eofv[[nb[1]]],v,c(newv[2],newv[1],v))
eofv[[nb[2]]] <- my.replace.nodes(eofv[[nb[2]]],v,c(newv[2]))
eofv[[nb[3]]] <- my.replace.nodes(eofv[[nb[3]]],v,c(v,newv[1],newv[2]))
eofv[[v]] <- my.replace.nodes(eofv[[v]],nb[2],c(newv[1]))
eofv[[newv[1]]] <- c(v,nb[1],newv[2],nb[3])
eofv[[newv[2]]] <- c(newv[1],nb)
# 三角形の色、辺の向きを定める
# 定めるにあたり、処理対称部分の色パターン、辺の向きパターンを確認
tmp.tri <- sort(c(v,nb[1],nb[2]))
tmp.tri.id <- which(tris[,1]==tmp.tri[1] & tris[,2] ==tmp.tri[2] & tris[,3]==tmp.tri[3])
tmp.col <- cols[tmp.tri.id]
add.tris <- matrix(c(v,nb[1],newv[1],v,nb[3],newv[1],newv[1],newv[2],nb[1],newv[1],newv[2],nb[3],newv[2],nb[1],nb[2],newv[2],nb[3],nb[2]),byrow=TRUE,ncol=3)
add.tris.ord <- add.tris
tmp.order.tri <- tris.ord[tmp.tri.id,]
tmp.order.tri2 <- c(tmp.order.tri,tmp.order.tri[1])
loc1 <- which(tmp.order.tri==v)
if(tmp.order.tri2[loc1+1]==nb[1]){
new.tris.ord <- rbind(tris.ord,add.tris.ord)
}else{
new.tris.ord <- rbind(tris.ord,cbind(add.tris.ord[,3],add.tris.ord[,2],add.tris.ord[,1]))
}
add.tris <- t(apply(add.tris,1,sort))
add.cols <- c(TRUE,FALSE,FALSE,TRUE,TRUE,FALSE)
if(!tmp.col){
add.cols <- !add.cols
}
newtris <- rbind(tris,add.tris)
newcols <- c(cols,add.cols)
rmtris1 <- sort(c(v,nb[1],nb[2]))
rmtris2 <- sort(c(v,nb[2],nb[3]))
rmtris1id <- which(newtris[,1]==rmtris1[1] & newtris[,2]==rmtris1[2] & newtris[,3]==rmtris1[3])
rmtris2id <- which(newtris[,1]==rmtris2[1] & newtris[,2]==rmtris2[2] & newtris[,3]==rmtris2[3])
newtris <- newtris[-c(rmtris1id,rmtris2id),]
newcols <- newcols[-c(rmtris1id,rmtris2id)]
new.tris.ord <- new.tris.ord[-c(rmtris1id,rmtris2id),]
newx <- rbind(x,matrix(0,2,3))
id1 <- which(eofv[[v]]==nb[1])
id2 <- which(eofv[[v]]==nb[2])
id3 <- which(eofv[[v]]==nb[3])
rest <- eofv[[v]][-(c(id1,id2,id3))]
newx[newv[1],] <- (2*newx[v,] + newx[nb[2],])/3
newx[newv[2],] <- (newx[v,] + 2 * newx[nb[2],])/3
newx <- newx/sqrt(apply(newx^2,1,sum))
return(list(eofv=eofv,x=newx,tris=newtris,cols=newcols,tris.ord=new.tris.ord))
}
```
同様に処理qを実施
```{r}
rule.q <- function(eofv,v,x,tris,cols,tris.ord){
loop <- TRUE
# 周囲ノードでnb[2],nb[5]として取り出された２頂点間にすでに辺があるときは
# 同一の処理ができないのでそうならない場合を取り出す
while(loop){
st.v <- sample(1:length(eofv[[v]]),1)
nb <- c(eofv[[v]],eofv[[v]])[st.v:(st.v+4)]
if(length(which(eofv[[nb[2]]]==nb[5])) == 0){
loop <- FALSE
}
}
newv <- 1 + length(eofv)
eofv[[nb[2]]] <- my.replace.nodes(eofv[[nb[2]]],v,c(newv,nb[5],v))
eofv[[nb[3]]] <- my.replace.nodes(eofv[[nb[3]]],v,c(newv))
eofv[[nb[4]]] <- my.replace.nodes(eofv[[nb[4]]],v,c(newv))
eofv[[nb[5]]] <- my.replace.nodes(eofv[[nb[5]]],v,c(v,nb[2],newv))
take1 <- which(eofv[[v]]==nb[3])
take2 <- which(eofv[[v]]==nb[4])
eofv[[v]] <- eofv[[v]][-c(take1,take2)]
eofv[[newv]] <- nb[2:5]
tmp.tri <- sort(c(v,nb[2],nb[3]))
tmp.tri.id <- which(tris[,1]==tmp.tri[1] & tris[,2] ==tmp.tri[2] & tris[,3]==tmp.tri[3])
tmp.col <- cols[tmp.tri.id]
add.tris <- matrix(c(v,nb[2],nb[5],nb[2],nb[5],newv,newv,nb[2],nb[3],newv,nb[4],nb[3],newv,nb[4],nb[5]),byrow=TRUE,ncol=3)
add.tris.ord <- add.tris
tmp.order.tri <- tris.ord[tmp.tri.id,]
tmp.order.tri2 <- c(tmp.order.tri,tmp.order.tri[1])
loc1 <- which(tmp.order.tri==v)
if(tmp.order.tri2[loc1+1]==nb[2]){
new.tris.ord <- rbind(tris.ord,add.tris.ord)
}else{
new.tris.ord <- rbind(tris.ord,cbind(add.tris.ord[,3],add.tris.ord[,2],add.tris.ord[,1]))
}
add.tris <- t(apply(add.tris,1,sort))
add.cols <- c(TRUE,FALSE,TRUE,FALSE,TRUE)
if(!tmp.col){
add.cols <- !add.cols
}
newtris <- rbind(tris,add.tris)
newcols <- c(cols,add.cols)
rmtris1 <- sort(c(v,nb[2],nb[3]))
rmtris2 <- sort(c(v,nb[3],nb[4]))
rmtris3 <- sort(c(v,nb[4],nb[5]))
rmtris1id <- which(newtris[,1]==rmtris1[1] & newtris[,2]==rmtris1[2] & newtris[,3]==rmtris1[3])
rmtris2id <- which(newtris[,1]==rmtris2[1] & newtris[,2]==rmtris2[2] & newtris[,3]==rmtris2[3])
rmtris3id <- which(newtris[,1]==rmtris3[1] & newtris[,2]==rmtris3[2] & newtris[,3]==rmtris3[3])
newtris <- newtris[-c(rmtris1id,rmtris2id,rmtris3id),]
newcols <- newcols[-c(rmtris1id,rmtris2id,rmtris3id)]
new.tris.ord <- new.tris.ord[-c(rmtris1id,rmtris2id,rmtris3id),]
newx <- rbind(x,matrix(0,1,3))
id1 <- which(eofv[[v]]==nb[2])
id2 <- which(eofv[[v]]==nb[3])
id3 <- which(eofv[[v]]==nb[4])
id4 <- which(eofv[[v]]==nb[5])
rest <- eofv[[v]][-(c(id2,id3))]
newx[newv[1],] <- (newx[nb[3],]+newx[nb[3],])/2
newx <- newx/sqrt(apply(newx^2,1,sum))
return(list(eofv=eofv,x=newx,tris=newtris,cols=newcols,tris.ord=new.tris.ord))
}
```
```{r}
my.process.choice <- function(eofv){
current.deg <- sapply(eofv,length)
p.cand <- which(current.deg>=4)
q.cand <- which(current.deg>=6)
s <- sample(1:length(c(p.cand,q.cand)),1)
type <- 0
v <- 0
if(s <= length(p.cand)){
type <- 1
v <- p.cand[s]
}else{
type <- 2
v <- q.cand[s - length(p.cand)]
}
return(list(type=type,v=v))
}
my.replace.nodes <- function(vec,v,u){
loc <- which(vec==v)
pre <- c()
post <- c()
if(loc==1){
post <- vec[-loc]
}else if(loc==length(vec)){
pre <- vec[-loc]
}else{
pre <- vec[1:(loc-1)]
post <- vec[(loc+1):length(vec)]
}
c(pre,u,post)
}
```
初期状態である正八面体から始めて、処理p,qが可能な頂点を選び逐次処理を続けていく。
引数は処理ステップ数、返り値は、三角形の構成頂点IDリスト(順序あり)と三角形の色を基本とし、それ以外のいくつかの情報が返る。
```{r}
my.EulTriSph <- function(n.step){
# エッジリスト
el <- matrix(c(1,2,1,3,1,4,1,5,2,3,3,4,4,5,5,2,2,6,3,6,4,6,5,6),byrow=TRUE,ncol=2)
# 頂点座標
x <- matrix(c(0,0,1,1,0,0,0,1,0,-1,0,0,0,-1,0,0,0,-1),byrow=TRUE,ncol=3)
x <- x/sqrt(apply(x^2,1,sum))
# 隣接頂点、順序考慮のリスト
eofv <- list()
eofv[[1]] <- c(2,3,4,5)
eofv[[2]] <- c(1,5,6,3)
eofv[[3]] <- c(1,2,6,4)
eofv[[4]] <- c(1,3,6,5)
eofv[[5]] <- c(1,4,6,2)
eofv[[6]] <- c(5,4,3,2)
# 三角形構成頂点の行列
# 三角形のノードのIDでソートしておく
tris <- matrix(c(1,2,3,1,3,4,1,4,5,1,5,2,6,2,3,6,3,4,6,4,5,6,5,2),byrow=TRUE,ncol=3)
tris <- t(apply(tris,1,sort))
# 三角形の向きをそろえておく
tris.ord <- matrix(c(1,2,3,1,4,3,1,4,5,1,2,5,6,2,3,6,4,3,6,4,5,6,2,5),byrow=TRUE,ncol=3)
# 三角形の色
cols <- c(FALSE,TRUE,FALSE,TRUE,TRUE,FALSE,TRUE,FALSE)
for(i in 1:n.step){
proc <- my.process.choice(eofv)
if(proc$type==1){
out <- rule.p(eofv,proc$v,x,tris,cols,tris.ord)
eofv <- out[[1]]
x <- out[[2]]
tris <- out[[3]]
cols <- out[[4]]
tris.ord <- out[[5]]
}else{
out <- rule.q(eofv,proc$v,x,tris,cols,tris.ord)
eofv <- out[[1]]
x <- out[[2]]
tris <- out[[3]]
cols <- out[[4]]
tris.ord <- out[[5]]
}
}
return(list(tris.ord=tris.ord,cols=cols,tris=tris,eofv=eofv,x=x))
}
```
ここで得られた三角化情報から有向グラフを作る。
引数は向き考慮の三角形構成頂点行列 tris、三角形の色ベクトル cols。
返り値はエッジリスト。
```{r}
my.dir.graph <- function(tris,cols){
tris.half <- tris[which(cols==0),]
rbind(cbind(tris.half[,1],tris.half[,2]),cbind(tris.half[,2],tris.half[,3]),cbind(tris.half[,3],tris.half[,1]))
}
```
描いてみる。
```{r}
n.step <- 20
euler.out <- my.EulTriSph(n.step)
el <- my.dir.graph(euler.out[[1]],euler.out[[2]])
plot(graph.edgelist(el))
```
球面オイラー三角化が得られたので、適当な辺をルート辺としてを選び、対応する木情報を取得する。
引数は三角形の辺のリストの行列と、三角形の色のベクトルとルート辺(始点・終点頂点からなる長さ2のベクトル。色がFALSEである三角形の辺であり、その三角形の頂点の並び順に合致していること)。
返り値は、２つの辺リストとルート辺。
２つの辺リストは、(1) 木のそれ、(2) 全体のそれ。
少し補足。
オイラー三角化されて2色塗りされた三角形のうち、片方の色の三角形の辺を取り出せば、全ての辺が対象になる。
ある、ルートエッジを採用することで、ルートエッジの始点ノードからのグラフ距離が決まる。
今、片方の色の三角形の三頂点のうち、ルートエッジの始点ノードからの距離が2番目に遠い点と、3番目に遠い点とを採用すると、それが求めるエッジ。
```{r}
my.euler.tree <- function(tris,cols,root.e){
el <- my.dir.graph(tris,cols)
root.st <- root.e[1]
root.end <- root.e[2]
rm.e <- which(el[,1]==root.st | el[,2]==root.st)
g <- graph.edgelist(el)
sh.dist <- shortest.paths(g,root.st,mode="out")
tris.half <- tris[which(cols==0),]
tris.label <- matrix(sh.dist[c(tris.half)],ncol=3)
ord.tris.label <- t(apply(tris.label,1,order))
selected.1 <- t(tris.half)[3*(0:(length(tris.half[,1])-1))+ord.tris.label[,2]]
selected.2 <- t(tris.half)[3*(0:(length(tris.half[,1])-1))+ord.tris.label[,3]]
selected.edges <- cbind(selected.1,selected.2)
selected.edges <- rbind(root.e,selected.edges)
list(tree.edge=selected.edges,whole.edge=el,root.e=root.e)
}
```
実行してみる。
```{r}
root.e <- el[1,]
euler.tree.out <- my.euler.tree(euler.out[[1]],euler.out[[2]],root.e)
g.tree <- graph.edgelist(euler.tree.out[[1]])
is.connected(g.tree)
g.tree
sum(euler.out[[2]])
nrow(euler.out[[1]])
```
また、グラフの辺数は頂点数より１多いことがわかる。これは木グラフの性質である。
また、nrow(euler.out[[1]])は三角形の数である。
そのうち半分がTRUE色の三角形、残りの半分がFALSE色の三角形であるが、それはsum(euler.out[[2]])によっても確認できる。
また、ある色の三角形の個数は、木の辺の数より１少なくなっている（各三角形から１辺を取り出し、ルート辺を加えることで木を作成したから）。
木のラベルにルート頂点からの距離も併せて表示すると、す隣接頂点の距離値は１だけ違うことが確かめられる。
```{r}
sh.dist <- shortest.paths(graph.edgelist(euler.tree.out[[2]]),root.e[1],mode="out")
plot(graph.edgelist(euler.tree.out[[1]]),vertex.label=sh.dist)
```
