曲面が酔歩で広がること

• 完全にランダムな酔歩様曲面も作れるが、広がり方にある偏りを入れることで、コーラルリーフ的な広がりが作れたりする

library(geometry)
library(gtools)
library(igraph)

CategoryVector<-function(d){
  df <- d - 1
  diagval <- 1:d
  diagval <- sqrt((d)/df) * sqrt((d - diagval)/(d - diagval + 1))
  others <- -diagval/(d - (1:d))
  m <- matrix(rep(others, d), nrow = d, byrow = TRUE)
  diag(m) <- diagval
  m[upper.tri(m)] <- 0
  as.matrix(m[, 1:df])
}

##########
pmt <- permutations(4,2)
E <- matrix(0,length(pmt[,1]),4)
for(i in 1:length(E[,1])){
	E[i,pmt[i,1]] <- 1
	E[i,pmt[i,2]] <- -1
}

my.make.edgelist <- function(edge){
	n <- length(edge)
	p <- which(edge==1)
	m <- which(edge==-1)
	s <- which(edge==0)
	ret <- list()
	for(i in 1:length(s)){
		tmp1 <- tmp2 <- rep(0,n)
		tmp1[m] <- 1
		tmp1[s[i]] <- -1
		tmp2[p] <- -1
		tmp2[s[i]] <- 1
		ret[[i]] <- rbind(tmp1,tmp2)
	}
	ret
}

my.make.triangle <- function(edge){
	edgelist <- my.make.edgelist(edge)
	ret <- list()
	cnt <- 1
	for(i in 1:length(edgelist)){
		ret[[cnt]] <- rbind(rep(0,length(edge)),edge,edge+edgelist[[i]][1,])
		cnt <- cnt+1
		ret[[cnt]] <- rbind(rep(0,length(edge)),edge,edge+edgelist[[i]][2,])
		cnt <- cnt+1
	}
	ret
}

# 三角形・辺の登録方法は、点をルールでソートしたものとする
my.v.sort <- function(X){
	X[with(as.data.frame(X),order(V1,V2,V3,V4)),]
}



my.check.same.tri <- function(tri,P){
	for(i in 1:length(P[1,1,])){
		tmp <- sum((tri-P[,,i])^2)
		if(tmp == 0){
			return(TRUE)
		}
	}
	return(FALSE)
}
my.check.inside <- function(tri,Insides){
	if(nrow(Insides)==0){
		return(FALSE)
	}
	for(i in 1:nrow(tri)){
		tmp <- apply((t(Insides)-tri[i,])^2,2,sum)
		if(min(tmp)==0){
			return(TRUE)
		}
	}
	return(FALSE)
}
my.check.on.border <- function(tri,B){
	Bs <- rbind(t(B[1,,]),t(B[2,,]))
	for(i in 1:nrow(tri)){
		tmp <- apply((t(Bs)-tri[i,])^2,2,sum)
		if(min(tmp)>0){
			return(FALSE)
		}
	}
	return(TRUE)
}
my.check.used.border <- function(tri,B){
	st <- rbind(tri[1,],tri[1,],tri[2,])
	end <- rbind(tri[2,],tri[3,],tri[3,])
	
	used.id <- rep(0,length(st[,1]))
	used.coords <- array(NA,dim=dim(B))
	B.st <- t(B[1,,])
	B.end <- t(B[2,,])
	
	for(i in 1:length(st[,1])){
		st.check <- apply((t(B.st)-st[i,])^2,2,sum)
		end.check <- apply((t(B.end)-end[i,])^2,2,sum)
		st.end.check <- st.check + end.check
		tmp <- which(st.end.check==0)
		if(length(tmp)>0){
			used.id[i] <- tmp
		}
		used.coords[1,,i] <- st[i,]
		used.coords[2,,i] <- end[i,]
	}
	return(list(used.id=used.id,used.corrds=used.coords))
}
my.check.shared.v <- function(b,tri){
	ret <- rep(0,length(tri[,1]))
	for(i in 1:length(ret)){
		for(j in 1:length(b[1,1,])){
			tmp1 <- sum((b[1,,j]-tri[i,])^2)
			tmp2 <- sum((b[2,,j]-tri[i,])^2)
			if(tmp1*tmp2==0){
				ret[i] <- ret[i]+1
			}
		}
	}
	ret
}
tris <- apply(E,1,my.make.triangle)

Ori <- c(0,0,0,0)

ijk <- sample(1:4,3)

ss <- sample(1:length(tris),1)

sss <- sample(1:4,1)

V <- tris[[ss]][[sss]]

P <- list()
P[[1]] <- my.v.sort(V)

P <-array(V,dim=c(3,4,1))
B <- array(0,dim=c(2,4,3))
B[1,,1] <- V[1,]
B[2,,1] <- V[2,]
B[1,,2] <- V[1,]
B[2,,2] <- V[3,]
B[1,,3] <- V[2,]
B[2,,3] <- V[3,]

P <- array(apply(P,3,my.v.sort),dim=dim(P))
B <- array(apply(B,3,my.v.sort),dim=dim(B))
Insides <- matrix(0,0,4)
n.step <- 500
cnt <- 1
while(cnt < n.step){
#for(i in 1:n.step){
	#print("step:")
	#print(cnt)
	b <- sample(1:length(B[1,1,]),1)
	#t <- sample(1:4,1)
	t <- sample(1:2,1)
	t <- sample(1:4,1,prob=c(0.5,0.45,0.05,0.00))
	tri <- my.make.triangle(B[2,,b]-B[1,,b])[[t]]
	tri2 <- t(t(tri) + B[1,,b])
	tri2.sorted <- my.v.sort(tri2)
	same.tri <- my.check.same.tri(tri2.sorted,P)
	if(same.tri){
		#print("same.tri")
		next
	}
	inside <- my.check.inside(tri2.sorted,Insides)
	if(inside){
		#print("inside")
		next
	}
	on.border <- my.check.on.border(tri2.sorted,B)
	#print(on.border)
	used.border <- my.check.used.border(tri2.sorted,B)
	# 先端点はボーダーに乗っているが点で接している
	#print(used.border[[1]])
	if(on.border & length(which(used.border[[1]]!=0))==1){
		#print("touch with a point")
		next
	}
	# 接続するボーダーの本数を返す
	shared.v <- my.check.shared.v(array(used.border[[2]][,,which(used.border[[1]]!=0)],dim=c(2,4,length(which(used.border[[1]]!=0)))),tri2.sorted)
	r <- runif(1)
	if(length(which(used.border[[1]]!=0))==1){
		if(r>0.1){
			#print("skip due to No. border=1")
			next
		}
	}else if(length(which(used.border[[1]]!=0))==3){
		next
	}
	#print("add")
	cnt <- cnt+1
	# used.borderをBから削り
	# !used.borderをBに加え
	# 三角形をPに加え
	# 複数のボーダーに接続する点をInsidesに加える
	#print(B[,,used.border[[1]]])
	#print(tri2.sorted)
	B <- B[,,-used.border[[1]]]
	B <- abind(B,used.border[[2]][,,which(used.border[[1]]==0)])
	P <- abind(P,tri2.sorted)
	Insides <- rbind(Insides,tri2.sorted[which(shared.v>1),])
	#print(P)
	#print(B)
	#print(Insides)
}

cv <- CategoryVector(4)

Ps <- matrix(0,0,4)
for(i in 1:length(P[1,1,])){
	Ps <- rbind(Ps,P[,,i])
}

Ps.3 <- Ps %*% cv
triangles3d(Ps.3,alpha=0.8,col=sample(c(4,5,6),length(P[,,1]),replace=TRUE))

• Rmdファイル

---
title: "酔歩曲面"
author: "ryamada"
date: "Monday, May 04, 2015"
output: html_document
---

この文書は３次元空間上の正四面体格子にSelf-avoidingにランダムな曲面をシミュレーション生成するための諸関数とその利用例に関するものである。

# Rのパッケージ読み込み
```{r}
library(knitr)
library(geometry)
library(gtools)
library(igraph)
library(rgl)
```

# 関数

```{r}
# 正四面体格子では、エッジは12軸あり、その正負方向を考慮すると、隣接格子点は24個ある
# その12通りの軸ベクトルを4次元表現する
# この基本エッジは、４成分のうちの２成分は0、残りの２成分は±1である
my.make.E <- function(){
  pmt <- permutations(4,2)
  E <- matrix(0,length(pmt[,1]),4)
  for(i in 1:length(E[,1])){
    E[i,pmt[i,1]] <- 1
    E[i,pmt[i,2]] <- -1
  }
  E
}

# 基本エッジを引数とする三角形は４つある
# 三角形の３辺ベクトルの和は(0,0,0,0)になるので
# 引数エッジと、返り値のエッジとの和は、原点からの単位エッジの位置に
# 来るはずである
# それを満足するような引数エッジに接続するエッジが４種類あり、それを返す
# それを作成するには、1,2,3,4から3つを取り、そのi,j,kの順列を考えることになる
# 今、４つのうちの２つは引数エッジの非ゼロ要素で決まっているので、kの取り方は２通り
# ijkの並べ方は順逆の２通りなので、計４通り
my.make.edgelist <- function(edge){
	n <- length(edge)
	p <- which(edge==1)
	m <- which(edge==-1)
	s <- which(edge==0)
	ret <- list()
	for(i in 1:length(s)){
		tmp1 <- tmp2 <- rep(0,n)
		tmp1[m] <- 1
		tmp1[s[i]] <- -1
		tmp2[p] <- -1
		tmp2[s[i]] <- 1
		ret[[i]] <- rbind(tmp1,tmp2)
	}
	ret
}
# my.make.edgelist()の出力を使い、実際に、原点(0,0,0,0)始まりで、引数edgeが
# 第２頂点を定めるとして、第３頂点座標を定める
# そのようにしてできる三角形の頂点座標を3x4行列で作成し
# その3x4行列を要素とする長さ４のリストを返す
my.make.triangle <- function(edge){
	edgelist <- my.make.edgelist(edge)
	ret <- list()
	cnt <- 1
	for(i in 1:length(edgelist)){
		ret[[cnt]] <- rbind(rep(0,length(edge)),edge,edge+edgelist[[i]][1,])
		cnt <- cnt+1
		ret[[cnt]] <- rbind(rep(0,length(edge)),edge,edge+edgelist[[i]][2,])
		cnt <- cnt+1
	}
	ret
}

# 三角形・辺の登録方法は、点をルールでソートしたものとする
# 4成分の第１成分を最優先、第４成分を最も後回しとして、全４列でソートする
my.v.sort <- function(X){
	X[with(as.data.frame(X),order(V1,V2,V3,V4)),]
}
# 三角形triに一致する三角形が、三角形のリストであるアレイオブジェクトPの中に
# 含まれるかどうかの判定
# 三角形は3x4行列で表す。３行はmy.v.sort()のルールでソートしてあるので
# 三角形の一致は、3x4成分のすべての一致によって確認できる
my.check.same.tri <- function(tri,P){
	for(i in 1:length(P[1,1,])){
		tmp <- sum((tri-P[,,i])^2)
		if(tmp == 0){
			return(TRUE)
		}
	}
	return(FALSE)
}
# 三角形を表すtriの頂点のうち一つでも、Insidesという頂点座標を納めた行列の中の座標と
# 一致すればTRUEを、そうでなければFALSEを返す
# この関数は、ボーダー上の辺からできた三角形をtriとして受け取るので
# 少なくとも、ボーダー辺にある２頂点はInsidesには含まれないから
# 結局、第３頂点がInsidesに含まれるかの判定となっている
my.check.inside <- function(tri,Insides){
	if(nrow(Insides)==0){
		return(FALSE)
	}
	for(i in 1:nrow(tri)){
		tmp <- apply((t(Insides)-tri[i,])^2,2,sum)
		if(min(tmp)==0){
			return(TRUE)
		}
	}
	return(FALSE)
}
# 作成した三角形の第３頂点が、ボーダー頂点であるかどうかを返す
# 実際、triの２頂点は必ずボーダー頂点であるので、一つでも、ボーダー頂点でない
# 頂点があれば、第３頂点はボーダー頂点ではない、という判定をしている
my.check.on.border <- function(tri,B){
	Bs <- rbind(t(B[1,,]),t(B[2,,]))
	for(i in 1:nrow(tri)){
		tmp <- apply((t(Bs)-tri[i,])^2,2,sum)
		if(min(tmp)>0){
			return(FALSE)
		}
	}
	return(TRUE)
}
# できた三角形の３辺が、ボーダー辺であるかどうかを判定する
# 1辺は必ず、ボーダー辺であり
# 残りの２辺はボーダー辺かもしれないし、そうでないかもしれない
# 3辺ともにボーダー辺であると、四面体が閉じてしまうので、そのような場合は
# 三角形を不採用にするためにも使う
# 実際には、Bの何行目のボーダーに一致したかの情報と
# その辺の始点終点とを返す(my.v.sort()の順を守る)
my.check.used.border <- function(tri,B){
	st <- rbind(tri[1,],tri[1,],tri[2,])
	end <- rbind(tri[2,],tri[3,],tri[3,])
	
	used.id <- rep(0,length(st[,1]))
	used.coords <- array(NA,dim=dim(B))
	B.st <- t(B[1,,])
	B.end <- t(B[2,,])
	
	for(i in 1:length(st[,1])){
		st.check <- apply((t(B.st)-st[i,])^2,2,sum)
		end.check <- apply((t(B.end)-end[i,])^2,2,sum)
		st.end.check <- st.check + end.check
		tmp <- which(st.end.check==0)
		if(length(tmp)>0){
			used.id[i] <- tmp
		}
		used.coords[1,,i] <- st[i,]
		used.coords[2,,i] <- end[i,]
	}
	return(list(used.id=used.id,used.corrds=used.coords))
}
# 三角形の３辺のうち、ボーダー辺が複数あるときに
# 三角形の頂点がいくつのボーダー辺に接続するかを数える
# ０個のボーダー辺に帰属する頂点は、ただの新規頂点
# １個のボーダー辺に帰属する場合は、内部化しないで、ボーダー辺上の点である
# ２個のボーダー辺に帰属する場合には、内在化する
# ３個のボーダー辺に帰属する場合には、四面体が閉じる
my.check.shared.v <- function(b,tri){
	ret <- rep(0,length(tri[,1]))
	for(i in 1:length(ret)){
		for(j in 1:length(b[1,1,])){
			tmp1 <- sum((b[1,,j]-tri[i,])^2)
			tmp2 <- sum((b[2,,j]-tri[i,])^2)
			if(tmp1*tmp2==0){
				ret[i] <- ret[i]+1
			}
		}
	}
	ret
}
# 正四面体の頂点座標を作る関数
# 格子座標を４次元表現して諸処理をするが、最後に３次元座標に変換するために使用する
CategoryVector<-function(d){
  df <- d - 1
  diagval <- 1:d
  diagval <- sqrt((d)/df) * sqrt((d - diagval)/(d - diagval + 1))
  others <- -diagval/(d - (1:d))
  m <- matrix(rep(others, d), nrow = d, byrow = TRUE)
  diag(m) <- diagval
  m[upper.tri(m)] <- 0
  as.matrix(m[, 1:df])
}
# 原点を含む１個の三角形を初期値として作る
my.initial.tri <- function(rand=FALSE){
  E <- my.make.E()
  tris <- apply(E,1,my.make.triangle)
  if(rand){
    ss <- sample(1:length(tris),1)
    sss <- sample(1:4,1)
  }else{
    ss <- 1
    sss <- 1
  }
  V <- tris[[ss]][[sss]]
  P <-array(V,dim=c(3,4,1))
  B <- array(0,dim=c(2,4,3))
  B[1,,1] <- V[1,]
  B[2,,1] <- V[2,]
  B[1,,2] <- V[1,]
  B[2,,2] <- V[3,]
  B[1,,3] <- V[2,]
  B[2,,3] <- V[3,]
    
  P <- array(apply(P,3,my.v.sort),dim=dim(P))
  B <- array(apply(B,3,my.v.sort),dim=dim(B))
  Insides <- matrix(0,0,4)
  return(list(P=P,B=B,Insides=Insides))
}
```

# シミュレーションする関数
```{r}
# 初期状態は原点を含む三角形１つ
# n.stepは１個の三角形から始める場合には、最終三角形の個数
# rは延長三角形の採用条件を決めるベクトル
# 延長三角形が含むボーダー辺の数1,2,3に応じて、三角形を採用する確率をr[1],r[2],r[3]とする
# r[3]=0は、四面体を閉じさせない条件
# r[1]を小さくし、r[2]=1とすれば、比較的曲面に切れ込みが少なくなる
# probは、選択したエッジから伸ばせる４つの三角形の採用確率を決めるベクトル
# 均等条件(0.25,0.25,0.25,0.25)からずらすと特定の方向への伸びが高確率となる
my.random.surface <- function(PBI=my.initial.tri(),n.step=100,r=c(1,1,0),prob=c(0.25,0.25,0.25,0.25)){
  P <- PBI[[1]]
  B <- PBI[[2]]
  Insides <- PBI[[3]]
  cnt <- 1
  while(cnt < n.step){
  #for(i in 1:n.step){
    #print("step:")
  	#print(cnt)
  	b <- sample(1:length(B[1,1,]),1)
  	#t <- sample(1:4,1)
  	#t <- sample(1:2,1)
  	t <- sample(1:4,1,prob=prob)
  	tri <- my.make.triangle(B[2,,b]-B[1,,b])[[t]]
  	tri2 <- t(t(tri) + B[1,,b])
  	tri2.sorted <- my.v.sort(tri2)
  	same.tri <- my.check.same.tri(tri2.sorted,P)
  	if(same.tri){
  		#print("same.tri")
  		next
  	}
  	inside <- my.check.inside(tri2.sorted,Insides)
  	if(inside){
  		#print("inside")
  		next
  	}
  	on.border <- my.check.on.border(tri2.sorted,B)
  	#print(on.border)
  	used.border <- my.check.used.border(tri2.sorted,B)
  	# 先端点はボーダーに乗っているが点で接している
  	#print(used.border[[1]])
  	if(on.border & length(which(used.border[[1]]!=0))==1){
  		#print("touch with a point")
  		next
  	}
  	# 接続するボーダーの本数を返す
  	shared.v <- my.check.shared.v(array(used.border[[2]][,,which(used.border[[1]]!=0)],dim=c(2,4,length(which(used.border[[1]]!=0)))),tri2.sorted)
  	R <- runif(1)
  	#if(length(which(used.border[[1]]!=0))==1){
  		if(R>r[length(which(used.border[[1]]!=0))]){
  			#print("skip due to No. border=1")
  			next
  		}
  	#}
    #else if(length(which(used.border[[1]]!=0))==2){
    #  if(R>r[2]){
    #    next
    #  }
  	#}else if(length(which(used.border[[1]]!=0))==3){
    #  if(R>r[3]){
    #    next
    #  }
  	#}
  	#print("add")
  	cnt <- cnt+1
  	# used.borderをBから削り
  	# !used.borderをBに加え
  	# 三角形をPに加え
  	# 複数のボーダーに接続する点をInsidesに加える
  	#print(B[,,used.border[[1]]])
  	#print(tri2.sorted)
  	B <- B[,,-used.border[[1]]]
  	B <- abind(B,used.border[[2]][,,which(used.border[[1]]==0)])
  	P <- abind(P,tri2.sorted)
  	Insides <- rbind(Insides,tri2.sorted[which(shared.v>1),])
  	#print(P)
  	#print(B)
  	#print(Insides)
  }
  
  cv <- CategoryVector(4)
  
  Ps <- matrix(0,0,4)
  for(i in 1:length(P[1,1,])){
  	Ps <- rbind(Ps,P[,,i])
  }
  
  Ps.3 <- Ps %*% cv
  return(list(P=P,B=B,Insides=Insides,Ps=Ps,Ps.3=Ps.3,PBI=PBI,prob=prob,r=r))
}
```
```{r setup}
knit_hooks$set(rgl = hook_rgl)
```
まずはデフォルト設定で、まったくのランダムな条件で実行する。

```{r,rgl=TRUE}
test.out <- my.random.surface()
triangles3d(test.out$Ps.3,alpha=0.8,col=sample(c(4,5,6),length(test.out$P[,,1]),replace=TRUE))
```

次に、ボーダーと２辺で接する三角形を優先して伸ばす条件で実行する
```{r,rgl=TRUE}
test.out <- my.random.surface(r=c(0.1,1,0))
plot3d(test.out$Ps.3)
triangles3d(test.out$Ps.3,alpha=0.8,col=sample(c(4,5,6),length(test.out$P[,,1]),replace=TRUE))
```

次に、ボーダーとの辺共有の条件は元に戻し、
```{r,rgl=TRUE}
test.out <- my.random.surface(r=c(0.1,1,0))
plot3d(test.out$Ps.3)
triangles3d(test.out$Ps.3,alpha=0.8,col=sample(c(4,5,6),length(test.out$P[,,1]),replace=TRUE))
```

```{r,rgl=TRUE}
test.out <- my.random.surface(r=c(0.1,1,0),prob=c(1,1,0.1,0),n.step=500)
plot3d(test.out$Ps.3)
triangles3d(test.out$Ps.3,alpha=0.8,col=sample(c(4,5,6),length(test.out$P[,,1]),replace=TRUE))
```

```{r}
my.random.surface.2 <- function(PBI=my.initial.tri(),n.step=100,r=c(1,1,0),prob=c(0.25,0.25,0.25,0.25)){
  P <- PBI[[1]]
  B <- PBI[[2]]
  Insides <- PBI[[3]]
  cnt <- 1
  while(cnt < n.step){
  #for(i in 1:n.step){
    #print("step:")
    #print(cnt)
    b <- sample(1:length(B[1,1,]),1,prob=(length(B[1,1,]):1)^2)
  	#t <- sample(1:4,1)
  	#t <- sample(1:2,1)
  	t <- sample(1:4,1,prob=prob)
  	tri <- my.make.triangle(B[2,,b]-B[1,,b])[[t]]
  	tri2 <- t(t(tri) + B[1,,b])
  	tri2.sorted <- my.v.sort(tri2)
  	same.tri <- my.check.same.tri(tri2.sorted,P)
  	if(same.tri){
  		#print("same.tri")
  		next
  	}
  	inside <- my.check.inside(tri2.sorted,Insides)
  	if(inside){
  		#print("inside")
  		next
  	}
  	on.border <- my.check.on.border(tri2.sorted,B)
  	#print(on.border)
  	used.border <- my.check.used.border(tri2.sorted,B)
  	# 先端点はボーダーに乗っているが点で接している
  	#print(used.border[[1]])
  	if(on.border & length(which(used.border[[1]]!=0))==1){
  		#print("touch with a point")
  		next
  	}
  	# 接続するボーダーの本数を返す
  	shared.v <- my.check.shared.v(array(used.border[[2]][,,which(used.border[[1]]!=0)],dim=c(2,4,length(which(used.border[[1]]!=0)))),tri2.sorted)
  	R <- runif(1)
  	#if(length(which(used.border[[1]]!=0))==1){
  		if(R>r[length(which(used.border[[1]]!=0))]){
  			#print("skip due to No. border=1")
  			next
  		}
  	#}
    #else if(length(which(used.border[[1]]!=0))==2){
    #  if(R>r[2]){
    #    next
    #  }
  	#}else if(length(which(used.border[[1]]!=0))==3){
    #  if(R>r[3]){
    #    next
    #  }
  	#}
  	#print("add")
  	cnt <- cnt+1
  	# used.borderをBから削り
  	# !used.borderをBに加え
  	# 三角形をPに加え
  	# 複数のボーダーに接続する点をInsidesに加える
  	#print(B[,,used.border[[1]]])
  	#print(tri2.sorted)
  	B <- B[,,-used.border[[1]]]
  	B <- abind(B,used.border[[2]][,,which(used.border[[1]]==0)])
  	P <- abind(P,tri2.sorted)
  	Insides <- rbind(Insides,tri2.sorted[which(shared.v>1),])
  	#print(P)
  	#print(B)
  	#print(Insides)
  }
  
  cv <- CategoryVector(4)
  
  Ps <- matrix(0,0,4)
  for(i in 1:length(P[1,1,])){
  	Ps <- rbind(Ps,P[,,i])
  }
  
  Ps.3 <- Ps %*% cv
  return(list(P=P,B=B,Insides=Insides,Ps=Ps,Ps.3=Ps.3,PBI=PBI,prob=prob,r=r))
}
```
```{r,rgl=TRUE}
test.out <- my.random.surface.2(r=c(1,1,0),n.step=30)
plot3d(test.out$Ps.3)
triangles3d(test.out$Ps.3,alpha=0.8,col=sample(c(4,5,6),length(test.out$P[,,1]),replace=TRUE))
```
```{r,rgl=TRUE}
test.out2 <- my.random.surface.2(PBI=list(test.out$P,test.out$B,test.out$Insides),r=c(0.1,1,0),prob=c(0.25,0.25,0.25,0.25),n.step=500)
plot3d(test.out2$Ps.3)
triangles3d(test.out2$Ps.3,alpha=0.8,col=sample(c(4,5,6),length(test.out2$P[,,1]),replace=TRUE))
```