3D Voxelデータから平滑化三角メッシュデータ

  • これがデモ
library(devtools)
install_github("ryamada22/Ronlyryamada")
library(Ronlyryamada)
example(my.catmull.clark.tri)
  • まずは、適当に3次元立体オブジェクトを作り、3次元アレイに0,1のデータとして与えよう
    • 適当なアレイ領域にいくつかの球のOR結合として作成する
my.3d.obj <- function(n,k=5,a=6,b=4){
	x <- 1:(2^k)
	xyz <- as.matrix(expand.grid(x,x,x))
	xyz.val <- rep(0,length(xyz[,1]))
	rs <- runif(n) * (max(x)/a)
	ctrs <- matrix(runif(n*3),ncol=3) * max(x)/b + max(x)*(b-1)/(b*2)
	for(i in 1:n){
		tmp <- which(apply((t(xyz) - ctrs[i,])^2,2,sum) < rs[i]^2)
		xyz.val[tmp] <- 1
	}
	return(list(xyz=xyz,xyz.val=xyz.val))
}
k <- 4
out <- my.3d.obj(3,k)
obj <- out$xyz[which(out$xyz.val==1),]
library(rgl)
plot3d(obj)
xyz.arr <- array(out$xyz.val,c(2^k,2^k,2^k))
  • このvoxelアレイ xyz.arrの輪郭を取り出して三角形メッシュ化する
    • その処理においては、marching cubeというアルゴリズムを用いる
    • Rではmisc3dパッケージのcomputeContour3d()関数があるのでそれを使う
    • 以下の関数では、computeContour3d()関数が返すつぶれた三角形(3点が同一の三角形や、長さが0の辺を含む三角形を除いている
    • この出力オブジェクトは、三角形の頂点座標を格納している
library(misc3d)
my.peri.tri <- function(v,h=max(v),smooth=0){
	#library(misc3d)
	con <- computeContour3d(v, h, 1)
	tris <- makeTriangles(con,smooth=smooth)

	V1 <- tris$v1
	V2 <- tris$v2
	V3 <- tris$v3

	L12 <- apply((V1-V2)^2,1,sum)
	L13 <- apply((V1-V3)^2,1,sum)
	L23 <- apply((V2-V3)^2,1,sum)
	L123 <- L12*L13*L23
	zeros <- (L123==0)

	tris$v1 <- V1[!zeros,]
	tris$v2 <- V2[!zeros,]
	tris$v3 <- V3[!zeros,]
	tris	
}
tris <- my.peri.tri(xyz.arr)
drawScene.rgl(tris)
  • これを少し使いやすいように、頂点座標を格納した行列と、三角形の3頂点IDを格納した行列のペアに作り直す
    • それを実行するにあたり、ベクトルの一致を探すユーティリティ関数も作っておく
#' v,uは行がベクトル
#' uの行ベクトルがvの行ベクトルの何行目かを返す
my.tri.v.id <- function(v,u){
	d <- matrix(0,length(v[,1]),length(u[,1]))
	for(i in 1:length(v[1,])){
		d <- d+abs(outer(v[,i],u[,i],"-"))
	}
	ret <- which(d==0,arr.ind=TRUE)
	ret[,c(2,1)]
}
# triangleオブジェクトを引数にとり
# 頂点にIDを振り、その座標vと
# 三角形の頂点座標セットと
# 三角形の頂点座標を納めたリストを返す

my.tri.vid <- function(tris){
	v <- unique(rbind(tris$v1,tris$v2,tris$v3))
	v1 <- my.tri.v.id(v,tris$v1)
	v2 <- my.tri.v.id(v,tris$v2)
	v3 <- my.tri.v.id(v,tris$v3)
	tri.vid <- cbind(v1[,2],v2[,2],v3[,2])
	tri.vid <- t(apply(tri.vid,1,sort))
	return(list(v=v,tri.vid=tri.vid,tri.coord=list(tris$v1,tris$v2,tris$v3)))
}
v <- matrix(c(1,2,3,2,3,1,3,3,3),byrow=TRUE,ncol=3)

u <- v[c(3,2,1,3,2,2,2,3),]
my.tri.v.id(v,u)
ttt <- my.tri.vid(tris)
  • さて、でこぼごが激しいので、これを平滑化する
    • 平滑化にあたっては、Catmull-Clark subdivisionを採用する
      • 簡単に言うと、各三角形を6つの三角形に分割しながら平滑化する方法
      • 三角形の重心を囲む6個の三角形に分割する
      • その6個を作るにあたり、3頂点を平滑化にあたって移動した3点と、3辺の「中央的な点」との計6点が重心点をぐるりと囲むようにして、そのぐるりを6個の三角形にする
      • 辺の中央的な点は、辺の両側の三角形の重心点と、辺の両端の2点との計4点の重心とする
      • 三角形の3点の移動は次のようにする
        • もとの頂点の座標P
        • もとの頂点を囲む三角形の重心の平均座標F(この三角形の個数をnとする)
        • もとの頂点に接続する辺(辺の数もnである)の「中央的な点」の平均座標R
        • \frac{F+2R+(n-2)P}{n}
    • この処理をするにあたり、点に接続する辺、点を持つ三角形、辺が帰属する三角形など、頂点、辺、三角形の相互帰属関係を取り出しておく
# 三角形の頂点IDから、ノード、エッジ、三角形の相互帰属関係を返す

# 三角形頂点からエッジリストを返す
my.tri.edge <- function(tri){
	ed <- rbind(tri[,1:2],tri[,2:3],tri[,c(1,3)])
	ed <- t(apply(ed,1,sort))
	ed <- unique(ed)
	ed
}
# 頂点 v、辺 e、三角形 vの相互の帰属関係を算出する
my.tri.vet <- function(tri){
	v.of.e <- my.tri.edge(tri)
	n.v <- length(unique(c(tri)))
	#v <- 1:length((x[,1]))
	v.of.t <- t(apply(tri,1,sort))
	e.of.v <- t.of.v <- list()
	e.of.t <- matrix(0,length(tri[,1]),3)
	t.of.e <- matrix(0,length(v.of.e[,1]),2)
	for(i in 1:n.v){
		e.of.v[[i]] <- sort(c(which(v.of.e[,1]==i),which(v.of.e[,2]==i)))
		t.of.v[[i]] <- sort(c(which(v.of.t[,1]==i),which(v.of.t[,2]==i),which(v.of.t[,3]==i)))
	}
	tmp.et.mat <- matrix(0,length(v.of.e[,1]),length(tri[,1]))
	for(i in 1:length(v.of.e[,1])){
		tmp1 <- t(tri[,1:2])-v.of.e[i,]
		tmp2 <- t(tri[,c(1,3)])-v.of.e[i,]
		tmp3 <- t(tri[,2:3])-v.of.e[i,]
		tmp123 <- apply(abs(tmp1),2,sum)*apply(abs(tmp2),2,sum)*apply(abs(tmp3),2,sum)
		tmp4 <- which(tmp123==0)
		tmp.et.mat[i,tmp4] <- 1		
	}
	for(i in 1:length(tmp.et.mat[,1])){
		t.of.e[i,] <- which(tmp.et.mat[i,]==1)
	}
	for(i in 1:length(tmp.et.mat[1,])){
		e.of.t[i,] <- which(tmp.et.mat[,i]==1)
	}
	return(list(e.of.v=e.of.v,t.of.v=t.of.v,v.of.e=v.of.e,t.of.e=t.of.e,v.of.t=v.of.t,e.of.t=e.of.t))
}
sss <- my.tri.vet(ttt[[2]])
    • そのうえで、Catmull-Clark subdivisionをする関数を書く
      • ちょっと遅いのだが…
my.catmull.clark.tri <- function(x.vet){
	x <- x.vet[[1]]
	v.of.t <- x.vet[[2]]$v.of.t
	t.of.e <- x.vet[[2]]$t.of.e
	v.of.e <- x.vet[[2]]$v.of.e
	face.pt <- (x[v.of.t[,1],] + x[v.of.t[,2],] + x[v.of.t[,3],])/3
	edge.pt <- (face.pt[t.of.e[,1],] + face.pt[t.of.e[,2],] + x[v.of.e[,1],] + x[v.of.e[,2],])/4
	new.v.pt <- matrix(0,length(x[,1]),3)
	for(i in 1:length(new.v.pt[,1])){
		ori.v <- x[i,]
		fs <- x.vet[[2]]$t.of.v[[i]]
		n <- length(fs)
		mean.f.pt <- apply(face.pt[fs,],2,mean)
		ed <- x.vet[[2]]$e.of.v[[i]]
		mean.e.pt <- apply(edge.pt[ed,],2,mean)
		new.v.pt[i,] <- (mean.f.pt + 2*mean.e.pt + (n-3)*ori.v)/n
	}
	new.v.x <- rbind(new.v.pt,face.pt,edge.pt)
	n.v <- length(new.v.pt[,1])
	n.f <- length(face.pt[,1])
	
	new.tris <- matrix(0,n.f*6,3)
	for(i in 1:length(v.of.t[,1])){
		e.id <- x.vet[[2]]$e.of.t[i,]
		v.of.e.of.t <- v.of.e[e.id,]
		
		new.f.id <- i + n.v
		new.e.id <- e.id + n.v + n.f
		
		tmp <- matrix(0,6,3)
		tmp[1:2,1] <- v.of.e.of.t[1,]
		tmp[3:4,1] <- v.of.e.of.t[2,]
		tmp[5:6,1] <- v.of.e.of.t[3,]
		tmp[1:2,2] <- new.e.id[1]
		tmp[3:4,2] <- new.e.id[2]
		tmp[5:6,2] <- new.e.id[3]
		tmp[,3] <- new.f.id
		
		new.tris[(1:6) + (i-1)*6,] <- tmp
	}
	new.tris <- t(apply(new.tris,1,sort))
	new.vet <- my.tri.vet(new.tris)
	return(list(new.v.x,new.vet))
}
x.vet <- list(ttt[[1]],sss)

new.vet <- my.catmull.clark.tri(x.vet)
plot3d(new.vet[[1]][new.vet[[2]]$v.of.t,])

x.vet <- new.vet

new.vet <- my.catmull.clark.tri(new.vet)
plot3d(new.vet[[1]][new.vet[[2]]$v.of.t,])