Nothing
#
# EM algorithm for Mixture of Multivariate Canonical Fundamental Skew t-distributioins
# Package: EMMIX-cskew
# Version: 0.9-1
#
# Code by S. X. Lee
# Updated on 07 Sep 2015
#
# Lee S.X. and Mclachlan, G.J. (2015) Finite mixtures of canonical fundamental
# skew t-distributions: the unification of the restricted and unrestricted
# skew t-mixture models. Statistics and Computing. doi:10.1007/s11222-015-9545-x
#
################################################################################
# contour.r
# contour plots in 2D and 3D
#
################################################################################
#2D Contours
fmcfust.contour.2d <- function(dat, model, grid=50, drawpoints=TRUE, clusters=NULL,
nlevels=10, map=c("scatter", "heat", "cluster"),
component=NULL, xlim, ylim, xlab, ylab, main, pcol=NULL, ccol=NULL, ...) {
levels <- nlevels
if(length(levels)==1) {
if(levels<0) stop("'levels' cannot have negative value")
if(is.whole(levels)) lflag<- T
else if(levels>1) stop("'levels' must be either an integer or a quantile value") else lflag<-F
}else if (all(levels > 0) && all(levels<1)) lflag <- F
else stop("'levels' must be either an integer or a vector of quantile values")
if(ncol(dat)<2) stop("Data must be 2D!")
if(ncol(dat)!=2) warning("Data should be 2-D. Only the first two columns will be used.")
if(is.null(clusters)) clusters <- rep(1, nrow(dat))
dat<-as.matrix(dat)[,1:2]
rx<-range(dat[,1])+c(-1,1)
ry<-range(dat[,2])+c(-1,1)
if (missing(xlim)) xlim<-rx
if (missing(ylim)) ylim<-ry
g<-length(model$dof) ;p<-2
if (missing(xlab)) xlab<-"X"
if (missing(ylab)) ylab<-"Y"
if (missing(main)) {
msg <- bquote(paste("Multivariate CFUST mixture Model (g = ", .(g), ")", sep=""))
} else msg=main
x <- seq(xlim[1], xlim[2], length=grid)
y <- seq(ylim[1], ylim[2], length=grid)
nx <- length(x)
ny <- length(y)
xoy <- cbind(rep(x,ny), as.vector(matrix(y,nx,ny,byrow=TRUE)))
X <- matrix(xoy, nx*ny, 2, byrow=FALSE)
ColClust <- c("lightyellow", "lightgreen", "lightblue3", "lightpink", "lightgrey",
"lightgreen", "lightsalmon", colors()[sample(1:length(colors()),g)])
ColPoint <- ColCont <- c("blue", "red", "black", "Brown", "green", colors()[sample(1:length(colors()),g)])
if(!is.null(pcol)) ColPoint <- array(pcol,g)
if(!is.null(ccol)) ColCont <- array(ccol,g)
if(map[1]=="heat") {
smoothScatter(dat, nrpoints=Inf, colramp=colorRampPalette(heatColor()),xlab=xlab,ylab=ylab,main=msg, ...)
drawpoints <- FALSE
}else if(map[1]=="cluster") {
cat("creating cluster map\n")
clustergrid <- cfustDA(g, X, model)
plot(x = 0, y = 0,type = "n", xlim = xlim, ylim = ylim,xlab=xlab,ylab=ylab, main=msg)
for(i in 1:g) points(X[clustergrid==i,], pch=22, lwd=grid/5, col=ColClust[i], bg=ColClust[i])
}else plot(x = 0, y = 0,type = "n", xlim = xlim, ylim = ylim, xlab=xlab, ylab=ylab, main=msg)
if(drawpoints) {for(i in 1:g) points(dat[clusters==i,],pch=19,col=ColPoint[i])}
if(is.null(component)) {
dens <- dfmcfust(X, model$mu, model$sigma, model$delta, model$dof, model$pro)
densMat <- matrix(dens,nx,ny)
if(lflag) contour(x, y, densMat, nlevels=levels, add=TRUE, drawlabels=FALSE, lty=1, col = ColCont[1],...)
else {
levels <- quantile(densMat, probs=levels)
contour(x, y, densMat, levels=levels, add=TRUE, drawlabels=FALSE, lty=1, col = ColCont[1],...)
}
} else {
ncomp <- length(as.numeric(component))
if(!all(component %in% c(1:g))) stop("argument 'component' is incorrectly specified!")
if(!lflag) levels <- quantile(densMat, probs=levels)
for(i in 1:ncomp) {
dens <- dcfust(X, model$mu[[component[i]]],
model$sigma[[component[i]]], model$delta[[component[i]]],
model$dof[component[i]])
densMat <- matrix(dens,nx,ny)
if(lflag) contour(x, y, densMat, nlevels=levels, add=TRUE,
drawlabels=FALSE, lty=1,col = ColCont[i],...)
else contour(x, y, densMat, levels=levels, add=TRUE,
drawlabels=FALSE, lty=1,col = ColCont[i],...)
}
}
}
cfustDA <- function(g, dat, model) {
if (g == 1) return(rep(1, nrow(dat)))
probMat <- matrix(0,g,nrow(dat))
for(i in 1:g) probMat[i,] <- model$pro[i]*dcfust(dat, model$mu[[i]], model$sigma[[i]], model$delta[[i]], model$dof[i])
clust <- apply(probMat,2,which.max)
return(clust)
}
#3D Contours
fmcfust.contour.3d <- function(dat, model, grid=20, drawpoints=TRUE, levels=0.9,
clusters=NULL, xlim, ylim, zlim, xlab, ylab, zlab, main, component=NULL, pcol=NULL, ccol=NULL,...) {
if(length(levels)==1) {
if(levels<0) stop("'levels' cannot have negative value")
if(is.whole(levels)) lflag<- T
else {if(levels>1) stop("'levels' must be either an integer or a quantile value") else lflag<-F}
}else if (all(levels > 0) && all(levels<1)) lflag <- F
else stop("'levels' must be either an integer or a vector of quantile values")
if(missing(dat)) {
if(drawpoints) stop("dat must be provided if drawpoints=TRUE.")
if(missing(xlim) || missing(ylim) || missing(zlim)) stop("x, y, z limits must be provided if dat is not provided.")
} else {
if(ncol(dat)<3) stop("Data must be 3D!")
if(ncol(dat)!=3) warning("Data should be 3-D. Only the first three columns will be used.")
if(is.null(clusters)) clusters <- rep(1, nrow(dat))
dat<-as.matrix(dat)[,1:3]
rx<-range(dat[,1])+c(-1,1)
ry<-range(dat[,2])+c(-1,1)
rz<-range(dat[,3])+c(-1,1)
if (missing(xlim)) xlim<-rx
if (missing(ylim)) ylim<-ry
if (missing(zlim)) zlim<-rz
}
g<-length(model$dof) ;p<-3
if (missing(xlab)) xlab<-"X"
if (missing(ylab)) ylab<-"Y"
if (missing(zlab)) zlab<-"Z"
#library(rgl)
ColPoint <- ColCont <- c("magenta", "green", "red", "blue", "grey", colors()[sample(1:length(colors()),g)])
if(!is.null(pcol)) ColPoint <- array(pcol,g)
if(!is.null(ccol)) ColCont <- array(ccol,g)
if(drawpoints) {
tdat <- dat; tclusters <- clusters #cut points that are out of range
Ind <- which(tdat[,1]>=xlim[1] & tdat[,1]<=xlim[2]); tdat <- tdat[Ind,]; tclusters <- tclusters[Ind]
Ind <- which(tdat[,2]>=ylim[1] & tdat[,2]<=ylim[2]); tdat <- tdat[Ind,]; tclusters <- tclusters[Ind]
Ind <- which(tdat[,3]>=zlim[1] & tdat[,3]<=zlim[2]); tdat <- tdat[Ind,]; tclusters <- tclusters[Ind]
plot3d(tdat[,1], tdat[,2], tdat[,3], col=ColPoint[tclusters], xlab=xlab, ylab=ylab, zlab=zlab, xlim=xlim, ylim=ylim, zlim=zlim)
# plot3d(tdat[tclusters==1,1], tdat[tclusters==1,2], tdat[tclusters==1,3],
# col=ColPoint[1], xlab=xlab, ylab=ylab, zlab=zlab, xlim=xlim, ylim=ylim, zlim=zlim)
# if(g>=2){
# for(i in 2:g)
# points3d(tdat[tclusters==i,1], tdat[tclusters==i,2], tdat[tclusters==i,3], col=ColPoint[i])
# }
} else {
plot3d(1,1,1,col=rgb(1,1,1,alpha=1), xlab=xlab, ylab=ylab, zlab=zlab, xlim=xlim, ylim=ylim, zlim=zlim)
}
XVec <- seq(xlim[1], xlim[2], length=grid)
YVec <- seq(ylim[1], ylim[2], length=grid)
ZVec <- seq(zlim[1], zlim[2], length=grid)
if(is.null(component)) {
mst.contour <- function(x, y, z) dfmcfust(cbind(x,y,z), model$mu,
model$sigma, model$delta, model$dof, model$pro)
contour3d(mst.contour, levels, XVec, YVec, ZVec, add=T, color=ColCont[1], fill=T, alpha=0.25, nflag=lflag)
} else{
ncomp <- length(as.numeric(component))
if(!all(component %in% c(1:g))) stop("argument 'component' is incorrectly specified!")
if(length(levels==1)) levels <- rep(levels,ncomp)
for(i in 1:ncomp){
mst.contour <- function(x, y, z) dcfust(cbind(x,y,z), model$mu[[component[i]]],
model$sigma[[component[i]]], model$delta[[component[i]]], model$dof[component[i]])
cat("Plotting component ", component[i] ,"\n")
contour3d(mst.contour, levels[i], XVec, YVec, ZVec, add=T, color=ColCont[component[i]], fill=T, alpha=0.25, nflag=lflag)
}
}
}
################################################################################
#other functions
heatColor <- function(){
rgb(c(255,255,254,254,253,252,227,189,128),
c(255,237,217,178,141,78,26,0,0),
c(204,160,118,76,60,42,28,38,38),maxColorValue=255)
}
#The following code is adapted from the package misc3d (with minor changes)
contour3d <- function(f, level, x = 1:dim(f)[1], y = 1:dim(f)[2], z = 1:dim(f)[3],
mask = NULL, color = "white", color2 = NA, alpha = 1, fill = TRUE,
col.mesh = if (fill) NA else color, material = "default", smooth = 4,
add = FALSE, draw = TRUE, engine = "rgl", separate=FALSE, nflag=F, ...){
if (! all(is.finite(x), is.finite(y), is.finite(z)))
stop("'x', 'y', and 'z' values must be finite and non-missing")
if (is.function(f) || is.array(f) && length(dim(f)) == 3){
if (is.function(f)){
if (length(formals(f)) < 3)
stop("The function must have at least 3 arguments.")
vol <- fgrid(f, x, y, z)
}else{
if (dim(f)[1] != length(x) || dim(f)[2] != length(y) || dim(f)[3] != length(z))
stop("dimensions of f do not match x, y, or z")
vol <- f
}
maxvol <- max(vol)
minvol <- min(vol)
if(missing(level)) {quan<-quantile(vol); level <- quan[4];}
else {
if(nflag) level <- quantile(vol, probs=seq(0,1,length.out=level))
else level <- quantile(vol, probs=level)
}
#cat("lflag is ", nflag, "\n")
#cat("Your quantile levels are ", level, "\n")
con <- which(! level <= maxvol & level >= minvol)
if (length(con) == length(level))
stop(paste("The 'level' has to be within the range of 'f' (between ", minvol, " and ", maxvol,").\n", sep=""))
else if (length(con) > 0){
warning(paste("The 'level' outside the range of 'f' (between ", minvol, " and ", maxvol, ") has been removed. \n", sep=""))
level <- level[-con]
if (is.list(mask)) mask <- mask[-con]
if (length(color) > 1) color <- color[-con]
if (length(color2) > 1) color2 <- color2[-con]
if (length(alpha) > 1) alpha <- alpha[-con]
if (length(fill) > 1) fill <- fill[-con]
if (length(col.mesh) > 1) col.mesh <- col.mesh[-con]
if (length(material) > 1) material <- material[-con]
if (length(smooth) > 1) smooth <- smooth[-con]
}
}
else stop("vol has to be a function or a 3-dimensional array")
scene <- contourTriangles(vol, maxvol, level, x, y, z, mask, color, color2,
alpha, fill, col.mesh, material, smooth)
if (! draw || engine == "none"){
if (! any(separate))
scene
else{
if (length(level)==1){
newScene <- separateTriangles(scene)
cat("Triangles are separated into ", length(newScene),
" chunks.", "\n", sep="")
}
else{
if (length(separate) < length(level))
separate <- c(separate, rep(FALSE, length(level)-length(separate)))
newScene <- NULL
for (i in 1:length(level)){
if (separate[i]){
new <- separateTriangles(scene[[i]])
newScene <- c(newScene, new)
cat("Triangles from level ", level[i],
" are separated into ", length(new), " chunks.",
"\n", sep="")
}
else
newScene <- c(newScene, list(scene[[i]]))
}
}
newScene
}
}
else {
scene <- colorScene(scene)
if (engine == "rgl")
drawScene.rgl(scene, add = add, ...)
else if (engine %in% c("standard", "grid"))
drawScene(scene, add = add, engine = engine, ...)
else stop(paste("unknown rendering engine:", engine))
}
}
PreProcessing <- local({
explode <- function(x)
floor(((x - 1) %% 2^(1:8))/2^(0:7))
BasicRotation <-
matrix(c(1,2,3,4,5,6,7,8,5,6,2,1,8,7,3,4,8,7,6,5,4,3,2,1,
4,3,7,8,1,2,6,5,2,6,7,3,1,5,8,4,6,5,8,7,2,1,4,3,
5,1,4,8,6,2,3,7,4,1,2,3,8,5,6,7,3,4,1,2,7,8,5,6,
2,3,4,1,6,7,8,5,6,7,3,2,5,8,4,1,7,8,4,3,6,5,1,2,
8,5,1,4,7,6,2,3,7,3,2,6,8,4,1,5,4,8,5,1,3,7,6,2,
3,2,6,7,4,1,5,8,2,1,5,6,3,4,8,7,1,4,8,5,2,3,7,6,
1,5,6,2,4,8,7,3,5,8,7,6,1,4,3,2,8,4,3,7,5,1,2,6,
3,7,8,4,2,6,5,1,7,6,5,8,3,2,1,4,6,2,1,5,7,3,4,8),
ncol=8, byrow=TRUE)
CaseRotation <-
matrix(c(1,24,2,19,2,17,3,17,2,24,4,24,3,24,6,10,2,15,3,19,
4,17,6,9,3,9,6,8,6,1,9,23,2,20,3,18,4,7,6,16,5,24,
7,5,7,24,12,9,4,20,6,22,8,24,10,24,7,9,15,24,13,20,
6,20,2,21,4,6,3,16,6,4,4,16,8,23,6,14,10,23,5,21,
7,10,7,16,15,9,7,2,13,8,12,23,6,6,3,6,6,17,6,18,
9,18,7,4,13,17,15,18,6,13,7,6,12,16,13,18,6,2,11,24,
7,3,7,12,3,12,2,23,5,23,4,23,7,1,3,14,7,14,6,21,
15,23,4,15,7,19,8,19,13,23,6,11,12,17,10,19,6,23,4,12,
7,18,8,22,13,16,7,13,11,23,13,21,7,15,8,21,13,22,14,24,
8,15,13,11,7,7,8,12,4,22,3,23,7,23,6,24,12,18,6,7,
13,19,9,24,6,19,7,21,11,18,13,24,7,20,15,16,7,22,6,15,
3,22,6,3,15,17,10,22,6,12,12,24,7,11,6,5,3,15,13,10,
7,8,8,20,4,9,7,17,5,22,4,18,2,22,2,22,4,18,5,22,
7,17,4,9,8,20,7,8,13,10,3,15,6,5,7,11,12,24,6,12,
10,22,15,17,6,3,3,22,6,15,7,22,15,16,7,20,13,24,11,18,
7,21,6,19,9,24,13,19,6,7,12,18,6,24,7,23,3,23,4,22,
8,12,7,7,13,11,8,15,14,24,13,22,8,21,7,15,13,21,11,23,
7,13,13,16,8,22,7,18,4,12,6,23,10,19,12,17,6,11,13,23,
8,19,7,19,4,15,15,23,6,21,7,14,3,14,7,1,4,23,5,23,
2,23,3,12,7,12,7,3,11,24,6,2,13,18,12,16,7,6,6,13,
15,18,13,17,7,4,9,18,6,18,6,17,3,6,6,6,12,23,13,8,
7,2,15,9,7,16,7,10,5,21,10,23,6,14,8,23,4,16,6,4,
3,16,4,6,2,21,6,20,13,20,15,24,7,9,10,24,8,24,6,22,
4,20,12,9,7,24,7,5,5,24,6,16,4,7,3,18,2,20,9,23,
6,1,6,8,3,9,6,9,4,17,3,19,2,15,6,10,3,24,4,24,
2,24,3,17,2,17,2,19,1,24),
ncol=2,byrow=TRUE)
CaseRotationFlip <-
cbind(CaseRotation,
c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,-1,1,1,1,1,1,1,1,1,1,
1,1,-1,1,-1,-1,-1,1,1,1,1,1,1,1,-1,1,1,1,1,1,1,-1,-1,1,1,
1,1,1,1,1,-1,1,1,1,-1,1,-1,-1,-1,1,1,1,1,1,1,1,-1,1,1,1,
-1,1,-1,-1,-1,1,1,1,1,1,1,1,-1,1,1,-1,-1,1,-1,-1,-1,1,1,1,1,
1,-1,1,-1,1,1,1,-1,-1,-1,-1,-1,1,1,-1,-1,1,-1,-1,-1,-1,-1,-1,
-1,-1,
-1,-1,-1,1,1,1,1,1,1,1,1,1,1,1,-1,1,1,-1,-1,1,1,1,1,1,-1,
-1,-1,1,-1,1,-1,-1,-1,-1,-1,1,1,1,-1,1,1,-1,-1,1,-1,-1,-1,-1,
-1,-1,
-1,1,1,1,-1,1,-1,-1,-1,1,-1,-1,-1,-1,-1,-1,-1,1,1,1,-1,1,-1,
-1,-1,
1,-1,-1,-1,-1,-1,-1,-1,1,1,-1,-1,-1,-1,-1,-1,1,-1,-1,-1,-1,-1,
-1,-1,1,
1,1,-1,1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,1,-1,-1,-1,-1,-1,-1,
-1,-1,-1,
-1,-1,-1,-1,-1,-1))
EdgePoints <-
matrix(c(1,1,2,2,2,3,3,3,4,4,4,1,5,5,6,6,6,7,
7,7,8,8,8,5,9,1,5,10,2,6,11,3,7,12,4,8,13,9,9),
ncol=3, byrow=TRUE)
BasicEdges <- list(c(1,4,9),
c(2,4,9,10),
c(1,4,5,6,9,10),
c(1,4,6,7,9,11),
c(1,4,10,11,12),
c(2,4,6,7,9,10,11),
c(1,2,5,6,7,8,9,10,11,13),
c(9,10,11,12),
c(1,2,7,8,9,11),
c(1,2,3,4,5,6,7,8,13),
c(1,2,6,7,9,12),
c(1,4,5,8,9,10,11,12,13),
c(1,2,3,4,5,6,7,8,9,10,11,12,13),
c(1,4,7,8,10,11))
EdgeSequence1 <-
list(c(1,2,3),
c(1,2,4,2,3,4),
list(c(1,2,5,3,4,6),
c(1,2,6,2,4,6,2,3,4,2,5,3)),
list(c(1,2,5,3,4,6),
c(1,3,5,3,4,5,2,5,4,1,2,6,1,6,3,2,4,6)),
c(1,3,2,2,3,5,3,4,5),
list(c(1,2,6,2,5,6,3,4,7),
c(1,2,7,2,4,7,2,5,4,3,4,5,3,5,6,1,6,3,1,3,7,7,3,4),
c(1,2,7,2,4,7,2,5,4,3,4,5,3,5,6)),
list(c(1,8,2,3,7,6,4,5,9),
c(1,8,2,3,9,4,3,7,9,5,9,7,5,7,6),
c(3,7,6,1,8,4,1,4,5,1,9,2,1,5,9),
c(4,5,9,1,7,6,1,6,2,2,6,3,2,3,8),
c(1,10,2,2,10,9,5,9,10,5,10,6,6,10,7,3,7,10,3,10,4,4,10,8,1,
8,10),
c(1,10,2,1,10,7,6,10,7,5,10,6,5,9,10,4,9,10,3,10,4,3,8,10,2,
10,8),
c(5,9,10,2,10,9,1,10,2,1,7,10,6,10,7,3,10,6,3,8,10,4,10,8,4,
5,10),
c(1,7,6,1,6,9,1,9,2,5,9,6,3,8,4),
c(1,7,8,3,8,7,3,7,6,3,6,5,3,5,4,4,5,9,4,9,8,2,8,9,1,8,2)),
c(1,2,3,1,3,4),
c(1,2,6,1,3,5,1,6,3,3,4,5),
list(c(1,4,5,2,6,3,3,6,7,4,8,5),
c(1,4,3,1,3,2,3,4,8,3,8,7,5,8,6,6,8,7,1,2,6,1,6,5),
c(1,2,9,1,9,5,5,9,8,4,8,9,3,4,9,3,9,7,6,7,9,2,6,9),
c(5,9,6,1,9,5,1,4,9,4,8,9,7,9,8,3,9,7,2,9,3,2,6,9),
c(1,2,5,2,6,5,3,4,8,3,8,7)),
c(1,2,3,1,3,6,1,6,5,3,4,6),
list(c(1,6,2,2,6,8,3,5,4,6,7,8),
c(1,5,3,1,3,6,3,6,7,3,7,4,2,4,8,2,5,4,1,5,2,4,7,8),
c(1,9,2,2,9,5,3,5,9,4,9,8,3,9,4,7,8,9,6,7,9,1,9,6),
c(4,9,5,1,5,9,1,9,2,2,9,8,7,8,9,6,9,7,3,6,9,3,9,4),
c(1,5,2,3,8,4,3,6,8,6,7,8)),
list(
##13.1
c(7,8,12,2,3,11,1,4,9,5,6,10),
##13.2
c(2, 3, 11, 7, 8, 12, 9, 5, 4, 5, 4, 6, 4, 6, 1,
6, 1, 10),
c(1, 4, 9, 7, 8, 12, 10, 2, 5, 2, 5, 3, 5, 3, 6,
3, 6, 11),
c(5, 6, 10, 1, 4, 9, 11, 7, 2, 7, 2, 8, 2, 8, 3,
8, 3, 12),
c(5, 6, 10, 2, 3, 11, 12, 4, 7, 4, 7, 1, 7, 1, 8,
1, 8, 9),
c(5, 6, 10, 7, 8, 12, 2, 11, 1, 11, 1, 9, 11, 9, 3,
9, 3, 4),
c(2, 3, 11, 4, 1, 9, 5, 10, 8, 10, 8, 12, 10, 12, 6,
12, 6, 7),
##13.3
c(7, 8, 12, 13, 3, 11, 13, 11, 6, 13, 6, 5, 13, 5, 9, 13,
9, 4, 13, 4, 1, 13, 1,10,13,10, 2,13,2,3),
c(2, 3, 11, 13, 6, 10, 13, 10, 1, 13, 1, 4, 13, 4, 12, 13,
12, 7, 13, 7, 8, 13, 8, 9,13,9, 5, 13,5,6),
c(7, 8, 12, 13, 6, 5, 13, 5, 9, 13, 9, 4, 13, 4, 3, 13,
3, 11, 13, 11, 2, 13, 2, 1,13,1,10,13,10,6),
c(2, 3, 11, 13, 4, 1, 13, 1, 10, 13, 10, 6, 13, 6, 7,13,
7, 12, 13, 12, 8, 13, 8, 5,13,5, 9,13, 9,4),
c(1, 4, 9, 13, 8, 12, 13, 12, 3, 13, 3, 2, 13, 2, 10, 13,
10, 5, 13, 5, 6, 13, 6,11,13,11,7,13, 7,8),
c(7, 8, 12, 13, 5, 6, 13, 6, 11, 13, 11, 3, 13, 3, 4, 13,
4, 9, 13, 9, 1, 13, 1, 2,13,2,10,13,10,5),
c(1, 4, 9, 13, 3, 2, 13, 2, 10, 13, 10, 5, 13, 5, 8, 13,
8, 12,13, 12, 7, 13, 7, 6,13,6,11,13,11,3),
c(5, 6, 10, 13, 1, 9, 13, 9, 8, 13, 8, 7, 13, 7, 11, 13,
11, 2,13, 2, 3, 13, 3, 12, 13, 12, 4,13, 4, 1),
c(5, 6, 10, 13, 8, 7, 13, 7, 11, 13, 11, 2, 13, 2, 1, 13,
1, 9,13, 9, 4, 13, 4, 3, 13, 3, 12,13, 12, 8),
c(1, 4, 9, 13, 2, 3, 13, 3, 12, 13, 12, 8, 13, 8, 5, 13,
5, 10, 13, 10, 6, 13, 6, 7, 13, 7, 11,13, 11, 2),
c(5, 6, 10, 13, 7, 8, 13, 8, 9, 13, 9, 1, 13, 1, 2, 13,
2, 11, 13, 11, 3, 13, 3, 4, 13, 4, 12,13, 12, 7),
c(2, 3, 11, 13, 1, 4, 13, 4, 12, 13, 12, 7, 13, 7, 6, 13,
6, 10, 13, 10, 5, 13, 5, 8, 13, 8, 9, 13, 9, 1),
##13.4
c(13, 3, 11, 13, 11, 6, 13, 6, 7, 13, 7, 12, 13, 12, 8, 13,
8, 5, 13, 5, 9, 13, 9, 4,
13, 4, 1, 13, 1, 10, 13,10, 2, 13, 2, 3),
c(13, 4, 12, 13, 12, 7, 13, 7, 8, 13, 8, 9, 13, 9, 5, 13,
5, 6, 13, 6, 10, 13, 10, 1, 13,
1, 2, 13, 2, 11, 13, 11, 3, 13, 3, 4),
c(13, 2, 10, 13, 10, 5, 13, 5, 6, 13, 6, 11, 13, 11, 7, 13,
7, 8, 13, 8, 12, 13,12, 3, 13,
3, 4, 13, 4, 9, 13, 9, 1, 13, 1, 2),
c(13, 1, 9, 13, 9, 8, 13, 8, 5, 13, 5, 10, 13, 10, 6, 13,
6, 7, 13, 7, 11, 13,11, 2, 13,
2, 3, 13, 3, 12, 13, 12, 4, 13, 4, 1),
##13.5.1
c(7, 8, 12, 2, 1, 10, 3, 4, 11, 4, 11, 6, 4, 6, 9,
6, 9, 5),
c(3, 2, 11, 8, 5, 9, 4, 1, 12, 1, 12, 7, 1, 7, 10,
7, 10, 6),
c(1, 4, 9, 6, 7, 11, 2, 3, 10, 3, 10, 5, 3, 5, 12,
5, 12, 8),
c(5, 6, 10, 4, 3, 12, 1, 2, 9, 2, 9, 8, 2, 8, 11,
8, 11, 7),
##13.5.2
c(1, 2, 10, 8, 5, 9, 8, 9, 4, 8, 4,12, 4,12, 3,12, 3,11,
12,11, 7,11, 7, 6, 7, 6, 8, 6, 8, 5),
c(8, 5, 9, 3, 4, 12, 3, 12, 7, 3, 7,11, 7,11, 6,11, 6,10,
11,10, 2,10, 2, 1, 2, 1, 3, 1, 3, 4),
c(6, 7, 11, 1, 2, 10, 1, 10, 5, 1, 5, 9, 5, 9, 8, 9, 8,12,
9,12, 4,12, 4, 3, 4, 3, 1, 3, 1, 2),
c(3, 4, 12, 6, 7, 11, 6, 11, 2, 6, 2, 10,2, 10, 1,10,1, 9,
10,9, 5, 5, 9, 8, 5, 8, 6, 8, 6, 7)),
c(1,4,2,1,6,4,1,5,6,3,6,4))
switch23 <- function(x){
num <- length(x) / 3
temp <- x[0: (num-1)*3 + 2]
x[0: (num-1)*3 + 2] <- x[0: (num-1)*3+ 3]
x[0: (num-1)*3 + 3] <- temp
x
}
SwitchSeq <- function(ed){
if (is.list(ed)){
lapply(1:length(ed), function(x) SwitchSeq(ed[[x]]))
}
else
switch23(ed)
}
EdgeSequence2 <- SwitchSeq(EdgeSequence1)
getedge <- function(x){
case <- x[1]
rotation <- x[2]
map <- rep(0,8)
for(i in 1:8){
temp <- as.integer(BasicRotation[rotation,][i])
map[temp] <- i
}
sapply(BasicEdges[[case-1]], function(x){
if (x!=13){
EndP1 <- EdgePoints[x,2]
EndP2 <- EdgePoints[x,3]
newEdge <- EdgePoints[(EdgePoints[,2]==map[EndP1]
&EdgePoints[,3]==map[EndP2])|
(EdgePoints[,3]==map[EndP1]
&EdgePoints[,2]==map[EndP2]),][1]
}
else newEdge <- 13
newEdge})
}
GetEdges <- local({
Edges <- apply(CaseRotationFlip[-c(1,256),], 1, getedge)
Case <- cbind(seq(1:256), CaseRotationFlip[,c(1,3)])
Edges <- apply(Case[-c(1,256),], 1, function(x){
case <- x[2]-1
EdgeNum <- x[1]-1
if (x[3]==1)
sapply(EdgeSequence1[[case]], function(x) Edges[[EdgeNum]][x])
else sapply(EdgeSequence2[[case]], function(x) Edges[[EdgeNum]][x])
})
Edges
})
BasicFace <- list(c(0),c(0),c(1),c(7),c(0),c(2,7),c(1,2,6,7),c(0),c(0),
c(5,6,7),c(0), c(1,4,7),c(1,2,3,4,5,6,7),c(0))
FacePoints <- matrix(c(seq(1,6),1,2,4,1,1,5,6,7,7,8,3,7,2,3,3,4,2,
6,5,6,8,5,4,8),
ncol=5)
FacePoints <- cbind(FacePoints, apply(FacePoints[,2:5],1,prod))
getface <- function(x) {
case <- x[1]
rotation <- x[2]
map <- rep(0,8)
for(i in 1:8){
temp <- as.integer(BasicRotation[rotation,][i])
map[temp] <- i
}
sapply(BasicFace[[case-1]], function(x){
EndP <- rep(0,4)
if (x==0) newFace <- 0
else if (x==7) newFace <- 7
else {
for (i in 1:4){
point <- FacePoints[x,i+1]
EndP[i] <- map[point]
}
newFace<- FacePoints[FacePoints[,6]==prod(EndP[1:4]),][1]
}
newFace})
}
flipface <- function(case, face){
if (face!=0){
index <- explode(case+1)
if (sum(index) > 4)
index <- ifelse(index==0,1,0)
if (face!=7 && index[FacePoints[face,2]]==0)
face <- -face
else if (face==7){
tcase <- CaseRotationFlip[case+1,1]-1
if ((tcase == 4 || tcase==6 || tcase==10 ||tcase==12)
&& !(index[1]+index[7]==2) && !(index[3]+index[5]==2))
face <- -face
else if (tcase==7
&& !(index[1]+index[7]==0) && !(index[3]+index[5]==0))
face <- -face
}
}
face
}
GetFaces <- local({
Faces <- apply(CaseRotationFlip[-c(1,256),], 1, getface)
for (i in 1:254)
for(j in 1:length(Faces[[i]]))
Faces[[i]][j] <- flipface(i, Faces[[i]][j])
Faces
})
special <- list(name = c(3, 4, 6, 7, 10,12,13),
nface= c(1, 1, 2, 4, 3, 3, 7),
sev = c(0, 1, 1, 1, 1, 1, 1),
nedge = c(18,24,48,177,96,96,816),
ind = list(c(0,1), c(0,1), c(0,2,1,3),
c(0,8,4,12,2,10,1,9,6,14,5,13,3,11,15,7),
c(0,4,1,5,2,6,3,7),c(0,4,2,6,1,5,3,7),
c(0,1,2,4,8,16,32,3,9,17,33,6,18,34,12,20,36,24,
40,35,25,22,44,19,41,38,28,83,105,102,92)),
position=list(list(c(1:6),c(7:18)),
list(c(1:6),c(7:24)),
list(c(1:9),c(10:33),c(34:48),c(34:48)),
list(c(1:9),c(1:9),c(10:24),c(10:24),
c(25:39),c(25:39),c(40:54), c(40:54),
c(55:81), c(55:81),c(82:108),c(82:108),
c(109:135),c(109:135),c(136:150),
c(151:177)),
list(c(1:12),c(13:36),c(37:60),c(37:60),
c(61:84), c(61:84),c(85:96),c(85:96)),
list(c(1:12), c(13:36), c(37:60), c(37:60),
c(61:84), c(61:84),c(85:96),c(85:96)),
list(c(1:12),
c(13:30),c(31:48),c(49:66),c(67:84),
c(85:102),c(103:120),
c(121:150),c(151:180),c(181:210),
c(211:240),c(241:270),c(271:300),
c(301:330),
c(331:360),c(361:390),c(391:420),
c(421:450),c(451:480),
c(481:516),c(517:552),c(553:588),
c(589:624),
c(625:642),c(643:660),c(661:678),
c(679:696),
c(697:726),c(727:756),c(757:786),
c(787:816))
)
)
list(Edges = GetEdges, Faces = GetFaces, EdgePoints = EdgePoints,
FacePoints = FacePoints, CaseRotationFlip = CaseRotationFlip,
special = special)
})
Faces <- PreProcessing$Faces
Edges <- PreProcessing$Edges
EdgePoints <- PreProcessing$EdgePoints
FacePoints <- PreProcessing$FacePoints
CaseRotationFlip <- PreProcessing$CaseRotationFlip
special <- PreProcessing$special
fgrid <- function(fun, x, y, z) {
g <- expand.grid(x = x, y = y, z = z)
array(fun(g$x, g$y, g$z), c(length(x), length(y), length(z)))
}
faceType <- function(v, nx, ny, level, maxvol) {
if(level==maxvol)
p <- v >= level
else p <- v > level
v[p] <- 1; v[! p] <- 0
v[-nx, -ny] + 2 * v[-1, -ny] + 4 * v[-1, -1] + 8 * v[-nx, -1]
}
levCells <- function(v, level, maxvol) {
nx <- dim(v)[1]
ny <- dim(v)[2]
nz <- dim(v)[3]
cells <- vector("list", nz - 1)
types <- vector("list", nz - 1)
bottomTypes <- faceType(v[,,1], nx, ny, level, maxvol)
for (k in 1 : (nz - 1)) {
topTypes <- faceType(v[,, k + 1], nx, ny, level, maxvol)
cellTypes <- bottomTypes + 16 * topTypes
contourCells <- which(cellTypes > 0 & cellTypes < 255)
cells[[k]] <- contourCells + (nx - 1) * (ny - 1) * (k - 1)
types[[k]] <- as.integer(cellTypes[contourCells])
bottomTypes <- topTypes
}
cells <- unlist(cells)
i <- as.integer((cells - 1) %% (nx - 1) + 1)
j <- as.integer(((cells - 1) %/% (nx - 1)) %% (ny - 1) + 1)
k <- as.integer((cells - 1) %/% ((nx - 1) * (ny - 1)) + 1)
t <- unlist(types)
list(i = i, j = j, k = k, t = t)
}
CalPoint <- function(x1,x2,y1,y2,z1,z2,v1,v2){
s <- v1 / (v1-v2)
x <- x1+s*(x2-x1)
y <- y1+s*(y2-y1)
z <- z1+s*(z2-z1)
c(x,y,z)
}
GetPoints<-function(edge, p1, info){
##**** need better name than info
## info is the output from GetBasic()
x1 <- EdgePoints[edge,2]
x2 <- EdgePoints[edge,3]
c((1-floor(x1/9))*info[p1+x1-1,1]+floor(x1/9)*info[p1,1],
(1-floor(x1/9))*info[p1+x2-1,1]+floor(x1/9)*info[p1+1,1],
(1-floor(x1/9))*info[p1+x1-1,2]+floor(x1/9)*info[p1+1,2],
(1-floor(x1/9))*info[p1+x2-1,2]+floor(x1/9)*info[p1+2,2],
(1-floor(x1/9))*info[p1+x1-1,3]+floor(x1/9)*info[p1+1,3],
(1-floor(x1/9))*info[p1+x2-1,3]+floor(x1/9)*info[p1+5,3],
(1-floor(x1/9))*info[p1+x1-1,4]+floor(x1/9)*(0*info[p1+1,3]+1),
(1-floor(x1/9))*info[p1+x2-1,4]+floor(x1/9)*(0*info[p1+1,3]-1))
}
FaceNo7 <- function(faces, p1, info){
##**** need better name than info
## info is the output from GetBasic()
index <- ifelse(faces > 0, 1, -1)
faces <- abs(faces)
e1 <- FacePoints[faces,2]
e2 <- FacePoints[faces,3]
e3 <- FacePoints[faces,4]
e4 <- FacePoints[faces,5]
A <- info[p1+e1-1,4]
B <- info[p1+e2-1,4]
C <- info[p1+e3-1,4]
D <- info[p1+e4-1,4]
index <- index*ifelse (A*B-C*D > 0, 1, -1)
ifelse(index==1, 1, 0)
}
Face7 <- function(faces, p1, info){
## info is the output from GetBasic()
index <- ifelse(faces > 0, 1, -1)
A0 <- info[p1,4]; B0 <- info[p1+3,4]
C0 <- info[p1+2,4]; D0 <- info[p1+1,4]
A1 <- info[p1+4,4]; B1 <- info[p1+7,4]
C1 <- info[p1+6,4]; D1 <- info[p1+5,4]
a <- (A1 - A0)*(C1 - C0) - (B1 - B0)*(D1 - D0)
b <- C0*(A1 - A0) + A0*(C1 - C0) - D0*(B1 - B0) - B0*(D1 - D0)
c <- A0*C0 - B0*D0
tmax <- -b/(2*a)
maximum <- a*tmax^2 + b*tmax + c
maximum <- ifelse(maximum=="NaN",-1,maximum)
cond1 <- ifelse (a < 0, 1 ,0)
cond2 <- ifelse (tmax > 0, 1 ,0)
cond3 <- ifelse (tmax < 1, 1, 0)
cond4 <- ifelse (maximum >0, 1, 0)
totalcond <- cond1 * cond2 * cond3 * cond4
index <- index*ifelse(totalcond==1, 1, -1)
ifelse(index==1, 1, 0)
}
## GetBasic()--- The output matrix "information" consists of 4 columns
## and #-of-cubes*8 rows. The first 3 columns tell the
## coordinate(x,y,z) of each vertex of the cube, the 4th gives the
## intensity minus the threshold, which actually makes the threshold
## eaqual to 0. This is convenient for further judgment of subcases
GetBasic <- function(R, vol, level, v) {
cube.1 <- cbind(v$i[R], v$j[R], v$k[R])
index <- matrix(c(0,1,1,0,0,1,1,0,
0,0,1,1,0,0,1,1,
0,0,0,0,1,1,1,1),
nrow=8)
ax.inc <- c(1,1,1)
ver.inc <- t(apply(index,1, function(x) x*ax.inc))
cube.co <-
kronecker(rep(1,nrow(cube.1)),ver.inc) + kronecker(cube.1,rep(1,8))
value <- vol[cube.co] - level
information <- cbind(cube.co, value)
information <- rbind(information, rep(0, 4))
p1 <- (1:length(R) - 1) * 8 + 1
cases <- v$t[R]
list(information=information, p1 = p1, cases=cases)
}
PreRender <- function(edges, p1, type, info) {
if(type==1){
if (typeof(edges)=="list"){
count <- sapply(edges, function(x) length(x))
edges <- cbind(unlist(edges), rep(p1,count))
}
else{
count <- nrow(edges)
edges <- cbind(as.vector(t(edges)), rep(p1,each=count))
}
}
else{
if (is.vector(edges))
edges <- matrix(edges, ncol = length(edges))
p1 <- edges[, 1]
count <- ncol(edges) - 1
edges <- cbind(as.vector(t(edges[, -1])), rep(p1, each = count))
}
##The output of GetPoints() are coordinates of cubes.
info <- GetPoints(edges[,1],edges[,2], info)
info <- matrix(info,ncol=8)
##The output of CalPoint() are coordinates of triangles.
info <- CalPoint(info[,1],info[,2],info[,3],info[,4],
info[,5],info[,6],info[,7],info[,8])
matrix(info,ncol=3)
}
rescale <- function(i, x) {
nx <- length(x)
low <- pmin(pmax(1, floor(i)), nx - 1)
x[low] + (i - low) * (x[low + 1] - x[low])
}
computeContour3d <- function (vol, maxvol, level,
x = 1:dim(vol)[1],
y = 1:dim(vol)[2],
z = 1:dim(vol)[3], mask) {
nx <- length(x)
ny <- length(y)
nz <- length(z)
if (is.function(mask)) mask <- fgrid(mask, x, y, z)
if (! all(mask)) vol[! mask] <- NA
v <- levCells(vol, level, maxvol)
tcase <- CaseRotationFlip[v$t+1,1]-1
R <- which(tcase %in% c(1,2,5,8,9,11,14))
if (length(R) > 0){
Basics <- GetBasic(R, vol, level, v)
information <- Basics$information
p1 <- Basics$p1
cases <- Basics$cases
edges <- Edges[cases]
triangles <- PreRender(edges, p1,type=1, information)
}
else triangles <- matrix(0, nrow=0,ncol=3) # emty contour, e.g.
for (i in 1:length(special$name)){
R <- which(tcase == special$name[i])
if (length(R) > 0) {
Basics <- GetBasic(R, vol, level, v)
information <- Basics$information
p1 <- Basics$p1
cases <- Basics$cases
nface <- special$nface[i]
nedge <- special$nedge[i]
faces <- matrix(unlist(Faces[cases]), ncol = nface, byrow = TRUE)
if (i==1)
index <- FaceNo7(faces[, 1], p1, information)
else if (i==2)
index <- Face7(faces[, 1], p1, information)
else{
index <- Face7(faces[, nface], p1, information)*2^(nface-1)
for(j in 1:(nface-1)){
temp <- FaceNo7(faces[, j], p1, information)
index <- index + temp * 2^(j-1)
}
}
edges <- matrix(unlist(Edges[cases]), ncol = nedge, byrow = TRUE)
edges <- cbind(edges, p1, index)
ind <- special$ind[[i]]
position <- special$position[[i]]
for (j in 1:length(ind)){
ed <- edges[which(index == ind[j]), c(nedge+1, position[[j]])]
if (length(ed) > 0) {
prtri <- PreRender(ed,nedge+1,type=2, information)
triangles <- rbind(triangles, prtri)
}
}
}
}
if (! identical(x, 1 : nx)) triangles[,1] <- rescale(triangles[,1], x)
if (! identical(y, 1 : ny)) triangles[,2] <- rescale(triangles[,2], y)
if (! identical(z, 1 : nz)) triangles[,3] <- rescale(triangles[,3], z)
triangles
}
################################################################################
contourTriangles <- function(vol, maxvol, level,
x = 1:dim(vol)[1],
y = 1:dim(vol)[2],
z = 1:dim(vol)[3],
mask = NULL, color = "white", color2 = NA,
alpha = 1, fill = TRUE,
col.mesh = if (fill) NA else color,
material = "default", smooth = 0) {
if (length(level) > 1) {
val <- vector("list", length(level))
for (i in seq(along = level)) {
m <- if (is.list(mask)) mask[[i]] else mask
col <- if (length(color) > 1) color[[i]] else color
col2 <- if (length(color2) > 1) color2[[i]] else color2
a <- if (length(alpha) > 1) alpha[[i]] else alpha
fl <- if (length(fill) > 1) fill[[i]] else fill
cm <- if (length(col.mesh) > 1) col.mesh[[i]] else col.mesh
mat <- if (length(material) > 1) material[[1]] else material
sm <- if (length(smooth) > 1) smooth[[1]] else smooth
val[[i]] <- contourTriangles(vol, maxvol, level[i], x, y, z, m,
col, col2, a, fl, cm, mat, sm)
}
val
}
else makeTriangles(computeContour3d(vol, maxvol, level, x, y, z, mask),
color = color, color2 = color2, alpha = alpha,
fill = fill, col.mesh = col.mesh,
material = material, smooth = smooth)
}
################################################################################
getlevel <- function(f, x, y, z) {
if (! all(is.finite(x), is.finite(y), is.finite(z)))
stop("'x', 'y', and 'z' values must be finite and non-missing")
if (is.function(f) || is.array(f) && length(dim(f)) == 3){
if (is.function(f)){
if (length(formals(f)) < 3)
stop("The function must have at least 3 arguments.")
vol <- fgrid(f, x, y, z)
}
else{
if (dim(f)[1] != length(x) || dim(f)[2] != length(y) || dim(f)[3] != length(z))
stop("dimensions of f do not match x, y, or z")
vol <- f
}
maxvol <- max(vol)
minvol <- min(vol)
print(paste("Your level is betwen ", minvol, " and ", maxvol, sep=""),quote=F)
}
}
################################################################################
makeTriangles <- function(v1, v2, v3,
color = "red", color2 = NA, alpha = 1,
fill = TRUE, col.mesh = if (fill) NA else color,
smooth = 0, material = "default") {
if (missing(v2) || missing(v3)) {
if (missing(v2) && missing(v3))
v <- unzipTriangleMatrix(v1)
else if (missing(v3))
v <- ve2t(list(vb = v1, ib = v2))
else stop("unknown form of triangle specification")
v1 <- v$v1
v2 <- v$v2
v3 <- v$v3
}
tris <- structure(list(v1 = v1, v2 = v2, v3 = v3,
color = color, color2 = color2, fill = fill,
material = material, col.mesh = col.mesh,
alpha = alpha, smooth = smooth),
class = "Triangles3D")
colorTriangles(tris)
}
is.Triangles3D <- function(x) identical(class(x), "Triangles3D")
updateTriangles <- function(triangles, color, color2, alpha, fill, col.mesh,
material, smooth) {
if (! missing(color)) triangles$color <- color
if (! missing(color2)) triangles$color2 <- color2
if (! missing(fill)) triangles$fill <- fill
if (! missing(col.mesh)) triangles$col.mesh <- col.mesh
if (! missing(material)) triangles$material <- material
if (! missing(alpha)) triangles$alpha <- alpha
if (! missing(smooth)) triangles$smooth <- smooth
colorTriangles(triangles)
}
#**** This assumes comparable scaling of dimensions
#**** 5 is the largest exponent for S that will work; smaller is OK
t2ve <- function (triangles)
{
vb <- rbind(triangles$v1, triangles$v2, triangles$v3)
vbmin <- min(vb)
vbmax <- max(vb)
S <- 10^5
score <- function(v, d) floor(as.vector(v %*% d))
scale <- function(v) (1 - 1 / S) * (v - vbmin) / (vbmax - vbmin)
d <- c(S, S^2, S^3)
scores <- score(scale(vb), d)
vb <- vb[! duplicated(scores),]
scores <- score(scale(vb), d)
ib <- rbind(match(score(scale(triangles$v1), d), scores),
match(score(scale(triangles$v2), d), scores),
match(score(scale(triangles$v3), d), scores))
list(vb = t(vb), ib = ib)
}
ve2t <- function(ve) {
list (v1 = t(ve$vb[,ve$ib[1,]]),
v2 = t(ve$vb[,ve$ib[2,]]),
v3 = t(ve$vb[,ve$ib[3,]]))
}
unzipTriangleMatrix <- function(tris) {
if (ncol(tris) != 3)
stop("triangle matrix must have three columns.")
if (nrow(tris) %% 3 != 0)
stop("number of rows in triangle matrix must be divisible by 3")
n <- nrow(tris) / 3
list(v1 = tris[3 * (1 : n) - 2,],
v2 = tris[3 * (1 : n) - 1,],
v3 = tris[3 * (1 : n),])
}
zipTriangles <- function(tris) {
n <- nrow(tris$v1)
if(is.null(n)) n <- 1
else if (nrow(tris$v2) != n || nrow(tris$v3) != n)
stop("vertex arrays must have the same number of rows")
v <- matrix(0, nrow = 3 * n, ncol = 3)
v[3 * (1 : n) - 2,] <- tris$v1
v[3 * (1 : n) - 1,] <- tris$v2
v[3 * (1 : n),] <- tris$v3
v
}
colorTriangles <- function(triangles) {
if (is.function(triangles$color) || is.function(triangles$color2)) {
v <- (triangles$v1 + triangles$v2 + triangles$v3) / 3
if (is.function(triangles$color))
triangles$color <- triangles$color(v[,1], v[,2], v[,3])
if (is.function(triangles$color2))
triangles$color2 <- triangles$color2(v[,1], v[,2], v[,3])
if (is.function(triangles$col.mesh))
triangles$col.mesh <- triangles$col.mesh(v[,1], v[,2], v[,3])
}
triangles
}
colorScene <- function(scene) {
if (is.Triangles3D(scene))
colorTriangles(scene)
else lapply(scene, colorTriangles)
}
## **** better to make new triangles including only requested components?
canonicalizeAndMergeScene <- function(scene, ...) {
which <- list(...)
if (is.Triangles3D(scene)) {
#cat("scene v1 is ", print(scene$v1), "\n")
n.tri <- nrow(scene$v1)
if(is.null(n.tri)) n.tri<-1
for (n in which)
if (length(scene[[n]]) != n.tri)
scene[[n]] <- rep(scene[[n]], length = n.tri)
scene
}
else {
scene <- lapply(scene, canonicalizeAndMergeScene, ...)
x <- scene[[1]]
x$v1 <- do.call(rbind, lapply(scene, function(x) x$v1))
x$v2 <- do.call(rbind, lapply(scene, function(x) x$v2))
x$v3 <- do.call(rbind, lapply(scene, function(x) x$v3))
for (n in which)
x[[n]] <- do.call(c, lapply(scene, function(x) x[[n]]))
x
}
}
expandTriangleGrid <- function(x, y) {
nx <- length(x) - 1
ny <- length(y) - 1
A <- c(0, 0)
B <- c(1, 0)
C <- c(1, 1)
D <- c(0, 1)
g <- expand.grid(x = 1 : nx, y = 1 : ny)
even <- (g$x + g$y) %% 2 == 0
gx11 <- ifelse(even, g$x + A[1], g$x + A[1])
gy11 <- ifelse(even, g$y + A[2], g$y + A[2])
gx12 <- ifelse(even, g$x + A[1], g$x + B[1])
gy12 <- ifelse(even, g$y + A[2], g$y + B[2])
i1 <- rbind(cbind(gx11, gy11), cbind(gx12, gy12))
gx21 <- ifelse(even, g$x + B[1], g$x + B[1])
gy21 <- ifelse(even, g$y + B[2], g$y + B[2])
gx22 <- ifelse(even, g$x + C[1], g$x + C[1])
gy22 <- ifelse(even, g$y + C[2], g$y + C[2])
i2 <- rbind(cbind(gx21, gy21), cbind(gx22, gy22))
gx31 <- ifelse(even, g$x + C[1], g$x + D[1])
gy31 <- ifelse(even, g$y + C[2], g$y + D[2])
gx32 <- ifelse(even, g$x + D[1], g$x + D[1])
gy32 <- ifelse(even, g$y + D[2], g$y + D[2])
i3 <- rbind(cbind(gx31, gy31), cbind(gx32, gy32))
v1 <- cbind(x[i1[,1]], y[i1[,2]])
v2 <- cbind(x[i2[,1]], y[i2[,2]])
v3 <- cbind(x[i3[,1]], y[i3[,2]])
list(v1 = v1, v2 = v2, v3 = v3)
}
## adapted from lattice ltransform3dto3d
trans3dto3d <- function (x, R.mat) {
if (length(x) == 0)
return(x)
val <- R.mat %*% rbind(t(x), 1)
val[1, ] <- val[1, ]/val[4, ]
val[2, ] <- val[2, ]/val[4, ]
val[3, ] <- val[3, ]/val[4, ]
t(val[1:3, , drop = FALSE])
}
transformTriangles <- function(triangles, R) {
tr <- function(v) trans3dto3d(v, R)
triangles$v1 <- tr(triangles$v1)
triangles$v2 <- tr(triangles$v2)
triangles$v3 <- tr(triangles$v3)
triangles
}
transformScene <- function(scene, rot.mat) {
if (is.Triangles3D(scene))
transformTriangles(scene, rot.mat)
else lapply(scene, transformTriangles, rot.mat)
}
translateTriangles <- function(triangles, x = 0, y = 0, z = 0) {
M <- diag(4)
M[1:3,4] <- c(x, y, z)
transformTriangles(triangles, M)
}
scaleTriangles <- function(triangles, x = 1, y = x, z = x) {
M <- diag(c(x, y, z, 1))
transformTriangles(triangles, M)
}
## triangleNormals computes the normal vectors to a collection of
## triangles as the vector crossprocuct of the direction from v1 to v2
## and the direction from v2 to v3. The result is an n by 3 matrix of
## unit representing the n unit normal vectors.
triangleNormals <- function(triangles) {
x <- triangles$v2 - triangles$v1
y <- triangles$v3 - triangles$v2
z <- cbind(x[,2]*y[,3] - x[,3]*y[,2],
x[,3]*y[,1] - x[,1]*y[,3],
x[,1]*y[,2] - x[,2]*y[,1])
z / sqrt(rowSums(z^2))
}
# adapted from lattice ltransform3dMatrix
trans3dMat <- function (screen, P = diag(4)) {
givens4 <- function(i, j, gamma) {
T <- diag(4)
cgamma <- cos(gamma)
sgamma <- sin(gamma)
T[c(i,j),c(i,j)] <- matrix(c(cgamma, sgamma, -sgamma, cgamma), 2, 2)
T
}
screen.names <- names(screen)
for (i in seq(along = screen.names)) {
if (screen.names[i] == "x")
P <- givens4(2, 3, screen[[i]] * pi/180) %*% P
else if (screen.names[i] == "y")
P <- givens4(1, 3, -screen[[i]] * pi/180) %*% P #**** whi negative?
else if (screen.names[i] == "z")
P <- givens4(1, 2, screen[[i]] * pi/180) %*% P
}
P
}
makeViewTransform <- function(ranges, scale, aspect, screen, R.mat) {
m <- c(mean(ranges$xlim), mean(ranges$ylim), mean(ranges$zlim))
s <- 0.5 * c(diff(ranges$xlim), diff(ranges$ylim), diff(ranges$zlim))
if (! scale) s <- rep(max(s), 3)
else s <- s / c(1, aspect)
A <- diag(1 / c(s, 1))
A[1:3, 4] <- -m / s
trans3dMat(screen, R.mat %*% A)
}
trianglesRanges <- function(triangles, xlim, ylim, zlim) {
v1 <- triangles$v1
v2 <- triangles$v2
v3 <- triangles$v3
if (is.null(xlim)) xlim <- range(v1[,1], v2[,1], v3[,1], na.rm = TRUE)
if (is.null(ylim)) ylim <- range(v1[,2], v2[,2], v3[,2], na.rm = TRUE)
if (is.null(zlim)) zlim <- range(v1[,3], v2[,3], v3[,3], na.rm = TRUE)
list(xlim = xlim, ylim = ylim, zlim = zlim)
}
sceneRanges <- function(scene, xlim, ylim, zlim) {
if (is.Triangles3D(scene))
trianglesRanges(scene, xlim, ylim, zlim)
else {
ranges <- lapply(scene, trianglesRanges, xlim, ylim, zlim)
list(xlim = range(sapply(ranges,function(x) x$xlim)),
ylim = range(sapply(ranges,function(x) x$ylim)),
zlim = range(sapply(ranges,function(x) x$zlim)))
}
}
addTrianglesPerspective <- function(triangles, distance) {
pt <- function(v) {
v[, 1] <- v[, 1] / (1 / distance - v[, 3])
v[, 2] <- v[, 2] / (1 / distance - v[, 3])
v
}
triangles$v1 <- pt(triangles$v1)
triangles$v2 <- pt(triangles$v2)
triangles$v3 <- pt(triangles$v3)
triangles
}
addPerspective <- function(scene, distance) {
if (is.Triangles3D(scene))
addTrianglesPerspective(scene, distance)
else lapply(scene, addTrianglesPerspective, distance)
}
screenRange <- function(v1, v2, v3)
range(v1[,1:2], v2[,1:2], v3[,1:2], na.rm = TRUE)
vertexTriangles <- function(ve) {
n.vert <- ncol(ve$vb)
ib <- ve$ib
vt <- function(i) which(ib[1,] == i | ib[2,] == i | ib[3,] == i)
lapply(1 : n.vert, vt)
}
# faster version
vertexTriangles <- function(ve) {
n.vert <- ncol(ve$vb)
val <- vector("list", n.vert)
ib <- ve$ib
for (i in 1 : ncol(ib)) {
val[[ib[1,i]]] <- c(val[[ib[1,i]]], i)
val[[ib[2,i]]] <- c(val[[ib[2,i]]], i)
val[[ib[3,i]]] <- c(val[[ib[3,i]]], i)
}
val
}
vertexNormals <- function(vt, N) {
vn <- function(tris) {
z <- apply(N[tris,,drop = FALSE], 2, mean, na.rm = TRUE);
z <- z / sqrt(sum(z^2))
if (any(is.na(z))) c(1,0,0) else z
}
t(sapply(vt, vn))
}
# faster version
vertexNormals <- function(vt, N) {
val <- matrix(0, nrow = length(vt), ncol = 3)
for (i in seq(along = vt)) {
Ni <- N[vt[[i]],,drop = FALSE]
Ni1 <- Ni[,1]
Ni2 <- Ni[,2]
Ni3 <- Ni[,3]
z1 <- if (any(is.na(Ni1))) mean(Ni1, na.rm = TRUE)
else sum(Ni1) / length(Ni1)
z2 <- if (any(is.na(Ni2))) mean(Ni2, na.rm = TRUE)
else sum(Ni2) / length(Ni2)
z3 <- if (any(is.na(Ni3))) mean(Ni3, na.rm = TRUE)
else sum(Ni3) / length(Ni3)
z <- c(z1, z2, z3)
z <- z / sqrt(sum(z^2))
val[i,] <- if (any(is.na(z))) c(1,0,0) else z
}
val
}
interpolateVertexNormals <- function(VN, ib) {
z <- (VN[ib[1,],] + VN[ib[2,],] + VN[ib[3,],]) / 3
z / sqrt(rowSums(z^2))
}
## triangleVertexNormals computes the normals at the vertices by
## averaging the normals of the incident triangles. This is used by
## the rgl engine. The result form is chosen so zipTriangles can be
## used on it.
triangleVertexNormals <- function(v) {
N <- triangleNormals(v)
ve <- t2ve(v)
vt <- vertexTriangles(ve)
VN <- vertexNormals(vt, N)
list(v1 = VN[ve$ib[1,],], v2 = VN[ve$ib[2,],], v3 = VN[ve$ib[3,],])
}
vertexColors <- function(vt, col) {
C <- t(col2rgb(col))
val <- matrix(0, nrow = length(vt), ncol = 3)
for (i in seq(along = vt)) {
vti <- vt[[i]]
nti <- length(vti)
Ci <- C[vti,,drop = FALSE]
Ci1 <- Ci[,1]
Ci2 <- Ci[,2]
Ci3 <- Ci[,3]
val[i,] <- c(sum(Ci1), sum(Ci2), sum(Ci3)) / nti
}
val
}
interpolateVertexColors <- function(VC, ib) {
TC <- (VC[ib[1,],] + VC[ib[2,],] + VC[ib[3,],]) / 3
rgb(TC[,1], TC[,2], TC[,3], maxColorValue=255)
}
triangleEdges <- function(vb, ib) {
edges <- cbind(ib[c(1,2),], ib[c(2,3),], ib[c(3,1),])
swap <- edges[1,] > edges[2,]
edges[,swap] <- edges[2:1,swap]
edges[,! duplicated(edges, MARGIN = 2)]
}
# faster version
triangleEdges <- function(vb, ib) {
n.vert <- ncol(vb)
edges <- cbind(ib[c(1,2),], ib[c(2,3),], ib[c(3,1),])
swap <- edges[1,] > edges[2,]
edges[,swap] <- edges[2:1,swap]
score <- as.vector(c(1 + n.vert, 1) %*% edges)
edges[,! duplicated(score)]
}
triangleMidTriangles <- function(vb, ib, VN) {
n.vert <- ncol(vb)
edges <- triangleEdges(vb, ib)
vb <- (vb[,edges[1,]] + vb[,edges[2,]]) / 2
d <- c(1 + n.vert, 1)
scores <- as.vector(d %*% edges)
mpi <- function(a, b) {
s <- d[1] * pmin(a, b) + d[2] * pmax(a, b)
match(s, scores)
}
mpi1 <- mpi(ib[1,], ib[2,])
mpi2 <- mpi(ib[2,], ib[3,])
mpi3 <- mpi(ib[3,], ib[1,])
ib <- rbind(mpi1, mpi2, mpi3)
z <- VN[edges[1,],] + VN[edges[2,],]
z <- z / sqrt(rowSums(z^2))
list(vb = vb, ib = ib, VN = z)
}
## surfaceTriangles creates a set of triangles for a grid specified by x,
## y and function falues computed with f if f is a function or taken
## from f if f is a matrix.
surfaceTriangles <- function(x, y, f,
color = "red", color2 = NA, alpha = 1,
fill = TRUE, col.mesh = if (fill) NA else color,
smooth = 0, material = "default") {
if (is.function(f))
ff <- function(ix, iy) f(x[ix], y[iy])
else
ff <- function(ix, iy) f[ix + length(x) * (iy - 1)]
i <- expandTriangleGrid(1 : length(x), 1 : length(y))
i1 <- i$v1
i2 <- i$v2
i3 <- i$v3
v1 <- cbind(x[i1[,1]], y[i1[,2]], ff(i1[,1], i1[,2]))
v2 <- cbind(x[i2[,1]], y[i2[,2]], ff(i2[,1], i2[,2]))
v3 <- cbind(x[i3[,1]], y[i3[,2]], ff(i3[,1], i3[,2]))
na1 <- is.na(v1[,1]) | is.na(v1[,2]) | is.na(v1[,3])
na2 <- is.na(v2[,1]) | is.na(v2[,2]) | is.na(v2[,3])
na3 <- is.na(v3[,1]) | is.na(v3[,2]) | is.na(v3[,3])
nna <- ! (na1 | na2 | na3)
makeTriangles(v1[nna,], v2[nna,], v3[nna,],
color = color, color2 = color2, fill = fill, smooth = smooth,
material = material, col.mesh = col.mesh, alpha = alpha)
}
## pointsTetrahedra computes a collection of tetrahedra centered at
## the specified point locations. This is useful, for example, for
## displaying raw data along with a density contour in a scene
## rendered with standard or grid graphics. Random orientation might
## be useful to avoid strange results at certain lighting angles.
pointsTetrahedra <- function(x, y, z, size = 0.01, color = "black", ...) {
n <- length(x)
if (length(y) != n || length(z) != n)
stop("coordinate vectors must be the same length.")
## Create a basic tetrahedron centered at the origin
a <- sqrt(3) / 2
b <- 1 / (2 * sqrt(3))
h <- sqrt(2 / 3)
mx <- 1 / 2
my <- (a + b) / 4
mz <- h / 4
A <- c( -mx, -my, -mz)
B <- c( 1 - mx, -my, -mz)
C <- c(1 / 2 - mx, a - my, -mz)
D <- c(1 / 2 - mx, b - my, h - mz)
v1 <- rbind(B, A, B, C)
v2 <- rbind(A, B, C, A)
v3 <- rbind(C, D, D, D)
## Scale the tetrahedron
if (length(size) < 3) size <- rep(size, len = 3)
if (n == 1) s <- diag(size)
else s <- diag(size * c(diff(range(x)), diff(range(y)), diff(range(z))))
sv1 <- v1 %*% s
sv2 <- v2 %*% s
sv3 <- v3 %*% s
## Compute the tetrahedra for the points, taking advantage of recycling
x4 <- rep(x, each = 4)
y4 <- rep(y, each = 4)
z4 <- rep(z, each = 4)
V1 <- cbind(x4 + sv1[,1], y4 + sv1[,2], z4 + sv1[,3])
V2 <- cbind(x4 + sv2[,1], y4 + sv2[,2], z4 + sv2[,3])
V3 <- cbind(x4 + sv3[,1], y4 + sv3[,2], z4 + sv3[,3])
makeTriangles(V1, V2, V3, color = color, ...)
}
## Compute for each triangle the indices of triangles that share an
## edge with it. This could be done more efficiently.
triangleNeighbors <- function(tris) {
ve <- t2ve(tris)
vt <- vertexTriangles(ve)
ib <- ve$ib
n.tri <- ncol(ib)
tn <- vector("list", n.tri)
for (i in 1 : n.tri) {
v1 <- unique(vt[[ib[1, i]]])
v2 <- unique(vt[[ib[2, i]]])
v3 <- unique(vt[[ib[3, i]]])
i12 <- intersect(v1, v2)
i23 <- intersect(v2, v3)
i31 <- intersect(v3, v1)
u <- union(union(i12, i23), i31)
tn[[i]] <- u[u != i]
}
tn
}
## 'unique' in unique(vt[[ib[1, i]]]) seems to be unnecessary
## unless a triangle has essentially two vertices or one vertex
triangleNeighbors <- function(tris) {
ve <- t2ve(tris)
vt <- vertexTriangles(ve)
ib <- ve$ib
n.tri <- ncol(ib)
tn <- vector("list", n.tri)
for (i in 1 : n.tri) {
v1 <- vt[[ib[1, i]]]
v2 <- vt[[ib[2, i]]]
v3 <- vt[[ib[3, i]]]
i12 <- intersect(v1, v2)
i23 <- intersect(v2, v3)
i31 <- intersect(v3, v1)
u <- union(union(i12, i23), i31)
tn[[i]] <- u[u != i]
}
tn
}
## Dijkstra's version of Rem's algorithm for computing equivalence
## classes based on a number of vertices 1:nvert and a set of N edges
## provided as an N x 2 matrix.
GetPatches <- function(nvert, edges) {
f <- 1:nvert
if (!(is.vector(edges)) && dim(edges)[1] != 0){
nedge <- nrow(edges)
for (e in 1:nedge) {
p0 <- edges[e, 1]
q0 <- edges[e, 2]
p1 <- f[p0]
q1 <- f[q0]
while (p1 != q1) {
if (q1 < p1) {
f[p0] <- q1
p0 <- p1
p1 <- f[p1]
}
else {
f[q0] <- p1
q0 <- q1
q1 <- f[q1]
}
}
}
}
if(is.vector(edges)){
if(edges[1] < edges[2])
f[edges[2]] <- edges[1]
else f[edges[1]] <- edges[2]
}
for (v in 1:nvert)
f[v] <- f[f[v]]
split(1:nvert,f)
}
## compute the edges to indicate which triangles share an edge -- this
## needs more error checking
triangleNeighborEdges <- function(tn) {
edges <- function(i) {
v <- tn[[i]]
if (length(v) > 0) cbind(i,v)
else numeric(0)
}
do.call(rbind, lapply(1:length(tn), edges))
}
## separate triangles into disconnected chunks
separateTriangles <- function(contour3dObj){
tn <- triangleNeighbors(contour3dObj)
edges <- triangleNeighborEdges(tn)
edges <- edges[edges[,1] < edges[,2],]
p <- GetPatches(length(tn), edges)
newContour3dObj <- vector("list", length(p))
for(i in 1:length(newContour3dObj)){
newContour3dObj[[i]] <- contour3dObj
newContour3dObj[[i]]$v1 <- contour3dObj$v1[p[[i]],]
newContour3dObj[[i]]$v2 <- contour3dObj$v2[p[[i]],]
newContour3dObj[[i]]$v3 <- contour3dObj$v3[p[[i]],]
}
newContour3dObj
}
drawScene.rgl <- function(scene, add = FALSE, ...) {
loadRGL()
if (! rgl.cur())
open3d()
if (!add)
clear3d()
scene <- colorScene(scene)
triangles <- canonicalizeAndMergeScene(scene, "color", "color2", "alpha",
"col.mesh", "fill", "smooth")
use.col2 <- ! all(is.na(triangles$color2))
if (use.col2 && any(is.na(triangles$color2)))
warning(paste("mixing surfaces with and without color2 may not",
"work properly in the rgl engine"))
col <- rep(triangles$color, each = 3)
if (use.col2)
col2 <- rep(triangles$color2, each=3)
alpha <- rep(triangles$alpha, each = 3)
fill <- rep(triangles$fill, each = 3)
col.mesh <- rep(triangles$col.mesh, each = 3)
data <- zipTriangles(triangles)
if (all(fill)) {
front <- "filled"
back <- "filled"
}
else if (any(fill))
##**** handle these by splitting; OK if no alpha < 1
stop(paste("for now rgl engine cannot handle mixed fill/wire",
"frame surfaces"))
else {
front <- "lines"
back <- "lines"
col <- col.mesh
}
oldstyle = getOption("old.misc3d.orientation")
if (! is.null(oldstyle) && oldstyle) {
data <- data[,c(1, 3, 2)]
data[,3] <- -data[,3]
}
if (any(triangles$smooth > 0)) {
if (any(triangles$smooth == 0))
stop(paste("for now for the rgl engine cannot handle mixed",
"smooth/non-smooth surfaces"))
normals <- zipTriangles(triangleVertexNormals(triangles))
}
else normals <- NULL
if (nrow(data) > 0) # to avoid a segfault in rgl
{
if (! use.col2)
triangles3d(data[,1], data[,2], data[,3],
col = col, alpha = alpha, normals = normals,
front = front, back = back, ...)
else {
triangles3d(data[,1], data[,2], data[,3],
col = col, alpha = alpha, normals = normals,
front = front, back = "cull", ...)
triangles3d(data[,1], data[,2], data[,3],
col = col2, alpha = alpha, normals = normals,
front = "cull", back = back, ...)
}
}
}
renderScene <- function(scene, box, fill, col.mesh, add, engine, polynum,
col.bg, depth, newpage) {
triangles <- canonicalizeAndMergeScene(scene, "color", "col.light",
"col.mesh", "fill")
v1 <- triangles$v1
v2 <- triangles$v2
v3 <- triangles$v3
n.tri <- nrow(v1)
if (fill) {
fill <- rep(triangles$fill, length = n.tri)
col.mesh <- rep(triangles$col.mesh, length = n.tri)
}
else col.mesh <- rep(col.mesh, length = n.tri)
col.fill <- ifelse(fill, triangles$col.light, NA)
z <- (v1[,3] + v2[,3] + v3[,3]) / 3
if (depth > 0) {
rgbcol <- col2rgb(col.fill, alpha = TRUE)/255
rgbcol.bg <- col2rgb(col.bg, alpha = TRUE)/255
s <- (1 + depth * z) / (1 + depth * max(z))
col.fill <- rgb(rgbcol[1,] * s + rgbcol.bg[1,] * (1 - s),
rgbcol[2,] * s + rgbcol.bg[2,] * (1 - s),
rgbcol[3,] * s + rgbcol.bg[3,] * (1 - s),
rgbcol[4,])
}
col.mesh <- ifelse(is.na(col.mesh), col.fill, col.mesh)
i <- order(z, na.last = NA)
if (engine == "grid")
render.grid(v1[i,], v2[i,], v3[i,], box, fill[i], col.fill[i],
col.mesh[i], add, polynum, newpage)
else render.standard(v1[i,], v2[i,], v3[i,], box, fill[i], col.fill[i],
col.mesh[i], add)
}
render.standard <- function(v1, v2, v3, box, fill, col.fill, col.mesh, add) {
if (! add) {
# rr <- screenRange(v1, v2, v3)
rr <- range(box)
plot(rr, rr,type="n", axes = FALSE, ann = FALSE)
}
xx <- as.vector(rbind(v1[,1], v2[,1], v3[,1], NA))
yy <- as.vector(rbind(v1[,2], v2[,2], v3[,2], NA))
polygon(xx, yy, col=col.fill, border=col.mesh)
}
render.grid <- function(v1, v2, v3, box, fill, col.fill, col.mesh,
add, polynum, newpage) {
if (! add) {
if (newpage) grid::grid.newpage()
rr <- range(box)
grid::pushViewport(grid::viewport(w = 0.8, h = 0.8,
xscale = rr, yscale = rr,
name = "misc3dScene"))
on.exit(grid::upViewport())
}
xx <- as.vector(rbind(v1[,1], v2[,1], v3[,1]))
yy <- as.vector(rbind(v1[,2], v2[,2], v3[,2]))
n.tri <- nrow(v1)
idlen <- rep(3, n.tri)
start <- 1
end <- start + polynum - 1
while (start <= n.tri) {
end <- min(end, n.tri)
j <- start : end
j3 <- (3 * start - 2) : (3 * end)
gp <- grid::gpar(fill = col.fill[j], col = col.mesh[j])
grid::grid.polygon(x = xx[j3], y = yy[j3], default.units = "native",
gp = gp, id.lengths = idlen[j])
start <- start + polynum
end <- start + polynum
}
}
makePerspMatrix <- function(d) {
rbind(c(1, 0, 0, 0),
c(0, 1, 0, 0),
c(0, 0, 1, 0),
c(0, 0, -1, 1 / d))
}
drawScene <- function(scene, light = c(0, 0, 1),
screen = list(z = 40, x = -60),
scale = TRUE,
R.mat = diag(4),
perspective = FALSE,
distance = if (perspective) 0.2 else 0,
fill = TRUE,
xlim = NULL, ylim = NULL, zlim = NULL,
aspect = c(1, 1),
col.mesh = if (fill) NA else "black",
polynum = 100,
lighting = phongLighting,
add = FALSE,
engine = "standard",
col.bg = "transparent", depth = 0, newpage = TRUE) {
scene <- colorScene(scene)
sr <- sceneRanges(scene, xlim, ylim, zlim)
if (add)
rot.mat <- R.mat
else
rot.mat <- makeViewTransform(sr, scale, aspect, screen, R.mat)
scene <- transformScene(scene, rot.mat)
scene <- lightScene(scene, lighting, light)
if (distance > 0) {
scene <- addPerspective(scene, distance)
rot.mat <- makePerspMatrix(distance) %*% rot.mat
}
box <- as.matrix(expand.grid(sr$xlim, sr$ylim, sr$zlim))
box <- trans3dto3d(box, rot.mat)
renderScene(scene, box, fill, col.mesh, add, engine, polynum,
col.bg, depth, newpage)
invisible(t(rot.mat))
}
loadRGL <- function() {
#if (! suppressWarnings(require(rgl,quietly=TRUE)))
# stop("rgl is not available")
}
phongLighting <- function(normals, view, light, color, color2, alpha,
material = "default") {
if (length(light) == 4) {
LI <- light[4]
light <- light[1:3]
}
else LI <- 1
if (is.character(material))
material <- getMaterial(material)
ambient <- material$ambient
diffuse <- material$diffuse
specular <- material$specular
exponent <- material$exponent
sr <- material$sr
V <- view / sqrt(sum(view^2))
L <- light / sqrt(sum(light^2))
H <- (L + V) / sqrt(sum((L + V)^2))
sgn <- as.vector(normals %*% V) > 0
N <- ifelse(sgn,1, -1) * normals
Is <- as.vector(specular * abs(N %*% H) ^ exponent)
Id <- as.vector(diffuse * pmax(N %*% L,0))
rgbcol <- t(col2rgb(ifelse(sgn, color, color2)) / 255)
Lrgbcol <- pmin(LI * ((ambient + Id + sr * Is) * rgbcol + (1 - sr) * Is),
1)
Lrgbcol[is.na(Lrgbcol)] <- 0
rgb(Lrgbcol[,1], Lrgbcol[,2], Lrgbcol[,3], alpha)
}
materials.database <- new.env(hash = TRUE)
registerMaterial <- function(name, ambient = 0.3, diffuse = 0.7, specular = 0.1, exponent = 10, sr = 0) {
value <- list(ambient = ambient, diffuse = diffuse, specular = specular, exponent = exponent, sr = sr)
assign(name, value, materials.database)
}
registerMaterial("shiny", ambient = 0.3, diffuse = 0.6, specular = 0.9, exponent = 20, sr = 0)
registerMaterial("dull", ambient = 0.3, diffuse = 0.8, specular = 0.0, exponent = 10, sr = 0)
registerMaterial("metal", ambient = 0.3, diffuse = 0.3, specular = 1.0, exponent = 25, sr = 0.5)
registerMaterial("default", ambient = 0.3, diffuse = 0.7, specular = 0.1, exponent = 10, sr = 0)
registerMaterial("metal", ambient = 0.45, diffuse = 0.45, specular = 1.5, exponent = 25, sr = 0.5)
registerMaterial("shiny", ambient = 0.36, diffuse = 0.72, specular = 1.08, exponent = 20, sr = 0)
getMaterial <- function(name) {
if (exists(name, materials.database, inherits = FALSE))
get(name, materials.database)
else get("default", materials.database, inherits = FALSE)
}
perspLighting <- function(normals, view, light, color, color2, alpha,
material = "default") {
if (length(light) == 4) {
LI <- light[4]
light <- light[1:3]
}
else LI <- 1
if (is.character(material))
material <- getMaterial(material)
exponent <- material$exponent
V <- view / sqrt(sum(view^2))
L <- light / sqrt(sum(light^2))
sgn <- as.vector(normals %*% V) > 0
N <- ifelse(sgn,1, -1) * normals
I <- (pmax(1 + as.vector(N %*% L), 0) / 2) ^ (exponent * (3 / 40))
Lrgbcol <- I * LI * t(col2rgb(ifelse(sgn, color, color2)) / 255)
rgb(Lrgbcol[,1], Lrgbcol[,2], Lrgbcol[,3], alpha)
}
triangleNormalsPhong <- function(triangles) {
N <- triangleNormals(triangles)
ve <- t2ve(triangles)
vt <- vertexTriangles(ve)
VN <- vertexNormals(vt, N)
interpolateVertexNormals(VN, ve$ib)
}
triangleNormalsPhongEX <- function(triangles, reps = 1) {
N <- triangleNormals(triangles)
ve <- t2ve(triangles)
vt <- vertexTriangles(ve)
VN <- vertexNormals(vt, N)
vb <- ve$vb
ib <- ve$ib
n.tri <- nrow(N)
while (reps > 0) {
reps <- reps - 1
n.ver <- nrow(VN)
mt <- triangleMidTriangles(vb, ib, VN)
vb <- cbind(vb, mt$vb)
VN <- rbind(VN, mt$VN)
mtib <- mt$ib + n.ver
ib <- matrix(rbind(ib[1,], mtib[1,],mtib[3,],
mtib[1,], ib[2,],mtib[2,],
mtib[2,], ib[3,], mtib[3,],
mtib), nrow = 3)
for (i in seq(along = triangles))
if (length(triangles[[i]]) == n.tri)
triangles[[i]] <- rep(triangles[[i]], each = 4)
n.tri <- 4 * n.tri
}
triangles$N <- interpolateVertexNormals(VN, ib)
triangles$v1 <- t(vb[,ib[1,]])
triangles$v2 <- t(vb[,ib[2,]])
triangles$v3 <- t(vb[,ib[3,]])
triangles
}
triangleNormalsPhongEX <- function(triangles, reps = 1) {
N <- triangleNormals(triangles)
ve <- t2ve(triangles)
vt <- vertexTriangles(ve)
VN <- vertexNormals(vt, N)
vb <- ve$vb
ib <- ve$ib
n.tri <- nrow(N)
color <- rep(triangles$color, length = n.tri)
color2 <- rep(triangles$color2, length = n.tri)
col.mesh <- rep(triangles$col.mesh, length = n.tri)
color2 <- ifelse(is.na(color2), color, color2)
col.mesh <- ifelse(is.na(col.mesh), color, col.mesh)
VC <- vertexColors(vt, color)
VC2 <- vertexColors(vt, color2)
VCm <- vertexColors(vt, col.mesh)
while (reps > 0) {
reps <- reps - 1
n.ver <- nrow(VN)
edges <- triangleEdges(vb, ib)
VC <- rbind(VC, (VC[edges[1,],] + VC[edges[2,],]) / 2)
VC2 <- rbind(VC2, (VC2[edges[1,],] + VC2[edges[2,],]) / 2)
VCm <- rbind(VCm, (VCm[edges[1,],] + VCm[edges[2,],]) / 2)
mt <- triangleMidTriangles(vb, ib, VN)
vb <- cbind(vb, mt$vb)
VN <- rbind(VN, mt$VN)
mtib <- mt$ib + n.ver
ib <- matrix(rbind(ib[1,], mtib[1,],mtib[3,],
mtib[1,], ib[2,],mtib[2,],
mtib[2,], ib[3,], mtib[3,],
mtib), nrow = 3)
for (i in seq(along = triangles))
if (length(triangles[[i]]) == n.tri)
triangles[[i]] <- rep(triangles[[i]], each = 4)
n.tri <- 4 * n.tri
}
triangles$color <- interpolateVertexColors(VC, ib)
triangles$color2 <- interpolateVertexColors(VC2, ib)
triangles$color.mesh <- interpolateVertexColors(VCm, ib)
triangles$N <- interpolateVertexNormals(VN, ib)
triangles$v1 <- t(vb[,ib[1,]])
triangles$v2 <- t(vb[,ib[2,]])
triangles$v3 <- t(vb[,ib[3,]])
triangles
}
lightTriangles <- function(triangles, lighting, light) {
view <- c(0, 0, 1)
normals <- triangleNormals(triangles)
smooth <- if (is.null(triangles$smooth)) 0 else triangles$smooth
if (smooth == 0)
normals <- triangleNormals(triangles)
else if (smooth == 1)
normals <- triangleNormalsPhong(triangles)
else {
triangles <- triangleNormalsPhongEX(triangles, reps = smooth - 1)
normals <- triangles$N
}
n.tri <- nrow(normals)
color <- rep(triangles$color, length = n.tri)
color2 <- rep(triangles$color2, length = n.tri)
color2 <- ifelse(is.na(color2), color, color2)
alpha <- rep(triangles$alpha, length = n.tri)
mat <- triangles$material
triangles$col.light <- lighting(normals, view, light, color, color2, alpha, mat)
triangles
}
lightScene <- function(scene, lighting, light) {
if (is.Triangles3D(scene)) lightTriangles(scene, lighting, light)
else lapply(scene, lightTriangles, lighting, light)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.