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){
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を格納した行列のペアに作り直す
- それを実行するにあたり、ベクトルの一致を探すユーティリティ関数も作っておく
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)]
}
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
- この処理をするにあたり、点に接続する辺、点を持つ三角形、辺が帰属する三角形など、頂点、辺、三角形の相互帰属関係を取り出しておく
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
}
my.tri.vet <- function(tri){
v.of.e <- my.tri.edge(tri)
n.v <- length(unique(c(tri)))
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,])