library(geometry) # convexhull, delaunay triangulation
library(Directional) # von-Mises
library(igraph)
library(rgl)
library(threejs)
library(dirichletprocess)
library(MCMCpack)
library(mvtnorm)

my.runifS2 <- function(n){
	X <- matrix(rnorm(n*3),ncol=3)
	X <- X/sqrt(apply(X^2,1,sum))
	return(X)
}
# こちらは、接点pが(0,0,1)とみなし、二次元平面上の点xをxy平面上の点とし
# 球面に写像した上で
# (0,0,1)を接点pに移す大円回転でxの対応点も移す
my.Tp2S2 <- function(x,p){
  tmp.p <- c(0,0,1)

  # 回転させたいが、まず、pのz = 1平面上の対応点をとる
  p. <- my.Riemannian.log(p,tmp.p)
  p. <- c(p.[1:2],0)
  # その方向の単位ベクトルにする
  p. <- p./sqrt(sum(p.^2))
  # p.と直交するベクトルを取る
  q. <- c(p.[2],-p.[1],0)

  # z = 1 平面上で (x-tmp.p)ベクトルを、(p.-tmp.p)と(q.-tmp.p)の線形和で表すこととし
  # その係数を求める
  M <- cbind(p.[1:2],q.[1:2])
  coefs <- solve(M) %*% x
  # tmp.pでの接平面をpでの接平面に移す
  # ベクトル p. - tmp. p が p.' - p に移るとき
  # p.' - pは以下に代わる
  # theta はtmp.pとpとの角
  theta <- acos(sum(tmp.p*p))
  p.._p <- c(cos(theta)*p.[1:2],-sin(theta))
  # q. - tmp.p = q.' - pとなる
  q.._p <- q.
  #M. <- cbind(p.._p,q.._q)
  x._p <- coefs[1] * p.._p + coefs[2] * q.._p
  x. <- c(x._p) + c(p)
  #return(x.)
  return(my.Riemannian.Exp(x._p,p))
}
# pが接点
# x が球面上の点
# 関数の帰り値にpを加えた座標が切平面に乗る
my.Riemannian.log <- function(x,p){
  cos.t <- sum(x*p)
  if(abs(cos.t)>1){
    cos.t <- sign(cos.t) * 1
  }
  t <- acos(cos.t)
  return((x-p*cos.t)*t/sin(t))
}
# 上記 log関数の逆関数
my.Riemannian.Exp <- function(x,p){
  r <- sqrt(sum(x^2))
  ret <- p * cos(r) + x/r*sin(r)
  return(ret)
}
my.DPMMS2 <- function(N,alpha=10,df=3,m=matrix(c(1,0,0,1),2,2)*0.01){
  s <- StickBreaking(alpha,N)
  t <- rmultinom(1,N,prob=s)
  mus <- my.runifS2(N)
  sigmas <- list()

  for(i in 1:N){
    sigmas[[i]] <- rwish(df,m)
  }
  x <- matrix(0,N,3)
  cnt <- 1
  for(i in 1:length(t)){
    if(t[i]!=0){
        x2d <- matrix(rmvnorm(t[i],c(0,0),sigmas[[i]]),ncol=2)
      x3d <- t(apply(x2d,1,my.Tp2S2,mus[i,]))
      x[cnt:(cnt+t[i]-1),] <- x3d
      cnt <- cnt + t[i]
    }
  }
  return(list(x=x,s=s,t=t,mus=mus,sigmas=sigmas))
}
my.plot3d.dpmms2 <- function(out,r=0.05){
  x <- out$x
  t <- out$t
  plot3d(x)
  cnt <- 1
  for(i in 1:length(t)){
    if(t[i]!=0){
        spheres3d(x[cnt:(cnt+t[i]-1),],color=i,radius=0.05)
    cnt <- cnt + t[i]
    }
  }
}

N <- 100
dpmms2out <- my.DPMMS2(N)
plot3d(dpmms2out$x)
my.plot3d.dpmms2(dpmms2out)


X <- dpmms2out$x

# 凸包
hull <- convhulln(X) # hullは三角形の頂点IDの３列行列

# エッジ
ed <- rbind(hull[,1:2],hull[,2:3],hull[,c(3,1)])
ed <- unique(t(apply(ed,1,sort)))

plot3d(X)
segments3d(X[t(ed),])

g <- graph.edgelist(ed,directed=FALSE)
adj <- as.matrix(get.adjacency(g))
gx <- layout_with_fr(g,dim=3,niter=2500)

dmat <- as.matrix(dist(gx))

hist(dmat[which(adj==1)])

graphjs(g,gx)

###########
# 正三角形でできた立体たち

# 正多面体３種
tetra <- rbind(c(1,2,3),c(1,3,4),c(1,4,2),c(3,2,4))
octa <- rbind(c(1,2,3),c(1,3,4),c(1,4,5),c(1,5,2),c(6,3,2),c(6,4,3),c(6,5,4),c(6,2,5))
icosa <- rbind(c(1,2,3),c(1,3,4),c(1,4,5),c(1,5,6),c(1,6,2),c(3,2,7),c(4,3,8),c(5,4,9),c(6,5,10),c(2,6,11),c(3,7,8),c(4,8,9),c(5,9,10),c(6,10,11),c(2,11,7),c(12,8,7),c(12,9,8),c(12,10,9),c(12,11,10),c(12,7,11))

# strictly-convex delta-hedron ５種

hexa <- rbind(c(1,2,3),c(1,3,4),c(1,4,2),c(5,3,2),c(5,4,3),c(5,2,4))
deca <- rbind(c(1,2,3),c(1,3,4),c(1,4,5),c(1,5,6),c(1,6,2),c(7,3,2),c(7,4,3),c(7,5,4),c(7,6,5),c(7,2,6))
# decaは５角錐の張り合わせ。dodecaはその張り合わせ部分に２つの三角形を割り込ませたもの
dodeca <- rbind(c(1,2,3),c(1,3,4),c(1,4,5),c(1,5,6),c(1,6,2),c(7,3,2),c(7,4,3),c(7,5,4),c(7,8,5),c(7,2,8),c(2,6,8),c(8,6,5))
# 正三角柱(側面は正方形)の側面に正四角錐を張り付け
tetracaideca <- rbind(c(1,2,3),c(7,8,9),c(4,2,1),c(4,1,7),c(4,7,8),c(4,8,2),c(5,3,2),c(5,2,8),c(5,8,9),c(5,9,3),c(6,1,3),c(6,3,9),c(6,9,7),c(6,7,1))
# 正四角錐で正四角反柱をサンドイッチ
heccaideca <- rbind(c(1,2,3),c(1,3,4),c(1,4,5),c(1,5,2),c(10,7,6),c(10,8,7),c(10,9,8),c(10,6,9),c(2,6,3),c(3,7,4),c(4,8,5),c(5,9,2),c(3,6,7),c(4,7,8),c(5,8,9),c(2,9,6))

# 半正多面体のうち、辺数=3k、頂点数=2kのものの双対は三角メッシュになる
# すべての頂点は合同なので、頂点次数は３であることになるから

# 半正多面体
trunc.tetra <- rbind(c(1,2,3),c(1,3,4),c(1,4,2),c(2,5,3),c(3,6,4),c(4,7,2),c(8,5,2),c(8,3,5),c(8,6,3),c(8,4,6),c(8,7,4),c(8,2,7))
trunc.cube <- rbind(c(1,2,6),c(1,6,3),c(1,3,7),c(1,7,4),c(1,4,8),c(1,8,5),c(1,5,9),c(1,9,2),c(3,6,7),c(4,7,8),c(5,8,9),c(2,9,6),c(6,11,7),c(7,12,8),c(8,13,9),c(9,10,6),c(6,14,11),c(11,14,7),c(7,14,12),c(12,14,8),c(8,14,13),c(13,14,9),c(9,14,10),c(10,14,8))
trunc.octa <- rbind(c(1,2,6),c(1,6,3),c(1,3,7),c(1,7,4),c(1,4,5),c(1,5,2),c(6,2,11),c(6,11,9),c(6,9,12),c(6,12,3),c(7,3,12),c(7,12,10),c(7,10,13),c(7,13,4),c(5,4,13),c(5,13,8),c(5,8,11),c(5,11,2),c(9,14,12),c(12,14,10),c(10,14,13),c(13,14,8),c(8,14,11),c(11,14,9))
trunc.dodeca <- rbind(c(1,2,7),c(1,7,3),c(1,3,8),c(1,8,4),c(1,4,9),c(1,9,5),c(1,5,10),c(1,10,6),c(1,6,11),c(1,11,2),c(2,11,7),c(3,7,8),c(4,8,9),c(5,9,10),c(6,10,11),c(7,12,8),c(8,14,9),c(9,16,10),c(10,18,11),c(11,20,7),c(7,20,26),c(7,26,21),c(7,21,22),c(7,22,12),c(8,12,22),c(8,22,13),c(8,13,23),c(8,23,14),c(9,14,23),c(9,23,15),c(9,15,24),c(9,24,16),c(10,16,22),c(10,22,17),c(10,17,25),c(10,25,18),c(11,18,25),c(11,25,19),c(11,19,26),c(11,26,20),c(21,26,22),c(13,22,23),c(15,23,24),c(17,24,25),c(19,25,26),c(22,26,31),c(22,31,32),c(22,32,27),c(22,27,23),c(23,27,32),c(23,32,28),c(23,28,24),c(24,28,32),c(24,32,29),c(24,20,25),c(25,29,32),c(25,32,30),c(25,30,26),c(26,30,32),c(26,32,31))

trunc.icosa <- rbind(c(1,7,4),c(4,7,8),c(1,9,7),c(1,2,9),c(2,23,9),c(2,3,23),c(3,29,23),c(3,19,29),c(3,13,19),c(2,13,3),c(2,10,13),c(1,10,2),c(1,5,10),c(1,4,5),c(4,6,5),c(4,8,6),c(6,8,17),c(8,26,17),c(17,26,31),c(17,31,16),c(6,17,16),c(11,6,16),c(5,6,11),c(5,11,12),c(10,5,12),c(10,12,13),c(13,12,14),c(13,14,19),c(8,21,26),c(8,7,21),c(7,22,21),c(7,9,22),c(9,23,22),c(22,23,25),c(23,29,25),c(16,31,32),c(16,32,15),c(16,15,11),c(11,15,12),c(12,15,14),c(15,32,18),c(15,18,14),c(14,18,19),c(20,19,18),c(29,19,20),c(26,27,31),c(27,26,24),c(26,21,24),c(21,22,24),c(24,22,25),c(27,24,28),c(24,25,28),c(25,29,28),c(28,29,20),c(31,27,30),c(27,28,30),c(28,20,30),c(32,31,30),c(18,32,30),c(20,18,30))

trunc.cubocuta <- rbind(c(1,2,3),c(1,3,4),c(1,4,5),c(1,5,2),c(2,5,9),c(2,9,8),c(2,8,6),c(2,6,3),c(3,6,10),c(3,10,16),c(3,16,11),c(3,11,7),c(3,7,4),c(4,7,12),c(4,12,8),c(4,8,5),c(5,8,13),c(5,13,19),c(5,19,14),c(5,14,9),c(6,15,10),c(7,11,12),c(8,12,13),c(9,14,15),c(10,15,21),c(10,21,22),c(10,22,16),c(16,22,11),c(11,22,17),c(11,17,12),c(12,17,23),c(12,23,18),c(12,18,13),c(13,18,24),c(13,24,19),c(19,24,14),c(14,24,20),c(14,20,15),c(15,20,25),c(15,25,21),c(21,25,22),c(22,25,26),c(22,26,23),c(22,23,17),c(23,26,24),c(23,24,18),c(24,26,25),c(20,24,25))

# truncated icosidodecahedron

# アルキメデスの角柱(プリズム)
# 上面・下面がk角形であり、側面が正方形であるような角柱を考える
# ３正則グラフなので、この双対グラフとしての三角メッシュを作る関数を作る
my.prism.tri <- function(k){
	top =1
	bottom = k+2
	sides = 2:(k+1)
	sides <- c(sides,2)
	ret <- matrix(0,0,3)
	for(i in 1:k){
		tmp1 <- c(top,sides[i],sides[i+1])
		tmp2 <- c(bottom,sides[i+1],sides[i])
		ret <- rbind(ret,tmp1,tmp2)
	}
	return(ret)
}
my.prism.tri(4)

# 上面・下面が３角形・６角形のそれぞれのアンチプリズム
antiprism3 <- octa # 正八面体のこと
antiprism6 <- rbind(c(1,2,3),c(1,3,4),c(1,4,5),c(1,5,6),c(1,6,7),c(1,7,2),c(14,9,8),c(14,10,9),c(14,11,10),c(14,12,11),c(14,13,12),c(14,8,13),c(2,8,3),c(3,9,4),c(4,10,5),c(5,11,6),c(6,12,7),c(7,13,2),c(3,8,9),c(4,9,10),c(5,10,11),c(6,11,12),c(7,12,13),c(2,13,8))


hull <- heccaideca
# エッジ
ed <- rbind(hull[,1:2],hull[,2:3],hull[,c(3,1)])
ed <- unique(t(apply(ed,1,sort)))


g <- graph.edgelist(ed,directed=FALSE)
adj <- as.matrix(get.adjacency(g))
gx <- layout_with_fr(g,dim=3,niter=500,start.temp=sqrt(vcount(g))*0.1)

dmat <- as.matrix(dist(gx))

plot(sort(dmat[which(adj==1)]))

graphjs(g,gx)

graphjs(g,
  layout=list(
    layout_randomly(g, dim=3),
    layout_on_sphere(g),
    layout_with_drl(g, dim=3),  # note! somewhat slow...
    layout_with_fr(g, dim=3, niter=30)),
  main=list("random layout", "sphere layout", "drl layout", "fr layout"),
  fpl=300)