#' Euler Triangulation on Sphere
#'
#' A set of functions to randomly generate euler triangulation on Sphere.
#' @examples
#' library(igraph)
#' n.step <- 20
#' euler.out <- my.EulTriSph(n.step)
#' el <- my.dir.graph(euler.out[[1]],euler.out[[2]])
#' plot(graph.edgelist(el))
#' 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]])
#' 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)
#' @import igraph
#' @export
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){
# p/q処理の選択と対象頂点を選択
proc <- my.process.choice(eofv)
if(proc$type==1){# p処理
# rule p
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{# q処理
# rule q
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))
}
#' @export
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]))
}
# 三角形構成頂点ID(値ソート)、色ベクトルの更新
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[v,] <- (apply(matrix(x[c(id1,id2,id3),],ncol=3),2,mean) + x[v,])/2
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))
}
#' @export
rule.q <- function(eofv,v,x,tris,cols,tris.ord){
loop <- TRUE
# 周囲ノードでnb[2],nb[5]として取り出された2頂点間にすでに辺があるときは
# 同一の処理ができないのでそうならない場合を取り出す
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))
#rest <- eofv[[v]][-((st.v+1):(st.v+4))]
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[v,] <- (apply(matrix(x[rest,],ncol=3),2,mean) + x[v,])/2
#newx[newv[1],] <- (apply(x[nb[3:4],],2,mean) + x[v,])/2
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))
}
#' @export
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))
}
#' @export
my.replace.nodes <- function(vec,v,u){
# 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)
}
#' @export
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]))
}
#' @export
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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.