SHARE
TWEET

Untitled

a guest Oct 21st, 2019 71 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. library(geometry) # convexhull, delaunay triangulation
  2. library(Directional) # von-Mises
  3. library(igraph)
  4. library(rgl)
  5. library(threejs)
  6. library(dirichletprocess)
  7. library(MCMCpack)
  8. library(mvtnorm)
  9.  
  10. my.runifS2 <- function(n){
  11.     X <- matrix(rnorm(n*3),ncol=3)
  12.     X <- X/sqrt(apply(X^2,1,sum))
  13.     return(X)
  14. }
  15. # こちらは、接点pが(0,0,1)とみなし、二次元平面上の点xをxy平面上の点とし
  16. # 球面に写像した上で
  17. # (0,0,1)を接点pに移す大円回転でxの対応点も移す
  18. my.Tp2S2 <- function(x,p){
  19.   tmp.p <- c(0,0,1)
  20.  
  21.   # 回転させたいが、まず、pのz = 1平面上の対応点をとる
  22.   p. <- my.Riemannian.log(p,tmp.p)
  23.   p. <- c(p.[1:2],0)
  24.   # その方向の単位ベクトルにする
  25.   p. <- p./sqrt(sum(p.^2))
  26.   # p.と直交するベクトルを取る
  27.   q. <- c(p.[2],-p.[1],0)
  28.  
  29.   # z = 1 平面上で (x-tmp.p)ベクトルを、(p.-tmp.p)と(q.-tmp.p)の線形和で表すこととし
  30.   # その係数を求める
  31.   M <- cbind(p.[1:2],q.[1:2])
  32.   coefs <- solve(M) %*% x
  33.   # tmp.pでの接平面をpでの接平面に移す
  34.   # ベクトル p. - tmp. p が p.' - p に移るとき
  35.   # p.' - pは以下に代わる
  36.   # theta はtmp.pとpとの角
  37.   theta <- acos(sum(tmp.p*p))
  38.   p.._p <- c(cos(theta)*p.[1:2],-sin(theta))
  39.   # q. - tmp.p = q.' - pとなる
  40.   q.._p <- q.
  41.   #M. <- cbind(p.._p,q.._q)
  42.   x._p <- coefs[1] * p.._p + coefs[2] * q.._p
  43.   x. <- c(x._p) + c(p)
  44.   #return(x.)
  45.   return(my.Riemannian.Exp(x._p,p))
  46. }
  47. # pが接点
  48. # x が球面上の点
  49. # 関数の帰り値にpを加えた座標が切平面に乗る
  50. my.Riemannian.log <- function(x,p){
  51.   cos.t <- sum(x*p)
  52.   if(abs(cos.t)>1){
  53.     cos.t <- sign(cos.t) * 1
  54.   }
  55.   t <- acos(cos.t)
  56.   return((x-p*cos.t)*t/sin(t))
  57. }
  58. # 上記 log関数の逆関数
  59. my.Riemannian.Exp <- function(x,p){
  60.   r <- sqrt(sum(x^2))
  61.   ret <- p * cos(r) + x/r*sin(r)
  62.   return(ret)
  63. }
  64. my.DPMMS2 <- function(N,alpha=10,df=3,m=matrix(c(1,0,0,1),2,2)*0.01){
  65.   s <- StickBreaking(alpha,N)
  66.   t <- rmultinom(1,N,prob=s)
  67.   mus <- my.runifS2(N)
  68.   sigmas <- list()
  69.  
  70.   for(i in 1:N){
  71.     sigmas[[i]] <- rwish(df,m)
  72.   }
  73.   x <- matrix(0,N,3)
  74.   cnt <- 1
  75.   for(i in 1:length(t)){
  76.     if(t[i]!=0){
  77.         x2d <- matrix(rmvnorm(t[i],c(0,0),sigmas[[i]]),ncol=2)
  78.       x3d <- t(apply(x2d,1,my.Tp2S2,mus[i,]))
  79.       x[cnt:(cnt+t[i]-1),] <- x3d
  80.       cnt <- cnt + t[i]
  81.     }
  82.   }
  83.   return(list(x=x,s=s,t=t,mus=mus,sigmas=sigmas))
  84. }
  85. my.plot3d.dpmms2 <- function(out,r=0.05){
  86.   x <- out$x
  87.   t <- out$t
  88.   plot3d(x)
  89.   cnt <- 1
  90.   for(i in 1:length(t)){
  91.     if(t[i]!=0){
  92.         spheres3d(x[cnt:(cnt+t[i]-1),],color=i,radius=0.05)
  93.     cnt <- cnt + t[i]
  94.     }
  95.   }
  96. }
  97.  
  98. N <- 100
  99. dpmms2out <- my.DPMMS2(N)
  100. plot3d(dpmms2out$x)
  101. my.plot3d.dpmms2(dpmms2out)
  102.  
  103.  
  104. X <- dpmms2out$x
  105.  
  106. # 凸包
  107. hull <- convhulln(X) # hullは三角形の頂点IDの3列行列
  108.  
  109. # エッジ
  110. ed <- rbind(hull[,1:2],hull[,2:3],hull[,c(3,1)])
  111. ed <- unique(t(apply(ed,1,sort)))
  112.  
  113. plot3d(X)
  114. segments3d(X[t(ed),])
  115.  
  116. g <- graph.edgelist(ed,directed=FALSE)
  117. adj <- as.matrix(get.adjacency(g))
  118. gx <- layout_with_fr(g,dim=3,niter=2500)
  119.  
  120. dmat <- as.matrix(dist(gx))
  121.  
  122. hist(dmat[which(adj==1)])
  123.  
  124. graphjs(g,gx)
  125.  
  126. ###########
  127. # 正三角形でできた立体たち
  128.  
  129. # 正多面体3種
  130. tetra <- rbind(c(1,2,3),c(1,3,4),c(1,4,2),c(3,2,4))
  131. 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))
  132. 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))
  133.  
  134. # strictly-convex delta-hedron 5種
  135.  
  136. 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))
  137. 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))
  138. # decaは5角錐の張り合わせ。dodecaはその張り合わせ部分に2つの三角形を割り込ませたもの
  139. 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))
  140. # 正三角柱(側面は正方形)の側面に正四角錐を張り付け
  141. 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))
  142. # 正四角錐で正四角反柱をサンドイッチ
  143. 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))
  144.  
  145. # 半正多面体のうち、辺数=3k、頂点数=2kのものの双対は三角メッシュになる
  146. # すべての頂点は合同なので、頂点次数は3であることになるから
  147.  
  148. # 半正多面体
  149. 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))
  150. 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))
  151. 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))
  152. 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))
  153.  
  154. 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))
  155.  
  156. 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))
  157.  
  158. # truncated icosidodecahedron
  159.  
  160. # アルキメデスの角柱(プリズム)
  161. # 上面・下面がk角形であり、側面が正方形であるような角柱を考える
  162. # 3正則グラフなので、この双対グラフとしての三角メッシュを作る関数を作る
  163. my.prism.tri <- function(k){
  164.     top =1
  165.     bottom = k+2
  166.     sides = 2:(k+1)
  167.     sides <- c(sides,2)
  168.     ret <- matrix(0,0,3)
  169.     for(i in 1:k){
  170.         tmp1 <- c(top,sides[i],sides[i+1])
  171.         tmp2 <- c(bottom,sides[i+1],sides[i])
  172.         ret <- rbind(ret,tmp1,tmp2)
  173.     }
  174.     return(ret)
  175. }
  176. my.prism.tri(4)
  177.  
  178. # 上面・下面が3角形・6角形のそれぞれのアンチプリズム
  179. antiprism3 <- octa # 正八面体のこと
  180. 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))
  181.  
  182.  
  183. hull <- heccaideca
  184. # エッジ
  185. ed <- rbind(hull[,1:2],hull[,2:3],hull[,c(3,1)])
  186. ed <- unique(t(apply(ed,1,sort)))
  187.  
  188.  
  189. g <- graph.edgelist(ed,directed=FALSE)
  190. adj <- as.matrix(get.adjacency(g))
  191. gx <- layout_with_fr(g,dim=3,niter=500,start.temp=sqrt(vcount(g))*0.1)
  192.  
  193. dmat <- as.matrix(dist(gx))
  194.  
  195. plot(sort(dmat[which(adj==1)]))
  196.  
  197. graphjs(g,gx)
  198.  
  199. graphjs(g,
  200.   layout=list(
  201.     layout_randomly(g, dim=3),
  202.     layout_on_sphere(g),
  203.     layout_with_drl(g, dim=3),  # note! somewhat slow...
  204.     layout_with_fr(g, dim=3, niter=30)),
  205.   main=list("random layout", "sphere layout", "drl layout", "fr layout"),
  206.   fpl=300)
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top